File: | liboctave/numeric/CmplxCHOL.cc |
Location: | line 123, 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 "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 | ||||
41 | extern "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 | ||||
88 | octave_idx_type | |||
89 | ComplexCHOL::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) | |||
| ||||
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; | |||
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 | 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) | |||
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) | |||
| ||||
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 | ||||
145 | static ComplexMatrix | |||
146 | chol2inv_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. | |||
181 | ComplexMatrix | |||
182 | ComplexCHOL::inverse (void) const | |||
183 | { | |||
184 | return chol2inv_internal (chol_mat); | |||
185 | } | |||
186 | ||||
187 | void | |||
188 | ComplexCHOL::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 | ||||
198 | void | |||
199 | ComplexCHOL::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 | ||||
216 | octave_idx_type | |||
217 | ComplexCHOL::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 | ||||
238 | octave_idx_type | |||
239 | ComplexCHOL::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 | ||||
264 | void | |||
265 | ComplexCHOL::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 | ||||
282 | void | |||
283 | ComplexCHOL::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 | ||||
301 | void | |||
302 | ComplexCHOL::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 | ||||
317 | static bool | |||
318 | singular (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 | ||||
325 | octave_idx_type | |||
326 | ComplexCHOL::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 | ||||
352 | octave_idx_type | |||
353 | ComplexCHOL::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 | ||||
393 | void | |||
394 | ComplexCHOL::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 | ||||
411 | void | |||
412 | ComplexCHOL::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 | ||||
442 | ComplexMatrix | |||
443 | chol2inv (const ComplexMatrix& r) | |||
444 | { | |||
445 | return chol2inv_internal (r); | |||
446 | } |