File: | liboctave/numeric/dbleCHOL.cc |
Location: | line 124, column 12 |
Description: | The left operand of '>' is a garbage value |
1 | /* | |||
2 | ||||
3 | Copyright (C) 1994-2013 John W. Eaton | |||
4 | Copyright (C) 2008-2009 Jaroslav Hajek | |||
5 | ||||
6 | This file is part of Octave. | |||
7 | ||||
8 | Octave is free software; you can redistribute it and/or modify it | |||
9 | under the terms of the GNU General Public License as published by the | |||
10 | Free Software Foundation; either version 3 of the License, or (at your | |||
11 | option) any later version. | |||
12 | ||||
13 | Octave is distributed in the hope that it will be useful, but WITHOUT | |||
14 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |||
15 | FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |||
16 | for more details. | |||
17 | ||||
18 | You should have received a copy of the GNU General Public License | |||
19 | along 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 | ||||
40 | extern "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 | ||||
89 | octave_idx_type | |||
90 | CHOL::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) | |||
| ||||
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; | |||
103 | ||||
104 | chol_mat.clear (n, n); | |||
105 | for (octave_idx_type j = 0; j < n; j++) | |||
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) | |||
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) | |||
| ||||
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 | ||||
146 | static Matrix | |||
147 | chol2inv_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. | |||
186 | Matrix | |||
187 | CHOL::inverse (void) const | |||
188 | { | |||
189 | return chol2inv_internal (chol_mat); | |||
190 | } | |||
191 | ||||
192 | void | |||
193 | CHOL::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 | ||||
203 | void | |||
204 | CHOL::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 | ||||
221 | octave_idx_type | |||
222 | CHOL::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 | ||||
243 | octave_idx_type | |||
244 | CHOL::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 | ||||
269 | void | |||
270 | CHOL::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 | ||||
287 | void | |||
288 | CHOL::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 | ||||
305 | void | |||
306 | CHOL::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 | ||||
321 | static bool | |||
322 | singular (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 | ||||
329 | octave_idx_type | |||
330 | CHOL::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 | ||||
355 | octave_idx_type | |||
356 | CHOL::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 | ||||
394 | void | |||
395 | CHOL::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 | ||||
412 | void | |||
413 | CHOL::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 | ||||
443 | Matrix | |||
444 | chol2inv (const Matrix& r) | |||
445 | { | |||
446 | return chol2inv_internal (r); | |||
447 | } |