Mercurial > hg > octave-nkf
annotate liboctave/numeric/CmplxCHOL.cc @ 20830:b65888ec820e draft default tip gccjit
dmalcom gcc jit import
author | Stefan Mahr <dac922@gmx.de> |
---|---|
date | Fri, 27 Feb 2015 16:59:36 +0100 |
parents | dcfbf4c1c3c8 |
children |
rev | line source |
---|---|
457 | 1 /* |
2 | |
19898
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
17769
diff
changeset
|
3 Copyright (C) 1994-2015 John W. Eaton |
11523 | 4 Copyright (C) 2008-2009 Jaroslav Hajek |
457 | 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 | |
7016 | 10 Free Software Foundation; either version 3 of the License, or (at your |
11 option) any later version. | |
457 | 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 | |
7016 | 19 along with Octave; see the file COPYING. If not, see |
20 <http://www.gnu.org/licenses/>. | |
457 | 21 |
22 */ | |
23 | |
24 #ifdef HAVE_CONFIG_H | |
1192 | 25 #include <config.h> |
457 | 26 #endif |
27 | |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
28 #include <vector> |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
29 |
6486 | 30 #include "dMatrix.h" |
31 #include "dRowVector.h" | |
457 | 32 #include "CmplxCHOL.h" |
1847 | 33 #include "f77-fcn.h" |
457 | 34 #include "lo-error.h" |
8377
25bc2d31e1bf
improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents:
7725
diff
changeset
|
35 #include "oct-locbuf.h" |
9862
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
8562
diff
changeset
|
36 #include "oct-norm.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 |
11518 | 44 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, |
45 const octave_idx_type&, Complex*, | |
46 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
47 F77_CHAR_ARG_LEN_DECL); |
5340 | 48 F77_RET_T |
11518 | 49 F77_FUNC (zpotri, ZPOTRI) (F77_CONST_CHAR_ARG_DECL, |
50 const octave_idx_type&, Complex*, | |
51 const octave_idx_type&, octave_idx_type& | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
52 F77_CHAR_ARG_LEN_DECL); |
6486 | 53 |
54 F77_RET_T | |
11518 | 55 F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL, |
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); | |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
60 #ifdef HAVE_QRUPDATE |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
61 |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
62 F77_RET_T |
11495 | 63 F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*, |
64 const octave_idx_type&, Complex*, double*); | |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
65 |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
66 F77_RET_T |
11495 | 67 F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*, |
68 const octave_idx_type&, Complex*, double*, | |
69 octave_idx_type&); | |
70 | |
71 F77_RET_T | |
72 F77_FUNC (zchinx, ZCHINX) (const octave_idx_type&, Complex*, | |
73 const octave_idx_type&, const octave_idx_type&, | |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
74 Complex*, double*, octave_idx_type&); |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
75 |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
76 F77_RET_T |
11495 | 77 F77_FUNC (zchdex, ZCHDEX) (const octave_idx_type&, Complex*, |
78 const octave_idx_type&, const octave_idx_type&, | |
79 double*); | |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
80 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
81 F77_RET_T |
11495 | 82 F77_FUNC (zchshx, ZCHSHX) (const octave_idx_type&, Complex*, |
83 const octave_idx_type&, const octave_idx_type&, | |
84 const octave_idx_type&, Complex*, double*); | |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
85 #endif |
457 | 86 } |
87 | |
5275 | 88 octave_idx_type |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
89 ComplexCHOL::init (const ComplexMatrix& a, bool upper, bool calc_cond) |
457 | 90 { |
5275 | 91 octave_idx_type a_nr = a.rows (); |
92 octave_idx_type a_nc = a.cols (); | |
457 | 93 |
1944 | 94 if (a_nr != a_nc) |
95 { | |
96 (*current_liboctave_error_handler) | |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
97 ("ComplexCHOL requires square matrix"); |
1944 | 98 return -1; |
99 } | |
457 | 100 |
5275 | 101 octave_idx_type n = a_nc; |
102 octave_idx_type info; | |
1944 | 103 |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
104 is_upper = upper; |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
105 |
9862
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
8562
diff
changeset
|
106 chol_mat.clear (n, n); |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
107 if (is_upper) |
9862
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
8562
diff
changeset
|
108 { |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
109 for (octave_idx_type j = 0; j < n; j++) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
110 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
111 for (octave_idx_type i = 0; i <= j; i++) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
112 chol_mat.xelem (i, j) = a (i, j); |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
113 for (octave_idx_type i = j + 1; i < n; i++) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
114 chol_mat.xelem (i, j) = 0.0; |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
115 } |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
116 } |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
117 else |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
118 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
119 for (octave_idx_type j = 0; j < n; j++) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
120 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
121 for (octave_idx_type i = 0; i < j; i++) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
122 chol_mat.xelem (i, j) = 0.0; |
20741
dcfbf4c1c3c8
eliminate trailing whitespace and tabs from sources
John W. Eaton <jwe@octave.org>
parents:
20672
diff
changeset
|
123 for (octave_idx_type i = j; i < n; i++) |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
124 chol_mat.xelem (i, j) = a (i, j); |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
125 } |
9862
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
8562
diff
changeset
|
126 } |
1944 | 127 Complex *h = chol_mat.fortran_vec (); |
457 | 128 |
6486 | 129 // Calculate the norm of the matrix, for later use. |
130 double anorm = 0; | |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
131 if (calc_cond) |
9862
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
8562
diff
changeset
|
132 anorm = xnorm (a, 1); |
6486 | 133 |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
134 if (is_upper) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
135 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
136 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
137 n, h, n, info |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
138 F77_CHAR_ARG_LEN (1))); |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
139 } |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
140 else |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
141 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
142 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
143 n, h, n, info |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
144 F77_CHAR_ARG_LEN (1))); |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
145 } |
457 | 146 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
147 xrcond = 0.0; |
9862
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
8562
diff
changeset
|
148 if (info > 0) |
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
8562
diff
changeset
|
149 chol_mat.resize (info - 1, info - 1); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
150 else if (calc_cond) |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
151 { |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
152 octave_idx_type zpocon_info = 0; |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
153 |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
154 // Now calculate the condition number for non-singular matrix. |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
155 Array<Complex> z (dim_vector (2*n, 1)); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
156 Complex *pz = z.fortran_vec (); |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
157 Array<double> rz (dim_vector (n, 1)); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
158 double *prz = rz.fortran_vec (); |
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
159 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
160 n, anorm, xrcond, pz, prz, zpocon_info |
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
161 F77_CHAR_ARG_LEN (1))); |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
162 |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
163 if (zpocon_info != 0) |
10314
07ebe522dac2
untabify liboctave C++ sources
John W. Eaton <jwe@octave.org>
parents:
10158
diff
changeset
|
164 info = -1; |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
165 } |
1944 | 166 |
167 return info; | |
457 | 168 } |
169 | |
5340 | 170 static ComplexMatrix |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
171 chol2inv_internal (const ComplexMatrix& r, bool is_upper = true) |
5340 | 172 { |
173 ComplexMatrix retval; | |
174 | |
175 octave_idx_type r_nr = r.rows (); | |
176 octave_idx_type r_nc = r.cols (); | |
177 | |
178 if (r_nr == r_nc) | |
179 { | |
180 octave_idx_type n = r_nc; | |
181 octave_idx_type info; | |
182 | |
183 ComplexMatrix tmp = r; | |
184 | |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
185 if (is_upper) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
186 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
187 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n, |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
188 tmp.fortran_vec (), n, info |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
189 F77_CHAR_ARG_LEN (1))); |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
190 } |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
191 else |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
192 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
193 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n, |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
194 tmp.fortran_vec (), n, info |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
195 F77_CHAR_ARG_LEN (1))); |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
196 } |
5340 | 197 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
198 // 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
|
199 // faster for that matter :-)), please let me know! |
5340 | 200 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
201 if (n > 1) |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
202 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
203 if (is_upper) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
204 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
205 for (octave_idx_type j = 0; j < r_nc; j++) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
206 for (octave_idx_type i = j+1; i < r_nr; i++) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
207 tmp.xelem (i, j) = tmp.xelem (j, i); |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
208 } |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
209 else |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
210 { |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
211 for (octave_idx_type j = 0; j < r_nc; j++) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
212 for (octave_idx_type i = j+1; i < r_nr; i++) |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
213 tmp.xelem (j, i) = tmp.xelem (i, j); |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
214 } |
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
215 } |
5340 | 216 |
7482
29980c6b8604
don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents:
7017
diff
changeset
|
217 retval = tmp; |
5340 | 218 } |
219 else | |
220 (*current_liboctave_error_handler) ("chol2inv requires square matrix"); | |
221 | |
222 return retval; | |
223 } | |
224 | |
225 // Compute the inverse of a matrix using the Cholesky factorization. | |
226 ComplexMatrix | |
227 ComplexCHOL::inverse (void) const | |
228 { | |
20672
5ce959c55cc0
Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents:
20442
diff
changeset
|
229 return chol2inv_internal (chol_mat, is_upper); |
5340 | 230 } |
231 | |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
232 void |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
233 ComplexCHOL::set (const ComplexMatrix& R) |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
234 { |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
235 if (R.is_square ()) |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
236 chol_mat = R; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
237 else |
7559
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
238 (*current_liboctave_error_handler) ("CHOL requires square matrix"); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
239 } |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
240 |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
241 #ifdef HAVE_QRUPDATE |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
242 |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
243 void |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
244 ComplexCHOL::update (const ComplexColumnVector& u) |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
245 { |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
246 octave_idx_type n = chol_mat.rows (); |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
247 |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19898
diff
changeset
|
248 if (u.numel () == n) |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
249 { |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
250 ComplexColumnVector utmp = u; |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
251 |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
252 OCTAVE_LOCAL_BUFFER (double, rw, n); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
253 |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
254 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
|
255 utmp.fortran_vec (), rw)); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
256 } |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
257 else |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
258 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
259 } |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
260 |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
261 octave_idx_type |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
262 ComplexCHOL::downdate (const ComplexColumnVector& u) |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
263 { |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
264 octave_idx_type info = -1; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
265 |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
266 octave_idx_type n = chol_mat.rows (); |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
267 |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19898
diff
changeset
|
268 if (u.numel () == n) |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
269 { |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
270 ComplexColumnVector utmp = u; |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
271 |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
272 OCTAVE_LOCAL_BUFFER (double, rw, n); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
273 |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
274 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
|
275 utmp.fortran_vec (), rw, info)); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
276 } |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
277 else |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
278 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch"); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
279 |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
280 return info; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
281 } |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7482
diff
changeset
|
282 |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
283 octave_idx_type |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
284 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
|
285 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
286 octave_idx_type info = -1; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
287 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
288 octave_idx_type n = chol_mat.rows (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
289 |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19898
diff
changeset
|
290 if (u.numel () != n + 1) |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
291 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch"); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
292 else if (j < 0 || j > n) |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
293 (*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
|
294 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
295 { |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
296 ComplexColumnVector utmp = u; |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
297 |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
298 OCTAVE_LOCAL_BUFFER (double, rw, n); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
299 |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
300 chol_mat.resize (n+1, n+1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
301 |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
302 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
|
303 j + 1, utmp.fortran_vec (), rw, info)); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
304 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
305 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
306 return info; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
307 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
308 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
309 void |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
310 ComplexCHOL::delete_sym (octave_idx_type j) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
311 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
312 octave_idx_type n = chol_mat.rows (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
313 |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
314 if (j < 0 || j > n-1) |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
315 (*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
|
316 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
317 { |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
318 OCTAVE_LOCAL_BUFFER (double, rw, n); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
319 |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
320 F77_XFCN (zchdex, ZCHDEX, (n, chol_mat.fortran_vec (), chol_mat.rows (), |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
321 j + 1, rw)); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
322 |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
323 chol_mat.resize (n-1, n-1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
324 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
325 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
326 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
327 void |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
328 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
|
329 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
330 octave_idx_type n = chol_mat.rows (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
331 |
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
332 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
|
333 (*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
|
334 else |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
335 { |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
336 OCTAVE_LOCAL_BUFFER (Complex, w, n); |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
337 OCTAVE_LOCAL_BUFFER (double, rw, n); |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
338 |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
339 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
|
340 i + 1, j + 1, w, rw)); |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
341 } |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
342 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7559
diff
changeset
|
343 |
8562
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
344 #else |
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 void |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
347 ComplexCHOL::update (const ComplexColumnVector& u) |
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 warn_qrupdate_once (); |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
350 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
351 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
|
352 |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19898
diff
changeset
|
353 if (u.numel () == n) |
8562
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
354 { |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
355 init (chol_mat.hermitian () * chol_mat |
8562
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
356 + ComplexMatrix (u) * ComplexMatrix (u).hermitian (), false); |
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 else |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
359 (*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
|
360 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
361 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
362 static bool |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
363 singular (const ComplexMatrix& a) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
364 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
365 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
|
366 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
|
367 return false; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
368 } |
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 octave_idx_type |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
371 ComplexCHOL::downdate (const ComplexColumnVector& u) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
372 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
373 warn_qrupdate_once (); |
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 octave_idx_type info = -1; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
376 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
377 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
|
378 |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19898
diff
changeset
|
379 if (u.numel () == n) |
8562
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
380 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
381 if (singular (chol_mat)) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
382 info = 2; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
383 else |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
384 { |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
385 info = init (chol_mat.hermitian () * chol_mat |
17769
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
386 - ComplexMatrix (u) * ComplexMatrix (u).hermitian (), |
49a5a4be04a1
maint: Use GNU style coding conventions for code in liboctave/
Rik <rik@octave.org>
parents:
17744
diff
changeset
|
387 false); |
8562
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
388 if (info) info = 1; |
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 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
391 else |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
392 (*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
|
393 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
394 return info; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
395 } |
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 octave_idx_type |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
398 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
|
399 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
400 warn_qrupdate_once (); |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
401 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
402 octave_idx_type info = -1; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
403 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
404 octave_idx_type n = chol_mat.rows (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
405 |
20442
a9574e3c6e9e
Deprecate Array::length() and Sparse::length() in favour of ::numel().
Carnë Draug <carandraug@octave.org>
parents:
19898
diff
changeset
|
406 if (u.numel () != n + 1) |
8562
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
407 (*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
|
408 else if (j < 0 || j > n) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
409 (*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
|
410 else |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
411 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
412 if (singular (chol_mat)) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
413 info = 2; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
414 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
|
415 info = 3; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
416 else |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
417 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
418 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
|
419 ComplexMatrix a1 (n+1, n+1); |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
420 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
|
421 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
|
422 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
423 if (l == j) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
424 a1(k, l) = u(k); |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
425 else if (k == j) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
426 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
|
427 else |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
428 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
|
429 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
430 info = init (a1, false); |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
431 if (info) info = 1; |
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 } |
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 return info; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
436 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
437 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
438 void |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
439 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
|
440 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
441 warn_qrupdate_once (); |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
442 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
443 octave_idx_type n = chol_mat.rows (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
444 |
8562
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
445 if (j < 0 || j > n-1) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
446 (*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
|
447 else |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
448 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
449 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
|
450 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
|
451 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
|
452 init (a, false); |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
453 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
454 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
455 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
456 void |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
457 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
|
458 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
459 warn_qrupdate_once (); |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
460 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
461 octave_idx_type n = chol_mat.rows (); |
11586
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
462 |
12df7854fa7c
strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents:
11570
diff
changeset
|
463 if (i < 0 || i > n-1 || j < 0 || j > n-1) |
8562
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
464 (*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
|
465 else |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
466 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
467 ComplexMatrix a = chol_mat.hermitian () * chol_mat; |
11570
57632dea2446
attempt better backward compatibility for Array constructors
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
468 Array<octave_idx_type> p (dim_vector (n, 1)); |
8562
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
469 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
|
470 if (i < j) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
471 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
472 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
|
473 p(j) = i; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
474 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
475 else if (j < i) |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
476 { |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
477 p(j) = i; |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
478 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
|
479 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
480 |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
481 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
|
482 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
483 } |
a6edd5c23cb5
use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents:
8547
diff
changeset
|
484 |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
485 #endif |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8377
diff
changeset
|
486 |
5340 | 487 ComplexMatrix |
488 chol2inv (const ComplexMatrix& r) | |
489 { | |
490 return chol2inv_internal (r); | |
491 } |