Bug Summary

File:liboctave/numeric/CmplxCHOL.cc
Location:line 123, column 12
Description:The left operand of '>' is a garbage value

Annotated Source Code

1/*
2
3Copyright (C) 1994-2013 John W. Eaton
4Copyright (C) 2008-2009 Jaroslav Hajek
5
6This file is part of Octave.
7
8Octave is free software; you can redistribute it and/or modify it
9under the terms of the GNU General Public License as published by the
10Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13Octave is distributed in the hope that it will be useful, but WITHOUT
14ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16for more details.
17
18You should have received a copy of the GNU General Public License
19along with Octave; see the file COPYING. If not, see
20<http://www.gnu.org/licenses/>.
21
22*/
23
24#ifdef HAVE_CONFIG_H1
25#include <config.h>
26#endif
27
28#include <vector>
29
30#include "dMatrix.h"
31#include "dRowVector.h"
32#include "CmplxCHOL.h"
33#include "f77-fcn.h"
34#include "lo-error.h"
35#include "oct-locbuf.h"
36#include "oct-norm.h"
37#ifndef HAVE_QRUPDATE1
38#include "dbleQR.h"
39#endif
40
41extern "C"
42{
43 F77_RET_Tint
44 F77_FUNC (zpotrf, ZPOTRF)zpotrf_ (F77_CONST_CHAR_ARG_DECLconst char *,
45 const octave_idx_type&, Complex*,
46 const octave_idx_type&, octave_idx_type&
47 F77_CHAR_ARG_LEN_DECL, long);
48 F77_RET_Tint
49 F77_FUNC (zpotri, ZPOTRI)zpotri_ (F77_CONST_CHAR_ARG_DECLconst char *,
50 const octave_idx_type&, Complex*,
51 const octave_idx_type&, octave_idx_type&
52 F77_CHAR_ARG_LEN_DECL, long);
53
54 F77_RET_Tint
55 F77_FUNC (zpocon, ZPOCON)zpocon_ (F77_CONST_CHAR_ARG_DECLconst char *,
56 const octave_idx_type&, Complex*,
57 const octave_idx_type&, const double&,
58 double&, Complex*, double*, octave_idx_type&
59 F77_CHAR_ARG_LEN_DECL, long);
60#ifdef HAVE_QRUPDATE1
61
62 F77_RET_Tint
63 F77_FUNC (zch1up, ZCH1UP)zch1up_ (const octave_idx_type&, Complex*,
64 const octave_idx_type&, Complex*, double*);
65
66 F77_RET_Tint
67 F77_FUNC (zch1dn, ZCH1DN)zch1dn_ (const octave_idx_type&, Complex*,
68 const octave_idx_type&, Complex*, double*,
69 octave_idx_type&);
70
71 F77_RET_Tint
72 F77_FUNC (zchinx, ZCHINX)zchinx_ (const octave_idx_type&, Complex*,
73 const octave_idx_type&, const octave_idx_type&,
74 Complex*, double*, octave_idx_type&);
75
76 F77_RET_Tint
77 F77_FUNC (zchdex, ZCHDEX)zchdex_ (const octave_idx_type&, Complex*,
78 const octave_idx_type&, const octave_idx_type&,
79 double*);
80
81 F77_RET_Tint
82 F77_FUNC (zchshx, ZCHSHX)zchshx_ (const octave_idx_type&, Complex*,
83 const octave_idx_type&, const octave_idx_type&,
84 const octave_idx_type&, Complex*, double*);
85#endif
86}
87
88octave_idx_type
89ComplexCHOL::init (const ComplexMatrix& a, bool calc_cond)
90{
91 octave_idx_type a_nr = a.rows ();
92 octave_idx_type a_nc = a.cols ();
93
94 if (a_nr != a_nc)
1
Taking false branch
95 {
96 (*current_liboctave_error_handler)
97 ("ComplexCHOL requires square matrix");
98 return -1;
99 }
100
101 octave_idx_type n = a_nc;
102 octave_idx_type info;
2
'info' declared without an initial value
103
104 chol_mat.clear (n, n);
105 for (octave_idx_type j = 0; j < n; j++)
3
Assuming 'j' is >= 'n'
4
Loop condition is false. Execution continues on line 112
106 {
107 for (octave_idx_type i = 0; i <= j; i++)
108 chol_mat.xelem (i, j) = a(i, j);
109 for (octave_idx_type i = j+1; i < n; i++)
110 chol_mat.xelem (i, j) = 0.0;
111 }
112 Complex *h = chol_mat.fortran_vec ();
113
114 // Calculate the norm of the matrix, for later use.
115 double anorm = 0;
116 if (calc_cond)
5
Assuming 'calc_cond' is 0
6
Taking false branch
117 anorm = xnorm (a, 1);
118
119 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ ("U", n, h, n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
120 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotrf_ ("U", n, h, n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
121
122 xrcond = 0.0;
123 if (info > 0)
7
The left operand of '>' is a garbage value
124 chol_mat.resize (info - 1, info - 1);
125 else if (calc_cond)
126 {
127 octave_idx_type zpocon_info = 0;
128
129 // Now calculate the condition number for non-singular matrix.
130 Array<Complex> z (dim_vector (2*n, 1));
131 Complex *pz = z.fortran_vec ();
132 Array<double> rz (dim_vector (n, 1));
133 double *prz = rz.fortran_vec ();
134 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ ("U", n, h, n, anorm, xrcond, pz, prz, zpocon_info
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
135 n, anorm, xrcond, pz, prz, zpocon_infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ ("U", n, h, n, anorm, xrcond, pz, prz, zpocon_info
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
136 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpocon_ ("U", n, h, n, anorm, xrcond, pz, prz, zpocon_info
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
137
138 if (zpocon_info != 0)
139 info = -1;
140 }
141
142 return info;
143}
144
145static ComplexMatrix
146chol2inv_internal (const ComplexMatrix& r)
147{
148 ComplexMatrix retval;
149
150 octave_idx_type r_nr = r.rows ();
151 octave_idx_type r_nc = r.cols ();
152
153 if (r_nr == r_nc)
154 {
155 octave_idx_type n = r_nc;
156 octave_idx_type info;
157
158 ComplexMatrix tmp = r;
159
160 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotri_ ("U", n, tmp.fortran_vec (), n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
161 tmp.fortran_vec (), n, infodo { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotri_ ("U", n, tmp.fortran_vec (), n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
162 F77_CHAR_ARG_LEN (1)))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zpotri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zpotri_ ("U", n, tmp.fortran_vec (), n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
163
164 // If someone thinks of a more graceful way of doing this (or
165 // faster for that matter :-)), please let me know!
166
167 if (n > 1)
168 for (octave_idx_type j = 0; j < r_nc; j++)
169 for (octave_idx_type i = j+1; i < r_nr; i++)
170 tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
171
172 retval = tmp;
173 }
174 else
175 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
176
177 return retval;
178}
179
180// Compute the inverse of a matrix using the Cholesky factorization.
181ComplexMatrix
182ComplexCHOL::inverse (void) const
183{
184 return chol2inv_internal (chol_mat);
185}
186
187void
188ComplexCHOL::set (const ComplexMatrix& R)
189{
190 if (R.is_square ())
191 chol_mat = R;
192 else
193 (*current_liboctave_error_handler) ("CHOL requires square matrix");
194}
195
196#ifdef HAVE_QRUPDATE1
197
198void
199ComplexCHOL::update (const ComplexColumnVector& u)
200{
201 octave_idx_type n = chol_mat.rows ();
202
203 if (u.length () == n)
204 {
205 ComplexColumnVector utmp = u;
206
207 OCTAVE_LOCAL_BUFFER (double, rw, n)octave_local_buffer<double> _buffer_rw (n); double *rw =
_buffer_rw
;
208
209 F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (), chol_mat.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zch1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zch1up_ (n, chol_mat.fortran_vec (), chol_mat.rows (), utmp
.fortran_vec (), rw); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
210 utmp.fortran_vec (), rw))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zch1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zch1up_ (n, chol_mat.fortran_vec (), chol_mat.rows (), utmp
.fortran_vec (), rw); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
211 }
212 else
213 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
214}
215
216octave_idx_type
217ComplexCHOL::downdate (const ComplexColumnVector& u)
218{
219 octave_idx_type info = -1;
220
221 octave_idx_type n = chol_mat.rows ();
222
223 if (u.length () == n)
224 {
225 ComplexColumnVector utmp = u;
226
227 OCTAVE_LOCAL_BUFFER (double, rw, n)octave_local_buffer<double> _buffer_rw (n); double *rw =
_buffer_rw
;
228
229 F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (), chol_mat.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zch1dn_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zch1dn_ (n, chol_mat.fortran_vec (), chol_mat.rows (), utmp
.fortran_vec (), rw, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
230 utmp.fortran_vec (), rw, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zch1dn_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zch1dn_ (n, chol_mat.fortran_vec (), chol_mat.rows (), utmp
.fortran_vec (), rw, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
231 }
232 else
233 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
234
235 return info;
236}
237
238octave_idx_type
239ComplexCHOL::insert_sym (const ComplexColumnVector& u, octave_idx_type j)
240{
241 octave_idx_type info = -1;
242
243 octave_idx_type n = chol_mat.rows ();
244
245 if (u.length () != n + 1)
246 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
247 else if (j < 0 || j > n)
248 (*current_liboctave_error_handler) ("cholinsert: index out of range");
249 else
250 {
251 ComplexColumnVector utmp = u;
252
253 OCTAVE_LOCAL_BUFFER (double, rw, n)octave_local_buffer<double> _buffer_rw (n); double *rw =
_buffer_rw
;
254
255 chol_mat.resize (n+1, n+1);
256
257 F77_XFCN (zchinx, ZCHINX, (n, chol_mat.fortran_vec (), chol_mat.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zchinx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zchinx_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, utmp.fortran_vec (), rw, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
258 j + 1, utmp.fortran_vec (), rw, info))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zchinx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zchinx_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, utmp.fortran_vec (), rw, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
259 }
260
261 return info;
262}
263
264void
265ComplexCHOL::delete_sym (octave_idx_type j)
266{
267 octave_idx_type n = chol_mat.rows ();
268
269 if (j < 0 || j > n-1)
270 (*current_liboctave_error_handler) ("choldelete: index out of range");
271 else
272 {
273 OCTAVE_LOCAL_BUFFER (double, rw, n)octave_local_buffer<double> _buffer_rw (n); double *rw =
_buffer_rw
;
274
275 F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zchdex_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zchdex_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, rw); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
276 j + 1, rw))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zchdex_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zchdex_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, rw); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
277
278 chol_mat.resize (n-1, n-1);
279 }
280}
281
282void
283ComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
284{
285 octave_idx_type n = chol_mat.rows ();
286
287 if (i < 0 || i > n-1 || j < 0 || j > n-1)
288 (*current_liboctave_error_handler) ("cholshift: index out of range");
289 else
290 {
291 OCTAVE_LOCAL_BUFFER (Complex, w, n)octave_local_buffer<Complex> _buffer_w (n); Complex *w =
_buffer_w
;
292 OCTAVE_LOCAL_BUFFER (double, rw, n)octave_local_buffer<double> _buffer_rw (n); double *rw =
_buffer_rw
;
293
294 F77_XFCN (zchshx, ZCHSHX, (n, chol_mat.fortran_vec (), chol_mat.rows (),do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zchshx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zchshx_ (n, chol_mat.fortran_vec (), chol_mat.rows (), i +
1, j + 1, w, rw); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
295 i + 1, j + 1, w, rw))do { octave_jmp_buf saved_context; sig_atomic_t saved_octave_interrupt_immediately
= octave_interrupt_immediately; f77_exception_encountered = 0
; octave_save_current_context (saved_context); if (_setjmp (current_context
)) { octave_interrupt_immediately = saved_octave_interrupt_immediately
; octave_restore_current_context (saved_context); if (f77_exception_encountered
) (*current_liboctave_error_handler) ("exception encountered in Fortran subroutine %s"
, "zchshx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; zchshx_ (n, chol_mat.fortran_vec (), chol_mat.rows (), i +
1, j + 1, w, rw); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
296 }
297}
298
299#else
300
301void
302ComplexCHOL::update (const ComplexColumnVector& u)
303{
304 warn_qrupdate_once ();
305
306 octave_idx_type n = chol_mat.rows ();
307
308 if (u.length () == n)
309 {
310 init (chol_mat.hermitian () * chol_mat
311 + ComplexMatrix (u) * ComplexMatrix (u).hermitian (), false);
312 }
313 else
314 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
315}
316
317static bool
318singular (const ComplexMatrix& a)
319{
320 for (octave_idx_type i = 0; i < a.rows (); i++)
321 if (a(i,i) == 0.0) return true;
322 return false;
323}
324
325octave_idx_type
326ComplexCHOL::downdate (const ComplexColumnVector& u)
327{
328 warn_qrupdate_once ();
329
330 octave_idx_type info = -1;
331
332 octave_idx_type n = chol_mat.rows ();
333
334 if (u.length () == n)
335 {
336 if (singular (chol_mat))
337 info = 2;
338 else
339 {
340 info = init (chol_mat.hermitian () * chol_mat
341 - ComplexMatrix (u) * ComplexMatrix (u).hermitian (),
342 false);
343 if (info) info = 1;
344 }
345 }
346 else
347 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
348
349 return info;
350}
351
352octave_idx_type
353ComplexCHOL::insert_sym (const ComplexColumnVector& u, octave_idx_type j)
354{
355 warn_qrupdate_once ();
356
357 octave_idx_type info = -1;
358
359 octave_idx_type n = chol_mat.rows ();
360
361 if (u.length () != n + 1)
362 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
363 else if (j < 0 || j > n)
364 (*current_liboctave_error_handler) ("cholinsert: index out of range");
365 else
366 {
367 if (singular (chol_mat))
368 info = 2;
369 else if (u(j).imag () != 0.0)
370 info = 3;
371 else
372 {
373 ComplexMatrix a = chol_mat.hermitian () * chol_mat;
374 ComplexMatrix a1 (n+1, n+1);
375 for (octave_idx_type k = 0; k < n+1; k++)
376 for (octave_idx_type l = 0; l < n+1; l++)
377 {
378 if (l == j)
379 a1(k, l) = u(k);
380 else if (k == j)
381 a1(k, l) = std::conj (u(l));
382 else
383 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
384 }
385 info = init (a1, false);
386 if (info) info = 1;
387 }
388 }
389
390 return info;
391}
392
393void
394ComplexCHOL::delete_sym (octave_idx_type j)
395{
396 warn_qrupdate_once ();
397
398 octave_idx_type n = chol_mat.rows ();
399
400 if (j < 0 || j > n-1)
401 (*current_liboctave_error_handler) ("choldelete: index out of range");
402 else
403 {
404 ComplexMatrix a = chol_mat.hermitian () * chol_mat;
405 a.delete_elements (1, idx_vector (j));
406 a.delete_elements (0, idx_vector (j));
407 init (a, false);
408 }
409}
410
411void
412ComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
413{
414 warn_qrupdate_once ();
415
416 octave_idx_type n = chol_mat.rows ();
417
418 if (i < 0 || i > n-1 || j < 0 || j > n-1)
419 (*current_liboctave_error_handler) ("cholshift: index out of range");
420 else
421 {
422 ComplexMatrix a = chol_mat.hermitian () * chol_mat;
423 Array<octave_idx_type> p (dim_vector (n, 1));
424 for (octave_idx_type k = 0; k < n; k++) p(k) = k;
425 if (i < j)
426 {
427 for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
428 p(j) = i;
429 }
430 else if (j < i)
431 {
432 p(j) = i;
433 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
434 }
435
436 init (a.index (idx_vector (p), idx_vector (p)), false);
437 }
438}
439
440#endif
441
442ComplexMatrix
443chol2inv (const ComplexMatrix& r)
444{
445 return chol2inv_internal (r);
446}