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