Bug Summary

File:liboctave/numeric/fCmplxCHOL.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 "fMatrix.h"
31#include "fRowVector.h"
32#include "fCmplxCHOL.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 (cpotrf, CPOTRF)cpotrf_ (F77_CONST_CHAR_ARG_DECLconst char *,
45 const octave_idx_type&, FloatComplex*,
46 const octave_idx_type&, octave_idx_type&
47 F77_CHAR_ARG_LEN_DECL, long);
48 F77_RET_Tint
49 F77_FUNC (cpotri, CPOTRI)cpotri_ (F77_CONST_CHAR_ARG_DECLconst char *,
50 const octave_idx_type&, FloatComplex*,
51 const octave_idx_type&, octave_idx_type&
52 F77_CHAR_ARG_LEN_DECL, long);
53
54 F77_RET_Tint
55 F77_FUNC (cpocon, CPOCON)cpocon_ (F77_CONST_CHAR_ARG_DECLconst char *,
56 const octave_idx_type&, FloatComplex*,
57 const octave_idx_type&, const float&,
58 float&, FloatComplex*, float*, octave_idx_type&
59 F77_CHAR_ARG_LEN_DECL, long);
60#ifdef HAVE_QRUPDATE1
61
62 F77_RET_Tint
63 F77_FUNC (cch1up, CCH1UP)cch1up_ (const octave_idx_type&, FloatComplex*,
64 const octave_idx_type&, FloatComplex*, float*);
65
66 F77_RET_Tint
67 F77_FUNC (cch1dn, CCH1DN)cch1dn_ (const octave_idx_type&, FloatComplex*,
68 const octave_idx_type&, FloatComplex*,
69 float*, octave_idx_type&);
70
71 F77_RET_Tint
72 F77_FUNC (cchinx, CCHINX)cchinx_ (const octave_idx_type&, FloatComplex*,
73 const octave_idx_type&, const octave_idx_type&,
74 FloatComplex*, float*, octave_idx_type&);
75
76 F77_RET_Tint
77 F77_FUNC (cchdex, CCHDEX)cchdex_ (const octave_idx_type&, FloatComplex*,
78 const octave_idx_type&, const octave_idx_type&,
79 float*);
80
81 F77_RET_Tint
82 F77_FUNC (cchshx, CCHSHX)cchshx_ (const octave_idx_type&, FloatComplex*,
83 const octave_idx_type&, const octave_idx_type&,
84 const octave_idx_type&, FloatComplex*, float*);
85#endif
86}
87
88octave_idx_type
89FloatComplexCHOL::init (const FloatComplexMatrix& 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 ("FloatComplexCHOL 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.0f;
111 }
112 FloatComplex *h = chol_mat.fortran_vec ();
113
114 // Calculate the norm of the matrix, for later use.
115 float 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 (cpotrf, CPOTRF, (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"
, "cpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cpotrf_ ("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"
, "cpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cpotrf_ ("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 cpocon_info = 0;
128
129 // Now calculate the condition number for non-singular matrix.
130 Array<FloatComplex> z (dim_vector (2*n, 1));
131 FloatComplex *pz = z.fortran_vec ();
132 Array<float> rz (dim_vector (n, 1));
133 float *prz = rz.fortran_vec ();
134 F77_XFCN (cpocon, CPOCON, (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"
, "cpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cpocon_ ("U", n, h, n, anorm, xrcond, pz, prz, cpocon_info
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
135 n, anorm, xrcond, pz, prz, cpocon_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"
, "cpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cpocon_ ("U", n, h, n, anorm, xrcond, pz, prz, cpocon_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"
, "cpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cpocon_ ("U", n, h, n, anorm, xrcond, pz, prz, cpocon_info
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
137
138 if (cpocon_info != 0)
139 info = -1;
140 }
141
142 return info;
143}
144
145static FloatComplexMatrix
146chol2inv_internal (const FloatComplexMatrix& r)
147{
148 FloatComplexMatrix 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 FloatComplexMatrix tmp = r;
159
160 F77_XFCN (cpotri, CPOTRI, (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"
, "cpotri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cpotri_ ("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"
, "cpotri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cpotri_ ("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"
, "cpotri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cpotri_ ("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.
181FloatComplexMatrix
182FloatComplexCHOL::inverse (void) const
183{
184 return chol2inv_internal (chol_mat);
185}
186
187void
188FloatComplexCHOL::set (const FloatComplexMatrix& 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
199FloatComplexCHOL::update (const FloatComplexColumnVector& u)
200{
201 octave_idx_type n = chol_mat.rows ();
202
203 if (u.length () == n)
204 {
205 FloatComplexColumnVector utmp = u;
206
207 OCTAVE_LOCAL_BUFFER (float, rw, n)octave_local_buffer<float> _buffer_rw (n); float *rw = _buffer_rw;
208
209 F77_XFCN (cch1up, CCH1UP, (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"
, "cch1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cch1up_ (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"
, "cch1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cch1up_ (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
217FloatComplexCHOL::downdate (const FloatComplexColumnVector& 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 FloatComplexColumnVector utmp = u;
226
227 OCTAVE_LOCAL_BUFFER (float, rw, n)octave_local_buffer<float> _buffer_rw (n); float *rw = _buffer_rw;
228
229 F77_XFCN (cch1dn, CCH1DN, (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"
, "cch1dn_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cch1dn_ (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"
, "cch1dn_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cch1dn_ (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
239FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u,
240 octave_idx_type j)
241{
242 octave_idx_type info = -1;
243
244 octave_idx_type n = chol_mat.rows ();
245
246 if (u.length () != n + 1)
247 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
248 else if (j < 0 || j > n)
249 (*current_liboctave_error_handler) ("cholinsert: index out of range");
250 else
251 {
252 FloatComplexColumnVector utmp = u;
253
254 OCTAVE_LOCAL_BUFFER (float, rw, n)octave_local_buffer<float> _buffer_rw (n); float *rw = _buffer_rw;
255
256 chol_mat.resize (n+1, n+1);
257
258 F77_XFCN (cchinx, CCHINX, (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"
, "cchinx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cchinx_ (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 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"
, "cchinx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cchinx_ (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)
;
260 }
261
262 return info;
263}
264
265void
266FloatComplexCHOL::delete_sym (octave_idx_type j)
267{
268 octave_idx_type n = chol_mat.rows ();
269
270 if (j < 0 || j > n-1)
271 (*current_liboctave_error_handler) ("choldelete: index out of range");
272 else
273 {
274 OCTAVE_LOCAL_BUFFER (float, rw, n)octave_local_buffer<float> _buffer_rw (n); float *rw = _buffer_rw;
275
276 F77_XFCN (cchdex, CCHDEX, (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"
, "cchdex_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cchdex_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, rw); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
277 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"
, "cchdex_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cchdex_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, rw); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
278
279 chol_mat.resize (n-1, n-1);
280 }
281}
282
283void
284FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
285{
286 octave_idx_type n = chol_mat.rows ();
287
288 if (i < 0 || i > n-1 || j < 0 || j > n-1)
289 (*current_liboctave_error_handler) ("cholshift: index out of range");
290 else
291 {
292 OCTAVE_LOCAL_BUFFER (FloatComplex, w, n)octave_local_buffer<FloatComplex> _buffer_w (n); FloatComplex
*w = _buffer_w
;
293 OCTAVE_LOCAL_BUFFER (float, rw, n)octave_local_buffer<float> _buffer_rw (n); float *rw = _buffer_rw;
294
295 F77_XFCN (cchshx, CCHSHX, (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"
, "cchshx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cchshx_ (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 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"
, "cchshx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; cchshx_ (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)
;
297 }
298}
299
300#else
301
302void
303FloatComplexCHOL::update (const FloatComplexColumnVector& u)
304{
305 warn_qrupdate_once ();
306
307 octave_idx_type n = chol_mat.rows ();
308
309 if (u.length () == n)
310 {
311 init (chol_mat.hermitian () * chol_mat
312 + FloatComplexMatrix (u) * FloatComplexMatrix (u).hermitian (),
313 false);
314 }
315 else
316 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
317}
318
319static bool
320singular (const FloatComplexMatrix& a)
321{
322 for (octave_idx_type i = 0; i < a.rows (); i++)
323 if (a(i,i) == 0.0f) return true;
324 return false;
325}
326
327octave_idx_type
328FloatComplexCHOL::downdate (const FloatComplexColumnVector& u)
329{
330 warn_qrupdate_once ();
331
332 octave_idx_type info = -1;
333
334 octave_idx_type n = chol_mat.rows ();
335
336 if (u.length () == n)
337 {
338 if (singular (chol_mat))
339 info = 2;
340 else
341 {
342 info = init (chol_mat.hermitian () * chol_mat
343 - FloatComplexMatrix (u)
344 * FloatComplexMatrix (u).hermitian (),
345 false);
346 if (info) info = 1;
347 }
348 }
349 else
350 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
351
352 return info;
353}
354
355octave_idx_type
356FloatComplexCHOL::insert_sym (const FloatComplexColumnVector& u,
357 octave_idx_type j)
358{
359 warn_qrupdate_once ();
360
361 octave_idx_type info = -1;
362
363 octave_idx_type n = chol_mat.rows ();
364
365 if (u.length () != n + 1)
366 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
367 else if (j < 0 || j > n)
368 (*current_liboctave_error_handler) ("cholinsert: index out of range");
369 else
370 {
371 if (singular (chol_mat))
372 info = 2;
373 else if (u(j).imag () != 0.0)
374 info = 3;
375 else
376 {
377 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
378 FloatComplexMatrix a1 (n+1, n+1);
379 for (octave_idx_type k = 0; k < n+1; k++)
380 for (octave_idx_type l = 0; l < n+1; l++)
381 {
382 if (l == j)
383 a1(k, l) = u(k);
384 else if (k == j)
385 a1(k, l) = std::conj (u(l));
386 else
387 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
388 }
389 info = init (a1, false);
390 if (info) info = 1;
391 }
392 }
393
394 return info;
395}
396
397void
398FloatComplexCHOL::delete_sym (octave_idx_type j)
399{
400 warn_qrupdate_once ();
401
402 octave_idx_type n = chol_mat.rows ();
403
404 if (j < 0 || j > n-1)
405 (*current_liboctave_error_handler) ("choldelete: index out of range");
406 else
407 {
408 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
409 a.delete_elements (1, idx_vector (j));
410 a.delete_elements (0, idx_vector (j));
411 init (a, false);
412 }
413}
414
415void
416FloatComplexCHOL::shift_sym (octave_idx_type i, octave_idx_type j)
417{
418 warn_qrupdate_once ();
419
420 octave_idx_type n = chol_mat.rows ();
421
422 if (i < 0 || i > n-1 || j < 0 || j > n-1)
423 (*current_liboctave_error_handler) ("cholshift: index out of range");
424 else
425 {
426 FloatComplexMatrix a = chol_mat.hermitian () * chol_mat;
427 Array<octave_idx_type> p (dim_vector (n, 1));
428 for (octave_idx_type k = 0; k < n; k++) p(k) = k;
429 if (i < j)
430 {
431 for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
432 p(j) = i;
433 }
434 else if (j < i)
435 {
436 p(j) = i;
437 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
438 }
439
440 init (a.index (idx_vector (p), idx_vector (p)), false);
441 }
442}
443
444#endif
445
446FloatComplexMatrix
447chol2inv (const FloatComplexMatrix& r)
448{
449 return chol2inv_internal (r);
450}