annotate liboctave/numeric/chol.cc @ 27919:1891570abac8

update Octave Project Developers copyright for the new year In files that have the "Octave Project Developers" copyright notice, update for 2020.
author John W. Eaton <jwe@octave.org>
date Mon, 06 Jan 2020 22:29:51 -0500
parents b442ec6dda5c
children bd51beb6205e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
1 /*
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
2
27919
1891570abac8 update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 27918
diff changeset
3 Copyright (C) 1994-2020 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27102
diff changeset
4
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27102
diff changeset
5 See the file COPYRIGHT.md in the top-level directory of this distribution
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27102
diff changeset
6 or <https://octave.org/COPYRIGHT.html/>.
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 27102
diff changeset
7
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
8
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
9 This file is part of Octave.
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
10
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23708
diff changeset
11 Octave is free software: you can redistribute it and/or modify it
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
12 under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23708
diff changeset
13 the Free Software Foundation, either version 3 of the License, or
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
14 (at your option) any later version.
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
15
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
16 Octave is distributed in the hope that it will be useful, but
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
17 WITHOUT ANY WARRANTY; without even the implied warranty of
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22402
diff changeset
19 GNU General Public License for more details.
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
20
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
21 You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6486
diff changeset
22 along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23708
diff changeset
23 <https://www.gnu.org/licenses/>.
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
24
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
25 */
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
26
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21301
diff changeset
27 #if defined (HAVE_CONFIG_H)
21301
40de9f8f23a6 Use '#include "config.h"' rather than <config.h>.
Rik <rik@octave.org>
parents: 21285
diff changeset
28 # include "config.h"
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
29 #endif
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
30
23475
d691ed308237 maint: Clean up #includes in liboctave/numeric directory.
Rik <rik@octave.org>
parents: 23220
diff changeset
31 #include "Array.h"
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
32 #include "CColVector.h"
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
33 #include "CMatrix.h"
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
34 #include "chol.h"
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
35 #include "dColVector.h"
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
36 #include "dMatrix.h"
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
37 #include "fCColVector.h"
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
38 #include "fCMatrix.h"
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
39 #include "fColVector.h"
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
40 #include "fMatrix.h"
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
41 #include "lo-error.h"
22322
93b3cdd36854 move most f77 function decls to separate header files
John W. Eaton <jwe@octave.org>
parents: 22317
diff changeset
42 #include "lo-lapack-proto.h"
93b3cdd36854 move most f77 function decls to separate header files
John W. Eaton <jwe@octave.org>
parents: 22317
diff changeset
43 #include "lo-qrupdate-proto.h"
8377
25bc2d31e1bf improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents: 7725
diff changeset
44 #include "oct-locbuf.h"
9862
c0aeedd8fb86 improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents: 8562
diff changeset
45 #include "oct-norm.h"
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
46
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
47 #if ! defined (HAVE_QRUPDATE)
21285
d3b265a83adc make building without qrupdate work again
John W. Eaton <jwe@octave.org>
parents: 21270
diff changeset
48 # include "qr.h"
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
49 #endif
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
50
22317
6ca3acf5fad8 move some new numeric classes to namespace octave::math
John W. Eaton <jwe@octave.org>
parents: 22274
diff changeset
51 namespace octave
6ca3acf5fad8 move some new numeric classes to namespace octave::math
John W. Eaton <jwe@octave.org>
parents: 22274
diff changeset
52 {
27101
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
53 static Matrix
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
54 chol2inv_internal (const Matrix& r, bool is_upper = true)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
55 {
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
56 Matrix retval;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
57
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
58 octave_idx_type r_nr = r.rows ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
59 octave_idx_type r_nc = r.cols ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
60
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
61 if (r_nr != r_nc)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
62 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
63
27102
84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 27101
diff changeset
64 F77_INT n = to_f77_int (r_nc);
27101
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
65 F77_INT info;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
66
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
67 Matrix tmp = r;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
68 double *v = tmp.fortran_vec ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
69
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
70 if (is_upper)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
71 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
72 v, n, info
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
73 F77_CHAR_ARG_LEN (1)));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
74 else
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
75 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
76 v, n, info
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
77 F77_CHAR_ARG_LEN (1)));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
78
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
79 // FIXME: Should we check info exit value and possibly report an error?
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
80
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
81 // If someone thinks of a more graceful way of doing this
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
82 // (or faster for that matter :-)), please let me know!
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
83
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
84 if (n > 1)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
85 {
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
86 if (is_upper)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
87 for (octave_idx_type j = 0; j < r_nc; j++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
88 for (octave_idx_type i = j+1; i < r_nr; i++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
89 tmp.xelem (i, j) = tmp.xelem (j, i);
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
90 else
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
91 for (octave_idx_type j = 0; j < r_nc; j++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
92 for (octave_idx_type i = j+1; i < r_nr; i++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
93 tmp.xelem (j, i) = tmp.xelem (i, j);
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
94 }
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
95
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
96 retval = tmp;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
97
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
98 return retval;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
99 }
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
100
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
101 static FloatMatrix
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
102 chol2inv_internal (const FloatMatrix& r, bool is_upper = true)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
103 {
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
104 FloatMatrix retval;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
105
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
106 octave_idx_type r_nr = r.rows ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
107 octave_idx_type r_nc = r.cols ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
108
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
109 if (r_nr != r_nc)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
110 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
111
27102
84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 27101
diff changeset
112 F77_INT n = to_f77_int (r_nc);
27101
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
113 F77_INT info;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
114
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
115 FloatMatrix tmp = r;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
116 float *v = tmp.fortran_vec ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
117
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
118 if (is_upper)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
119 F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
120 v, n, info
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
121 F77_CHAR_ARG_LEN (1)));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
122 else
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
123 F77_XFCN (spotri, SPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
124 v, n, info
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
125 F77_CHAR_ARG_LEN (1)));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
126
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
127 // FIXME: Should we check info exit value and possibly report an error?
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
128
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
129 // If someone thinks of a more graceful way of doing this (or
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
130 // faster for that matter :-)), please let me know!
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
131
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
132 if (n > 1)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
133 {
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
134 if (is_upper)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
135 for (octave_idx_type j = 0; j < r_nc; j++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
136 for (octave_idx_type i = j+1; i < r_nr; i++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
137 tmp.xelem (i, j) = tmp.xelem (j, i);
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
138 else
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
139 for (octave_idx_type j = 0; j < r_nc; j++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
140 for (octave_idx_type i = j+1; i < r_nr; i++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
141 tmp.xelem (j, i) = tmp.xelem (i, j);
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
142 }
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
143
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
144 retval = tmp;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
145
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
146 return retval;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
147 }
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
148
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
149 static ComplexMatrix
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
150 chol2inv_internal (const ComplexMatrix& r, bool is_upper = true)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
151 {
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
152 ComplexMatrix retval;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
153
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
154 octave_idx_type r_nr = r.rows ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
155 octave_idx_type r_nc = r.cols ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
156
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
157 if (r_nr != r_nc)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
158 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
159
27102
84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 27101
diff changeset
160 F77_INT n = to_f77_int (r_nc);
27101
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
161 F77_INT info;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
162
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
163 ComplexMatrix tmp = r;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
164
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
165 if (is_upper)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
166 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
167 F77_DBLE_CMPLX_ARG (tmp.fortran_vec ()), n, info
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
168 F77_CHAR_ARG_LEN (1)));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
169 else
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
170 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
171 F77_DBLE_CMPLX_ARG (tmp.fortran_vec ()), n, info
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
172 F77_CHAR_ARG_LEN (1)));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
173
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
174 // If someone thinks of a more graceful way of doing this (or
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
175 // faster for that matter :-)), please let me know!
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
176
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
177 if (n > 1)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
178 {
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
179 if (is_upper)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
180 for (octave_idx_type j = 0; j < r_nc; j++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
181 for (octave_idx_type i = j+1; i < r_nr; i++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
182 tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
183 else
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
184 for (octave_idx_type j = 0; j < r_nc; j++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
185 for (octave_idx_type i = j+1; i < r_nr; i++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
186 tmp.xelem (j, i) = std::conj (tmp.xelem (i, j));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
187 }
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
188
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
189 retval = tmp;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
190
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
191 return retval;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
192 }
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
193
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
194 static FloatComplexMatrix
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
195 chol2inv_internal (const FloatComplexMatrix& r, bool is_upper = true)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
196 {
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
197 FloatComplexMatrix retval;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
198
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
199 octave_idx_type r_nr = r.rows ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
200 octave_idx_type r_nc = r.cols ();
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
201
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
202 if (r_nr != r_nc)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
203 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
204
27102
84ff9953faa1 where possible, eliminate octave:: namespace qualifier inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 27101
diff changeset
205 F77_INT n = to_f77_int (r_nc);
27101
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
206 F77_INT info;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
207
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
208 FloatComplexMatrix tmp = r;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
209
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
210 if (is_upper)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
211 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
212 F77_CMPLX_ARG (tmp.fortran_vec ()), n, info
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
213 F77_CHAR_ARG_LEN (1)));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
214 else
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
215 F77_XFCN (cpotri, CPOTRI, (F77_CONST_CHAR_ARG2 ("L", 1), n,
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
216 F77_CMPLX_ARG (tmp.fortran_vec ()), n, info
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
217 F77_CHAR_ARG_LEN (1)));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
218
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
219 // If someone thinks of a more graceful way of doing this (or
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
220 // faster for that matter :-)), please let me know!
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
221
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
222 if (n > 1)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
223 {
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
224 if (is_upper)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
225 for (octave_idx_type j = 0; j < r_nc; j++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
226 for (octave_idx_type i = j+1; i < r_nr; i++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
227 tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
228 else
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
229 for (octave_idx_type j = 0; j < r_nc; j++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
230 for (octave_idx_type i = j+1; i < r_nr; i++)
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
231 tmp.xelem (j, i) = std::conj (tmp.xelem (i, j));
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
232 }
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
233
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
234 retval = tmp;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
235
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
236 return retval;
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
237 }
097774bed4ed move some static functions inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 26385
diff changeset
238
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
239 namespace math
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
240 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
241 template <typename T>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
242 T
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
243 chol2inv (const T& r)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
244 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
245 return chol2inv_internal (r);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
246 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
247
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
248 // Compute the inverse of a matrix using the Cholesky factorization.
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
249 template <typename T>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
250 T
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
251 chol<T>::inverse (void) const
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
252 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
253 return chol2inv_internal (chol_mat, is_upper);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
254 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
255
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
256 template <typename T>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
257 void
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
258 chol<T>::set (const T& R)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
259 {
23593
a8361bc2361a maint: Deprecate is_square and replace with issquare.
Rik <rik@octave.org>
parents: 23475
diff changeset
260 if (! R.issquare ())
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
261 (*current_liboctave_error_handler) ("chol: requires square matrix");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
262
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
263 chol_mat = R;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
264 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
265
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
266 #if ! defined (HAVE_QRUPDATE)
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
267
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
268 template <typename T>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
269 void
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
270 chol<T>::update (const VT& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
271 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
272 warn_qrupdate_once ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
273
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
274 octave_idx_type n = chol_mat.rows ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
275
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
276 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
277 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
278
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
279 init (chol_mat.hermitian () * chol_mat + T (u) * T (u).hermitian (),
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
280 true, false);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
281 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
282
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
283 template <typename T>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
284 bool
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
285 singular (const T& a)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
286 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
287 static typename T::element_type zero (0);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
288 for (octave_idx_type i = 0; i < a.rows (); i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
289 if (a(i,i) == zero) return true;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
290 return false;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
291 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
292
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
293 template <typename T>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
294 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
295 chol<T>::downdate (const VT& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
296 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
297 warn_qrupdate_once ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
298
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
299 octave_idx_type info = -1;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
300
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
301 octave_idx_type n = chol_mat.rows ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
302
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
303 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
304 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
305
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
306 if (singular (chol_mat))
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
307 info = 2;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
308 else
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
309 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
310 info = init (chol_mat.hermitian () * chol_mat
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
311 - T (u) * T (u).hermitian (), true, false);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
312 if (info) info = 1;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
313 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
314
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
315 return info;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
316 }
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
317
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
318 template <typename T>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
319 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
320 chol<T>::insert_sym (const VT& u, octave_idx_type j)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
321 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
322 static typename T::element_type zero (0);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
323
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
324 warn_qrupdate_once ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
325
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
326 octave_idx_type info = -1;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
327
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
328 octave_idx_type n = chol_mat.rows ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
329
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
330 if (u.numel () != n + 1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
331 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
332 if (j < 0 || j > n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
333 (*current_liboctave_error_handler) ("cholinsert: index out of range");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
334
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
335 if (singular (chol_mat))
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
336 info = 2;
23708
750e42a35adc Use imag, real, arg from std library for Complex types.
Rik <rik@octave.org>
parents: 23696
diff changeset
337 else if (std::imag (u(j)) != zero)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
338 info = 3;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
339 else
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
340 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
341 T a = chol_mat.hermitian () * chol_mat;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
342 T a1 (n+1, n+1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
343 for (octave_idx_type k = 0; k < n+1; k++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
344 for (octave_idx_type l = 0; l < n+1; l++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
345 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
346 if (l == j)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
347 a1(k, l) = u(k);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
348 else if (k == j)
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
349 a1(k, l) = math::conj (u(l));
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
350 else
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
351 a1(k, l) = a(k < j ? k : k-1, l < j ? l : l-1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
352 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
353 info = init (a1, true, false);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
354 if (info) info = 1;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
355 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
356
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
357 return info;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
358 }
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
359
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
360 template <typename T>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
361 void
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
362 chol<T>::delete_sym (octave_idx_type j)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
363 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
364 warn_qrupdate_once ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
365
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
366 octave_idx_type n = chol_mat.rows ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
367
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
368 if (j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
369 (*current_liboctave_error_handler) ("choldelete: index out of range");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
370
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
371 T a = chol_mat.hermitian () * chol_mat;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
372 a.delete_elements (1, idx_vector (j));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
373 a.delete_elements (0, idx_vector (j));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
374 init (a, true, false);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
375 }
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
376
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
377 template <typename T>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
378 void
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
379 chol<T>::shift_sym (octave_idx_type i, octave_idx_type j)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
380 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
381 warn_qrupdate_once ();
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
382
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
383 octave_idx_type n = chol_mat.rows ();
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
384
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
385 if (i < 0 || i > n-1 || j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
386 (*current_liboctave_error_handler) ("cholshift: index out of range");
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
387
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
388 T a = chol_mat.hermitian () * chol_mat;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
389 Array<octave_idx_type> p (dim_vector (n, 1));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
390 for (octave_idx_type k = 0; k < n; k++) p(k) = k;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
391 if (i < j)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
392 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
393 for (octave_idx_type k = i; k < j; k++) p(k) = k+1;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
394 p(j) = i;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
395 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
396 else if (j < i)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
397 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
398 p(j) = i;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
399 for (octave_idx_type k = j+1; k < i+1; k++) p(k) = k-1;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
400 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
401
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
402 init (a.index (idx_vector (p), idx_vector (p)), true, false);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
403 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
404
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
405 #endif
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
406
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
407 // Specializations.
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
408
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
409 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
410 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
411 chol<Matrix>::init (const Matrix& a, bool upper, bool calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
412 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
413 octave_idx_type a_nr = a.rows ();
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
414 octave_idx_type a_nc = a.cols ();
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
415
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
416 if (a_nr != a_nc)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
417 (*current_liboctave_error_handler) ("chol: requires square matrix");
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
418
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
419 F77_INT n = to_f77_int (a_nc);
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
420 F77_INT info;
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
421
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
422 is_upper = upper;
20462
5ce959c55cc0 Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents: 20232
diff changeset
423
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
424 chol_mat.clear (n, n);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
425 if (is_upper)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
426 for (octave_idx_type j = 0; j < n; j++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
427 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
428 for (octave_idx_type i = 0; i <= j; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
429 chol_mat.xelem (i, j) = a(i, j);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
430 for (octave_idx_type i = j+1; i < n; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
431 chol_mat.xelem (i, j) = 0.0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
432 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
433 else
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
434 for (octave_idx_type j = 0; j < n; j++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
435 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
436 for (octave_idx_type i = 0; i < j; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
437 chol_mat.xelem (i, j) = 0.0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
438 for (octave_idx_type i = j; i < n; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
439 chol_mat.xelem (i, j) = a(i, j);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
440 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
441 double *h = chol_mat.fortran_vec ();
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
442
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
443 // Calculate the norm of the matrix, for later use.
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
444 double anorm = 0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
445 if (calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
446 anorm = xnorm (a, 1);
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
447
20462
5ce959c55cc0 Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents: 20232
diff changeset
448 if (is_upper)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
449 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
20628
48fedd8fbff7 maint: Apply Octave coding style to Cholesky classes
Mike Miller <mtmiller@octave.org>
parents: 20531
diff changeset
450 F77_CHAR_ARG_LEN (1)));
20462
5ce959c55cc0 Propagate 'lower' in chol(a, 'lower') to underlying library function.
PrasannaKumar Muralidharan <prasannatsmkumar@gmail.com>
parents: 20232
diff changeset
451 else
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
452 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
20628
48fedd8fbff7 maint: Apply Octave coding style to Cholesky classes
Mike Miller <mtmiller@octave.org>
parents: 20531
diff changeset
453 F77_CHAR_ARG_LEN (1)));
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
454
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
455 xrcond = 0.0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
456 if (info > 0)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
457 chol_mat.resize (info - 1, info - 1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
458 else if (calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
459 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
460 F77_INT dpocon_info = 0;
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
461
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
462 // Now calculate the condition number for non-singular matrix.
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
463 Array<double> z (dim_vector (3*n, 1));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
464 double *pz = z.fortran_vec ();
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
465 OCTAVE_LOCAL_BUFFER (F77_INT, iz, n);
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
466 if (is_upper)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
467 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
468 n, anorm, xrcond, pz, iz, dpocon_info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
469 F77_CHAR_ARG_LEN (1)));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
470 else
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
471 F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
472 n, anorm, xrcond, pz, iz, dpocon_info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
473 F77_CHAR_ARG_LEN (1)));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
474
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
475 if (dpocon_info != 0)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
476 info = -1;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
477 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
478
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
479 return info;
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
480 }
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
481
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
482 #if defined (HAVE_QRUPDATE)
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
483
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
484 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
485 void
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
486 chol<Matrix>::update (const ColumnVector& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
487 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
488 F77_INT n = to_f77_int (chol_mat.rows ());
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
489
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
490 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
491 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
492
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
493 ColumnVector utmp = u;
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
494
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
495 OCTAVE_LOCAL_BUFFER (double, w, n);
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
496
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
497 F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (), n,
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
498 utmp.fortran_vec (), w));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
499 }
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
500
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
501 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
502 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
503 chol<Matrix>::downdate (const ColumnVector& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
504 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
505 F77_INT info = -1;
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
506
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
507 F77_INT n = to_f77_int (chol_mat.rows ());
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
508
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
509 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
510 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
511
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
512 ColumnVector utmp = u;
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
513
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
514 OCTAVE_LOCAL_BUFFER (double, w, n);
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
515
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
516 F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (), n,
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
517 utmp.fortran_vec (), w, info));
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
518
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
519 return info;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
520 }
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
521
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
522 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
523 octave_idx_type
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
524 chol<Matrix>::insert_sym (const ColumnVector& u, octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
525 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
526 F77_INT info = -1;
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
527
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
528 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
529 F77_INT j = to_f77_int (j_arg);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
530
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
531 if (u.numel () != n + 1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
532 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
533 if (j < 0 || j > n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
534 (*current_liboctave_error_handler) ("cholinsert: index out of range");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
535
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
536 ColumnVector utmp = u;
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
537
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
538 OCTAVE_LOCAL_BUFFER (double, w, n);
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
539
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
540 chol_mat.resize (n+1, n+1);
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
541 F77_INT ldcm = to_f77_int (chol_mat.rows ());
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
542
22969
023026207910 * chol.cc: Fix thinko in previous change.
John W. Eaton <jwe@octave.org>
parents: 22942
diff changeset
543 F77_XFCN (dchinx, DCHINX, (n, chol_mat.fortran_vec (), ldcm,
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
544 j + 1, utmp.fortran_vec (), w, info));
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
545
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
546 return info;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
547 }
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
548
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
549 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
550 void
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
551 chol<Matrix>::delete_sym (octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
552 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
553 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
554 F77_INT j = to_f77_int (j_arg);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
555
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
556 if (j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
557 (*current_liboctave_error_handler) ("choldelete: index out of range");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
558
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
559 OCTAVE_LOCAL_BUFFER (double, w, n);
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
560
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
561 F77_XFCN (dchdex, DCHDEX, (n, chol_mat.fortran_vec (), n, j + 1, w));
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
562
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
563 chol_mat.resize (n-1, n-1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
564 }
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
565
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
566 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
567 void
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
568 chol<Matrix>::shift_sym (octave_idx_type i_arg, octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
569 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
570 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
571 F77_INT i = to_f77_int (i_arg);
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
572 F77_INT j = to_f77_int (j_arg);
11586
12df7854fa7c strip trailing whitespace from source files
John W. Eaton <jwe@octave.org>
parents: 11570
diff changeset
573
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
574 if (i < 0 || i > n-1 || j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
575 (*current_liboctave_error_handler) ("cholshift: index out of range");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
576
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
577 OCTAVE_LOCAL_BUFFER (double, w, 2*n);
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
578
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
579 F77_XFCN (dchshx, DCHSHX, (n, chol_mat.fortran_vec (), n,
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
580 i + 1, j + 1, w));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
581 }
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
582
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
583 #endif
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
584
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
585 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
586 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
587 chol<FloatMatrix>::init (const FloatMatrix& a, bool upper, bool calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
588 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
589 octave_idx_type a_nr = a.rows ();
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
590 octave_idx_type a_nc = a.cols ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
591
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
592 if (a_nr != a_nc)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
593 (*current_liboctave_error_handler) ("chol: requires square matrix");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
594
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
595 F77_INT n = to_f77_int (a_nc);
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
596 F77_INT info;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
597
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
598 is_upper = upper;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
599
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
600 chol_mat.clear (n, n);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
601 if (is_upper)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
602 for (octave_idx_type j = 0; j < n; j++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
603 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
604 for (octave_idx_type i = 0; i <= j; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
605 chol_mat.xelem (i, j) = a(i, j);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
606 for (octave_idx_type i = j+1; i < n; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
607 chol_mat.xelem (i, j) = 0.0f;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
608 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
609 else
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
610 for (octave_idx_type j = 0; j < n; j++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
611 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
612 for (octave_idx_type i = 0; i < j; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
613 chol_mat.xelem (i, j) = 0.0f;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
614 for (octave_idx_type i = j; i < n; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
615 chol_mat.xelem (i, j) = a(i, j);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
616 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
617 float *h = chol_mat.fortran_vec ();
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
618
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
619 // Calculate the norm of the matrix, for later use.
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
620 float anorm = 0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
621 if (calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
622 anorm = xnorm (a, 1);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
623
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
624 if (is_upper)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
625 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
626 F77_CHAR_ARG_LEN (1)));
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
627 else
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
628 F77_XFCN (spotrf, SPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n, h, n, info
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
629 F77_CHAR_ARG_LEN (1)));
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
630
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
631 xrcond = 0.0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
632 if (info > 0)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
633 chol_mat.resize (info - 1, info - 1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
634 else if (calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
635 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
636 F77_INT spocon_info = 0;
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
637
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
638 // Now calculate the condition number for non-singular matrix.
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
639 Array<float> z (dim_vector (3*n, 1));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
640 float *pz = z.fortran_vec ();
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
641 OCTAVE_LOCAL_BUFFER (F77_INT, iz, n);
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
642 if (is_upper)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
643 F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
644 n, anorm, xrcond, pz, iz, spocon_info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
645 F77_CHAR_ARG_LEN (1)));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
646 else
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
647 F77_XFCN (spocon, SPOCON, (F77_CONST_CHAR_ARG2 ("L", 1), n, h,
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
648 n, anorm, xrcond, pz, iz, spocon_info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
649 F77_CHAR_ARG_LEN (1)));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
650
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
651 if (spocon_info != 0)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
652 info = -1;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
653 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
654
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
655 return info;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
656 }
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
657
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21301
diff changeset
658 #if defined (HAVE_QRUPDATE)
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
659
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
660 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
661 void
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
662 chol<FloatMatrix>::update (const FloatColumnVector& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
663 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
664 F77_INT n = to_f77_int (chol_mat.rows ());
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
665
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
666 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
667 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
668
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
669 FloatColumnVector utmp = u;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
670
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
671 OCTAVE_LOCAL_BUFFER (float, w, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
672
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
673 F77_XFCN (sch1up, SCH1UP, (n, chol_mat.fortran_vec (), n,
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
674 utmp.fortran_vec (), w));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
675 }
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
676
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
677 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
678 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
679 chol<FloatMatrix>::downdate (const FloatColumnVector& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
680 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
681 F77_INT info = -1;
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
682
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
683 F77_INT n = to_f77_int (chol_mat.rows ());
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
684
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
685 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
686 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
687
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
688 FloatColumnVector utmp = u;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
689
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
690 OCTAVE_LOCAL_BUFFER (float, w, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
691
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
692 F77_XFCN (sch1dn, SCH1DN, (n, chol_mat.fortran_vec (), n,
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
693 utmp.fortran_vec (), w, info));
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
694
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
695 return info;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
696 }
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
697
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
698 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
699 octave_idx_type
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
700 chol<FloatMatrix>::insert_sym (const FloatColumnVector& u,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
701 octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
702 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
703 F77_INT info = -1;
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
704
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
705 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
706 F77_INT j = to_f77_int (j_arg);
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
707
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
708 if (u.numel () != n + 1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
709 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
710 if (j < 0 || j > n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
711 (*current_liboctave_error_handler) ("cholinsert: index out of range");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
712
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
713 FloatColumnVector utmp = u;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
714
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
715 OCTAVE_LOCAL_BUFFER (float, w, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
716
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
717 chol_mat.resize (n+1, n+1);
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
718 F77_INT ldcm = to_f77_int (chol_mat.rows ());
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
719
22969
023026207910 * chol.cc: Fix thinko in previous change.
John W. Eaton <jwe@octave.org>
parents: 22942
diff changeset
720 F77_XFCN (schinx, SCHINX, (n, chol_mat.fortran_vec (), ldcm,
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
721 j + 1, utmp.fortran_vec (), w, info));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
722
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
723 return info;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
724 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
725
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
726 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
727 void
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
728 chol<FloatMatrix>::delete_sym (octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
729 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
730 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
731 F77_INT j = to_f77_int (j_arg);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
732
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
733 if (j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
734 (*current_liboctave_error_handler) ("choldelete: index out of range");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
735
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
736 OCTAVE_LOCAL_BUFFER (float, w, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
737
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
738 F77_XFCN (schdex, SCHDEX, (n, chol_mat.fortran_vec (), n,
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
739 j + 1, w));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
740
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
741 chol_mat.resize (n-1, n-1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
742 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
743
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
744 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
745 void
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
746 chol<FloatMatrix>::shift_sym (octave_idx_type i_arg, octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
747 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
748 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
749 F77_INT i = to_f77_int (i_arg);
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
750 F77_INT j = to_f77_int (j_arg);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
751
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
752 if (i < 0 || i > n-1 || j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
753 (*current_liboctave_error_handler) ("cholshift: index out of range");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
754
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
755 OCTAVE_LOCAL_BUFFER (float, w, 2*n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
756
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
757 F77_XFCN (schshx, SCHSHX, (n, chol_mat.fortran_vec (), n,
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
758 i + 1, j + 1, w));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
759 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
760
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
761 #endif
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
762
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
763 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
764 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
765 chol<ComplexMatrix>::init (const ComplexMatrix& a, bool upper, bool calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
766 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
767 octave_idx_type a_nr = a.rows ();
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
768 octave_idx_type a_nc = a.cols ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
769
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
770 if (a_nr != a_nc)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
771 (*current_liboctave_error_handler) ("chol: requires square matrix");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
772
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
773 F77_INT n = to_f77_int (a_nc);
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
774 F77_INT info;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
775
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
776 is_upper = upper;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
777
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
778 chol_mat.clear (n, n);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
779 if (is_upper)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
780 for (octave_idx_type j = 0; j < n; j++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
781 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
782 for (octave_idx_type i = 0; i <= j; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
783 chol_mat.xelem (i, j) = a(i, j);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
784 for (octave_idx_type i = j+1; i < n; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
785 chol_mat.xelem (i, j) = 0.0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
786 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
787 else
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
788 for (octave_idx_type j = 0; j < n; j++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
789 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
790 for (octave_idx_type i = 0; i < j; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
791 chol_mat.xelem (i, j) = 0.0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
792 for (octave_idx_type i = j; i < n; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
793 chol_mat.xelem (i, j) = a(i, j);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
794 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
795 Complex *h = chol_mat.fortran_vec ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
796
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
797 // Calculate the norm of the matrix, for later use.
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
798 double anorm = 0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
799 if (calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
800 anorm = xnorm (a, 1);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
801
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
802 if (is_upper)
22402
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
803 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n,
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
804 F77_DBLE_CMPLX_ARG (h), n, info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
805 F77_CHAR_ARG_LEN (1)));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
806 else
22402
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
807 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1), n,
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
808 F77_DBLE_CMPLX_ARG (h), n, info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
809 F77_CHAR_ARG_LEN (1)));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
810
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
811 xrcond = 0.0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
812 if (info > 0)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
813 chol_mat.resize (info - 1, info - 1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
814 else if (calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
815 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
816 F77_INT zpocon_info = 0;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
817
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
818 // Now calculate the condition number for non-singular matrix.
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
819 Array<Complex> z (dim_vector (2*n, 1));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
820 Complex *pz = z.fortran_vec ();
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
821 Array<double> rz (dim_vector (n, 1));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
822 double *prz = rz.fortran_vec ();
22402
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
823 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n,
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
824 F77_DBLE_CMPLX_ARG (h), n, anorm, xrcond,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
825 F77_DBLE_CMPLX_ARG (pz), prz, zpocon_info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
826 F77_CHAR_ARG_LEN (1)));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
827
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
828 if (zpocon_info != 0)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
829 info = -1;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
830 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
831
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
832 return info;
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
833 }
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
834
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21301
diff changeset
835 #if defined (HAVE_QRUPDATE)
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
836
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
837 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
838 void
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
839 chol<ComplexMatrix>::update (const ComplexColumnVector& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
840 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
841 F77_INT n = to_f77_int (chol_mat.rows ());
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
842
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
843 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
844 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
845
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
846 ComplexColumnVector utmp = u;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
847
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
848 OCTAVE_LOCAL_BUFFER (double, rw, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
849
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
850 F77_XFCN (zch1up, ZCH1UP, (n,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
851 F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
852 n,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
853 F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
854 rw));
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
855 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
856
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
857 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
858 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
859 chol<ComplexMatrix>::downdate (const ComplexColumnVector& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
860 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
861 F77_INT info = -1;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
862
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
863 F77_INT n = to_f77_int (chol_mat.rows ());
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
864
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
865 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
866 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
867
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
868 ComplexColumnVector utmp = u;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
869
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
870 OCTAVE_LOCAL_BUFFER (double, rw, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
871
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
872 F77_XFCN (zch1dn, ZCH1DN, (n,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
873 F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
874 n,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
875 F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
876 rw, info));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
877
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
878 return info;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
879 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
880
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
881 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
882 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
883 chol<ComplexMatrix>::insert_sym (const ComplexColumnVector& u,
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
884 octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
885 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
886 F77_INT info = -1;
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
887
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
888 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
889 F77_INT j = to_f77_int (j_arg);
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
890
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
891 if (u.numel () != n + 1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
892 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
893 if (j < 0 || j > n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
894 (*current_liboctave_error_handler) ("cholinsert: index out of range");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
895
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
896 ComplexColumnVector utmp = u;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
897
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
898 OCTAVE_LOCAL_BUFFER (double, rw, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
899
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
900 chol_mat.resize (n+1, n+1);
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
901 F77_INT ldcm = to_f77_int (chol_mat.rows ());
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
902
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
903 F77_XFCN (zchinx, ZCHINX, (n,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
904 F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
22969
023026207910 * chol.cc: Fix thinko in previous change.
John W. Eaton <jwe@octave.org>
parents: 22942
diff changeset
905 ldcm, j + 1,
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
906 F77_DBLE_CMPLX_ARG (utmp.fortran_vec ()),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
907 rw, info));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
908
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
909 return info;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
910 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
911
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
912 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
913 void
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
914 chol<ComplexMatrix>::delete_sym (octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
915 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
916 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
917 F77_INT j = to_f77_int (j_arg);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
918
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
919 if (j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
920 (*current_liboctave_error_handler) ("choldelete: index out of range");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
921
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
922 OCTAVE_LOCAL_BUFFER (double, rw, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
923
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
924 F77_XFCN (zchdex, ZCHDEX, (n,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
925 F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
926 n, j + 1, rw));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
927
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
928 chol_mat.resize (n-1, n-1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
929 }
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
930
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
931 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
932 void
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
933 chol<ComplexMatrix>::shift_sym (octave_idx_type i_arg,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
934 octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
935 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
936 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
937 F77_INT i = to_f77_int (i_arg);
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
938 F77_INT j = to_f77_int (j_arg);
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
939
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
940 if (i < 0 || i > n-1 || j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
941 (*current_liboctave_error_handler) ("cholshift: index out of range");
21136
7cac4e7458f2 maint: clean up code around calls to current_liboctave_error_handler.
Rik <rik@octave.org>
parents: 20629
diff changeset
942
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
943 OCTAVE_LOCAL_BUFFER (Complex, w, n);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
944 OCTAVE_LOCAL_BUFFER (double, rw, n);
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
945
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
946 F77_XFCN (zchshx, ZCHSHX, (n,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
947 F77_DBLE_CMPLX_ARG (chol_mat.fortran_vec ()),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
948 n, i + 1, j + 1,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
949 F77_DBLE_CMPLX_ARG (w), rw));
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
950 }
8562
a6edd5c23cb5 use replacement methods if qrupdate is not available
Jaroslav Hajek <highegg@gmail.com>
parents: 8547
diff changeset
951
8547
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
952 #endif
d66c9b6e506a imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents: 8377
diff changeset
953
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
954 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
955 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
956 chol<FloatComplexMatrix>::init (const FloatComplexMatrix& a, bool upper,
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
957 bool calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
958 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
959 octave_idx_type a_nr = a.rows ();
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
960 octave_idx_type a_nc = a.cols ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
961
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
962 if (a_nr != a_nc)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
963 (*current_liboctave_error_handler) ("chol: requires square matrix");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
964
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
965 F77_INT n = to_f77_int (a_nc);
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
966 F77_INT info;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
967
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
968 is_upper = upper;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
969
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
970 chol_mat.clear (n, n);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
971 if (is_upper)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
972 for (octave_idx_type j = 0; j < n; j++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
973 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
974 for (octave_idx_type i = 0; i <= j; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
975 chol_mat.xelem (i, j) = a(i, j);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
976 for (octave_idx_type i = j+1; i < n; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
977 chol_mat.xelem (i, j) = 0.0f;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
978 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
979 else
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
980 for (octave_idx_type j = 0; j < n; j++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
981 {
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
982 for (octave_idx_type i = 0; i < j; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
983 chol_mat.xelem (i, j) = 0.0f;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
984 for (octave_idx_type i = j; i < n; i++)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
985 chol_mat.xelem (i, j) = a(i, j);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
986 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
987 FloatComplex *h = chol_mat.fortran_vec ();
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
988
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
989 // Calculate the norm of the matrix, for later use.
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
990 float anorm = 0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
991 if (calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
992 anorm = xnorm (a, 1);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
993
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
994 if (is_upper)
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
995 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
996 n, F77_CMPLX_ARG (h), n, info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
997 F77_CHAR_ARG_LEN (1)));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
998 else
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
999 F77_XFCN (cpotrf, CPOTRF, (F77_CONST_CHAR_ARG2 ("L", 1),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1000 n, F77_CMPLX_ARG (h), n, info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1001 F77_CHAR_ARG_LEN (1)));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1002
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1003 xrcond = 0.0;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1004 if (info > 0)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1005 chol_mat.resize (info - 1, info - 1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1006 else if (calc_cond)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1007 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1008 F77_INT cpocon_info = 0;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1009
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1010 // Now calculate the condition number for non-singular matrix.
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1011 Array<FloatComplex> z (dim_vector (2*n, 1));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1012 FloatComplex *pz = z.fortran_vec ();
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1013 Array<float> rz (dim_vector (n, 1));
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1014 float *prz = rz.fortran_vec ();
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1015 F77_XFCN (cpocon, CPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1016 F77_CMPLX_ARG (h), n, anorm, xrcond,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1017 F77_CMPLX_ARG (pz), prz, cpocon_info
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1018 F77_CHAR_ARG_LEN (1)));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1019
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1020 if (cpocon_info != 0)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1021 info = -1;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1022 }
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1023
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1024 return info;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1025 }
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1026
21724
aba2e6293dd8 use "#if ..." consistently instead of "#ifdef" and "#ifndef"
John W. Eaton <jwe@octave.org>
parents: 21301
diff changeset
1027 #if defined (HAVE_QRUPDATE)
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1028
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1029 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1030 void
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1031 chol<FloatComplexMatrix>::update (const FloatComplexColumnVector& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1032 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1033 F77_INT n = to_f77_int (chol_mat.rows ());
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1034
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1035 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1036 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1037
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1038 FloatComplexColumnVector utmp = u;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1039
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1040 OCTAVE_LOCAL_BUFFER (float, rw, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1041
22402
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
1042 F77_XFCN (cch1up, CCH1UP, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1043 n, F77_CMPLX_ARG (utmp.fortran_vec ()), rw));
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1044 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1045
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1046 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1047 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1048 chol<FloatComplexMatrix>::downdate (const FloatComplexColumnVector& u)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1049 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1050 F77_INT info = -1;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1051
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1052 F77_INT n = to_f77_int (chol_mat.rows ());
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1053
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1054 if (u.numel () != n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1055 (*current_liboctave_error_handler) ("cholupdate: dimension mismatch");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1056
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1057 FloatComplexColumnVector utmp = u;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1058
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1059 OCTAVE_LOCAL_BUFFER (float, rw, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1060
22402
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
1061 F77_XFCN (cch1dn, CCH1DN, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1062 n, F77_CMPLX_ARG (utmp.fortran_vec ()),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1063 rw, info));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1064
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1065 return info;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1066 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1067
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1068 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1069 octave_idx_type
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1070 chol<FloatComplexMatrix>::insert_sym (const FloatComplexColumnVector& u,
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1071 octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1072 {
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1073 F77_INT info = -1;
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1074 F77_INT j = to_f77_int (j_arg);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1075
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1076 F77_INT n = to_f77_int (chol_mat.rows ());
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1077
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1078 if (u.numel () != n + 1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1079 (*current_liboctave_error_handler) ("cholinsert: dimension mismatch");
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1080 if (j < 0 || j > n)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1081 (*current_liboctave_error_handler) ("cholinsert: index out of range");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1082
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1083 FloatComplexColumnVector utmp = u;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1084
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1085 OCTAVE_LOCAL_BUFFER (float, rw, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1086
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1087 chol_mat.resize (n+1, n+1);
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1088 F77_INT ldcm = to_f77_int (chol_mat.rows ());
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1089
22402
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
1090 F77_XFCN (cchinx, CCHINX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
22969
023026207910 * chol.cc: Fix thinko in previous change.
John W. Eaton <jwe@octave.org>
parents: 22942
diff changeset
1091 ldcm, j + 1,
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1092 F77_CMPLX_ARG (utmp.fortran_vec ()),
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1093 rw, info));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1094
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1095 return info;
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1096 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1097
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1098 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1099 void
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1100 chol<FloatComplexMatrix>::delete_sym (octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1101 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1102 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1103 F77_INT j = to_f77_int (j_arg);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1104
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1105 if (j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1106 (*current_liboctave_error_handler) ("choldelete: index out of range");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1107
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1108 OCTAVE_LOCAL_BUFFER (float, rw, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1109
22402
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
1110 F77_XFCN (cchdex, CCHDEX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1111 n, j + 1, rw));
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1112
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1113 chol_mat.resize (n-1, n-1);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1114 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1115
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1116 template <>
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1117 void
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1118 chol<FloatComplexMatrix>::shift_sym (octave_idx_type i_arg,
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1119 octave_idx_type j_arg)
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1120 {
23696
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1121 F77_INT n = to_f77_int (chol_mat.rows ());
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1122 F77_INT i = to_f77_int (i_arg);
08036a7f3660 remove octave:: namespace tag from symbols used inside octave namespace
John W. Eaton <jwe@octave.org>
parents: 23593
diff changeset
1123 F77_INT j = to_f77_int (j_arg);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1124
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1125 if (i < 0 || i > n-1 || j < 0 || j > n-1)
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1126 (*current_liboctave_error_handler) ("cholshift: index out of range");
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1127
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1128 OCTAVE_LOCAL_BUFFER (FloatComplex, w, n);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1129 OCTAVE_LOCAL_BUFFER (float, rw, n);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1130
22402
4caa7b28d183 maint: Style check C++ code in liboctave/
Rik <rik@octave.org>
parents: 22329
diff changeset
1131 F77_XFCN (cchshx, CCHSHX, (n, F77_CMPLX_ARG (chol_mat.fortran_vec ()),
22942
7627ebcf194d use F77_INT instead of octave_idx_type for liboctave chol classes
John W. Eaton <jwe@octave.org>
parents: 22927
diff changeset
1132 n, i + 1, j + 1, F77_CMPLX_ARG (w), rw));
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1133 }
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1134
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1135 #endif
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1136
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1137 // Instantiations we need.
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1138
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1139 template class chol<Matrix>;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1140
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1141 template class chol<FloatMatrix>;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1142
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1143 template class chol<ComplexMatrix>;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1144
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1145 template class chol<FloatComplexMatrix>;
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1146
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1147 template Matrix
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1148 chol2inv<Matrix> (const Matrix& r);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1149
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1150 template ComplexMatrix
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1151 chol2inv<ComplexMatrix> (const ComplexMatrix& r);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1152
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1153 template FloatMatrix
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1154 chol2inv<FloatMatrix> (const FloatMatrix& r);
21269
3c8a3d35661a better use of templates for Cholesky factorization
John W. Eaton <jwe@octave.org>
parents: 21202
diff changeset
1155
22329
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1156 template FloatComplexMatrix
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1157 chol2inv<FloatComplexMatrix> (const FloatComplexMatrix& r);
7f3c7a8bd131 maint: Indent namespaces in liboctave/numeric files.
John W. Eaton <jwe@octave.org>
parents: 22323
diff changeset
1158 }
22317
6ca3acf5fad8 move some new numeric classes to namespace octave::math
John W. Eaton <jwe@octave.org>
parents: 22274
diff changeset
1159 }