Bug Summary

File:liboctave/numeric/dbleCHOL.cc
Location:line 124, 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 "dRowVector.h"
31#include "dbleCHOL.h"
32#include "f77-fcn.h"
33#include "lo-error.h"
34#include "oct-locbuf.h"
35#include "oct-norm.h"
36#ifndef HAVE_QRUPDATE1
37#include "dbleQR.h"
38#endif
39
40extern "C"
41{
42 F77_RET_Tint
43 F77_FUNC (dpotrf, DPOTRF)dpotrf_ (F77_CONST_CHAR_ARG_DECLconst char *,
44 const octave_idx_type&, double*,
45 const octave_idx_type&, octave_idx_type&
46 F77_CHAR_ARG_LEN_DECL, long);
47
48 F77_RET_Tint
49 F77_FUNC (dpotri, DPOTRI)dpotri_ (F77_CONST_CHAR_ARG_DECLconst char *,
50 const octave_idx_type&, double*,
51 const octave_idx_type&, octave_idx_type&
52 F77_CHAR_ARG_LEN_DECL, long);
53
54 F77_RET_Tint
55 F77_FUNC (dpocon, DPOCON)dpocon_ (F77_CONST_CHAR_ARG_DECLconst char *,
56 const octave_idx_type&, double*,
57 const octave_idx_type&, const double&,
58 double&, double*, octave_idx_type*,
59 octave_idx_type&
60 F77_CHAR_ARG_LEN_DECL, long);
61#ifdef HAVE_QRUPDATE1
62
63 F77_RET_Tint
64 F77_FUNC (dch1up, DCH1UP)dch1up_ (const octave_idx_type&, double*,
65 const octave_idx_type&, double*, double*);
66
67 F77_RET_Tint
68 F77_FUNC (dch1dn, DCH1DN)dch1dn_ (const octave_idx_type&, double*,
69 const octave_idx_type&, double*, double*,
70 octave_idx_type&);
71
72 F77_RET_Tint
73 F77_FUNC (dchinx, DCHINX)dchinx_ (const octave_idx_type&, double*,
74 const octave_idx_type&, const octave_idx_type&,
75 double*, double*, octave_idx_type&);
76
77 F77_RET_Tint
78 F77_FUNC (dchdex, DCHDEX)dchdex_ (const octave_idx_type&, double*,
79 const octave_idx_type&, const octave_idx_type&,
80 double*);
81
82 F77_RET_Tint
83 F77_FUNC (dchshx, DCHSHX)dchshx_ (const octave_idx_type&, double*,
84 const octave_idx_type&, const octave_idx_type&,
85 const octave_idx_type&, double*);
86#endif
87}
88
89octave_idx_type
90CHOL::init (const Matrix& a, bool calc_cond)
91{
92 octave_idx_type a_nr = a.rows ();
93 octave_idx_type a_nc = a.cols ();
94
95 if (a_nr != a_nc)
1
Taking false branch
96 {
97 (*current_liboctave_error_handler) ("CHOL 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 double *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 (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 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"
, "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dpotrf_ ("U", n, h, n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
120 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"
, "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dpotrf_ ("U", n, h, n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
121 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"
, "dpotrf_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dpotrf_ ("U", n, h, n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
122
123 xrcond = 0.0;
124 if (info > 0)
7
The left operand of '>' is a garbage value
125 chol_mat.resize (info - 1, info - 1);
126 else if (calc_cond)
127 {
128 octave_idx_type dpocon_info = 0;
129
130 // Now calculate the condition number for non-singular matrix.
131 Array<double> z (dim_vector (3*n, 1));
132 double *pz = z.fortran_vec ();
133 Array<octave_idx_type> iz (dim_vector (n, 1));
134 octave_idx_type *piz = iz.fortran_vec ();
135 F77_XFCN (dpocon, DPOCON, (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"
, "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dpocon_ ("U", n, h, n, anorm, xrcond, pz, piz, dpocon_info
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
136 n, anorm, xrcond, pz, piz, dpocon_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"
, "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dpocon_ ("U", n, h, n, anorm, xrcond, pz, piz, dpocon_info
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
137 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"
, "dpocon_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dpocon_ ("U", n, h, n, anorm, xrcond, pz, piz, dpocon_info
, 1); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
138
139 if (dpocon_info != 0)
140 info = -1;
141 }
142
143 return info;
144}
145
146static Matrix
147chol2inv_internal (const Matrix& r)
148{
149 Matrix retval;
150
151 octave_idx_type r_nr = r.rows ();
152 octave_idx_type r_nc = r.cols ();
153
154 if (r_nr == r_nc)
155 {
156 octave_idx_type n = r_nc;
157 octave_idx_type info = 0;
158
159 Matrix tmp = r;
160 double *v = tmp.fortran_vec ();
161
162 if (info == 0)
163 {
164 F77_XFCN (dpotri, DPOTRI, (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"
, "dpotri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dpotri_ ("U", n, v, n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
165 v, 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"
, "dpotri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dpotri_ ("U", n, v, n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
166 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"
, "dpotri_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dpotri_ ("U", n, v, n, info , 1); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
167
168 // If someone thinks of a more graceful way of doing this (or
169 // faster for that matter :-)), please let me know!
170
171 if (n > 1)
172 for (octave_idx_type j = 0; j < r_nc; j++)
173 for (octave_idx_type i = j+1; i < r_nr; i++)
174 tmp.xelem (i, j) = tmp.xelem (j, i);
175
176 retval = tmp;
177 }
178 }
179 else
180 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
181
182 return retval;
183}
184
185// Compute the inverse of a matrix using the Cholesky factorization.
186Matrix
187CHOL::inverse (void) const
188{
189 return chol2inv_internal (chol_mat);
190}
191
192void
193CHOL::set (const Matrix& R)
194{
195 if (R.is_square ())
196 chol_mat = R;
197 else
198 (*current_liboctave_error_handler) ("CHOL requires square matrix");
199}
200
201#ifdef HAVE_QRUPDATE1
202
203void
204CHOL::update (const ColumnVector& u)
205{
206 octave_idx_type n = chol_mat.rows ();
207
208 if (u.length () == n)
209 {
210 ColumnVector utmp = u;
211
212 OCTAVE_LOCAL_BUFFER (double, w, n)octave_local_buffer<double> _buffer_w (n); double *w = _buffer_w;
213
214 F77_XFCN (dch1up, DCH1UP, (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"
, "dch1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dch1up_ (n, chol_mat.fortran_vec (), chol_mat.rows (), utmp
.fortran_vec (), w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
215 utmp.fortran_vec (), w))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"
, "dch1up_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dch1up_ (n, chol_mat.fortran_vec (), chol_mat.rows (), utmp
.fortran_vec (), w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
216 }
217 else
218 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
219}
220
221octave_idx_type
222CHOL::downdate (const ColumnVector& u)
223{
224 octave_idx_type info = -1;
225
226 octave_idx_type n = chol_mat.rows ();
227
228 if (u.length () == n)
229 {
230 ColumnVector utmp = u;
231
232 OCTAVE_LOCAL_BUFFER (double, w, n)octave_local_buffer<double> _buffer_w (n); double *w = _buffer_w;
233
234 F77_XFCN (dch1dn, DCH1DN, (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"
, "dch1dn_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dch1dn_ (n, chol_mat.fortran_vec (), chol_mat.rows (), utmp
.fortran_vec (), w, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
235 utmp.fortran_vec (), w, 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"
, "dch1dn_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dch1dn_ (n, chol_mat.fortran_vec (), chol_mat.rows (), utmp
.fortran_vec (), w, info); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
236 }
237 else
238 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
239
240 return info;
241}
242
243octave_idx_type
244CHOL::insert_sym (const ColumnVector& u, octave_idx_type j)
245{
246 octave_idx_type info = -1;
247
248 octave_idx_type n = chol_mat.rows ();
249
250 if (u.length () != n + 1)
251 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
252 else if (j < 0 || j > n)
253 (*current_liboctave_error_handler) ("cholinsert: index out of range");
254 else
255 {
256 ColumnVector utmp = u;
257
258 OCTAVE_LOCAL_BUFFER (double, w, n)octave_local_buffer<double> _buffer_w (n); double *w = _buffer_w;
259
260 chol_mat.resize (n+1, n+1);
261
262 F77_XFCN (dchinx, DCHINX, (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"
, "dchinx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dchinx_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, utmp.fortran_vec (), w, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
263 j + 1, utmp.fortran_vec (), w, 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"
, "dchinx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dchinx_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, utmp.fortran_vec (), w, info); octave_interrupt_immediately
--; octave_restore_current_context (saved_context); } } while
(0)
;
264 }
265
266 return info;
267}
268
269void
270CHOL::delete_sym (octave_idx_type j)
271{
272 octave_idx_type n = chol_mat.rows ();
273
274 if (j < 0 || j > n-1)
275 (*current_liboctave_error_handler) ("choldelete: index out of range");
276 else
277 {
278 OCTAVE_LOCAL_BUFFER (double, w, n)octave_local_buffer<double> _buffer_w (n); double *w = _buffer_w;
279
280 F77_XFCN (dchdex, DCHDEX, (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"
, "dchdex_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dchdex_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
281 j + 1, w))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"
, "dchdex_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dchdex_ (n, chol_mat.fortran_vec (), chol_mat.rows (), j +
1, w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
282
283 chol_mat.resize (n-1, n-1);
284 }
285}
286
287void
288CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
289{
290 octave_idx_type n = chol_mat.rows ();
291
292 if (i < 0 || i > n-1 || j < 0 || j > n-1)
293 (*current_liboctave_error_handler) ("cholshift: index out of range");
294 else
295 {
296 OCTAVE_LOCAL_BUFFER (double, w, 2*n)octave_local_buffer<double> _buffer_w (2*n); double *w =
_buffer_w
;
297
298 F77_XFCN (dchshx, DCHSHX, (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"
, "dchshx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dchshx_ (n, chol_mat.fortran_vec (), chol_mat.rows (), i +
1, j + 1, w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
299 i + 1, j + 1, w))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"
, "dchshx_"); else octave_rethrow_exception (); } else { octave_interrupt_immediately
++; dchshx_ (n, chol_mat.fortran_vec (), chol_mat.rows (), i +
1, j + 1, w); octave_interrupt_immediately--; octave_restore_current_context
(saved_context); } } while (0)
;
300 }
301}
302
303#else
304
305void
306CHOL::update (const ColumnVector& u)
307{
308 warn_qrupdate_once ();
309
310 octave_idx_type n = chol_mat.rows ();
311
312 if (u.length () == n)
313 {
314 init (chol_mat.transpose () * chol_mat
315 + Matrix (u) * Matrix (u).transpose (), false);
316 }
317 else
318 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
319}
320
321static bool
322singular (const Matrix& a)
323{
324 for (octave_idx_type i = 0; i < a.rows (); i++)
325 if (a(i,i) == 0.0) return true;
326 return false;
327}
328
329octave_idx_type
330CHOL::downdate (const ColumnVector& u)
331{
332 warn_qrupdate_once ();
333
334 octave_idx_type info = -1;
335
336 octave_idx_type n = chol_mat.rows ();
337
338 if (u.length () == n)
339 {
340 if (singular (chol_mat))
341 info = 2;
342 else
343 {
344 info = init (chol_mat.transpose () * chol_mat
345 - Matrix (u) * Matrix (u).transpose (), 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
356CHOL::insert_sym (const ColumnVector& u, octave_idx_type j)
357{
358 warn_qrupdate_once ();
359
360 octave_idx_type info = -1;
361
362 octave_idx_type n = chol_mat.rows ();
363
364 if (u.length () != n + 1)
365 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
366 else if (j < 0 || j > n)
367 (*current_liboctave_error_handler) ("cholinsert: index out of range");
368 else
369 {
370 if (singular (chol_mat))
371 info = 2;
372 else
373 {
374 Matrix a = chol_mat.transpose () * chol_mat;
375 Matrix a1 (n+1, n+1);
376 for (octave_idx_type k = 0; k < n+1; k++)
377 for (octave_idx_type l = 0; l < n+1; l++)
378 {
379 if (l == j)
380 a1(k, l) = u(k);
381 else if (k == j)
382 a1(k, l) = u(l);
383 else
384 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
385 }
386 info = init (a1, false);
387 if (info) info = 1;
388 }
389 }
390
391 return info;
392}
393
394void
395CHOL::delete_sym (octave_idx_type j)
396{
397 warn_qrupdate_once ();
398
399 octave_idx_type n = chol_mat.rows ();
400
401 if (j < 0 || j > n-1)
402 (*current_liboctave_error_handler) ("choldelete: index out of range");
403 else
404 {
405 Matrix a = chol_mat.transpose () * chol_mat;
406 a.delete_elements (1, idx_vector (j));
407 a.delete_elements (0, idx_vector (j));
408 init (a, false);
409 }
410}
411
412void
413CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
414{
415 warn_qrupdate_once ();
416
417 octave_idx_type n = chol_mat.rows ();
418
419 if (i < 0 || i > n-1 || j < 0 || j > n-1)
420 (*current_liboctave_error_handler) ("cholshift: index out of range");
421 else
422 {
423 Matrix a = chol_mat.transpose () * chol_mat;
424 Array<octave_idx_type> p (dim_vector (n, 1));
425 for (octave_idx_type k = 0; k < n; k++) p(k) = k;
426 if (i < j)
427 {
428 for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
429 p(j) = i;
430 }
431 else if (j < i)
432 {
433 p(j) = i;
434 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
435 }
436
437 init (a.index (idx_vector (p), idx_vector (p)), false);
438 }
439}
440
441#endif
442
443Matrix
444chol2inv (const Matrix& r)
445{
446 return chol2inv_internal (r);
447}