Mercurial > hg > octave-lyh

annotate liboctave/CmplxCHOL.cc @ 7554:40574114c514

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
implement Cholesky factorization updating
author Jaroslav Hajek <highegg@gmail.com>
date Tue, 04 Mar 2008 22:25:50 -0500
parents 29980c6b8604
children 07522d7dcdf8
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
7017
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
3 Copyright (C) 1994, 1995, 1996, 1997, 2002, 2003, 2004, 2005, 2007
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
4 John W. Eaton
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
5
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
6 This file is part of Octave.
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
7
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
8 Octave is free software; you can redistribute it and/or modify it
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
9 under the terms of the GNU General Public License as published by the
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6486
diff changeset
10 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6486
diff changeset
11 option) any later version.
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
12
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
13 Octave is distributed in the hope that it will be useful, but WITHOUT
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
16 for more details.
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
17
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
18 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
19 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6486
diff changeset
20 <http://www.gnu.org/licenses/>.
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
21
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
22 */
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
23
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
24 // updating/downdating by Jaroslav Hajek 2008
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
25
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
26 #ifdef HAVE_CONFIG_H
1192
b6360f2d4fa6 [project @ 1995-03-30 21:38:35 by jwe]
jwe
parents: 1011
diff changeset
27 #include <config.h>
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
28 #endif
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
29
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
30 #include <vector>
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
31
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
32 #include "dMatrix.h"
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
33 #include "dRowVector.h"
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
34 #include "CmplxCHOL.h"
1847
2ffe49eb95a5 [project @ 1996-02-03 12:47:55 by jwe]
jwe
parents: 1368
diff changeset
35 #include "f77-fcn.h"
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
36 #include "lo-error.h"
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
37
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
38 extern "C"
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
39 {
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
40 F77_RET_T
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
41 F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
42 Complex*, const octave_idx_type&, octave_idx_type&
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
43 F77_CHAR_ARG_LEN_DECL);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
44 F77_RET_T
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
45 F77_FUNC (zpotri, ZPOTRI) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
46 Complex*, const octave_idx_type&, octave_idx_type&
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
47 F77_CHAR_ARG_LEN_DECL);
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
48
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
49 F77_RET_T
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
50 F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
51 Complex*, const octave_idx_type&, const double&,
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
52 double&, Complex*, double*,
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
53 octave_idx_type& F77_CHAR_ARG_LEN_DECL);
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
54 F77_RET_T
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
55 F77_FUNC (zch1up, ZCH1UP) (const octave_idx_type&, Complex*, Complex*, double*);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
56
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
57 F77_RET_T
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
58 F77_FUNC (zch1dn, ZCH1DN) (const octave_idx_type&, Complex*, Complex*, double*,
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
59 octave_idx_type&);
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
60 }
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
61
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
62 octave_idx_type
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
63 ComplexCHOL::init (const ComplexMatrix& a, bool calc_cond)
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
64 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
65 octave_idx_type a_nr = a.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
66 octave_idx_type a_nc = a.cols ();
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
67
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
68 if (a_nr != a_nc)
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
69 {
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
70 (*current_liboctave_error_handler)
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
71 ("ComplexCHOL requires square matrix");
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
72 return -1;
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
73 }
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
74
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
75 octave_idx_type n = a_nc;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
76 octave_idx_type info;
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
77
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
78 chol_mat = a;
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
79 Complex *h = chol_mat.fortran_vec ();
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
80
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
81 // Calculate the norm of the matrix, for later use.
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
82 double anorm = 0;
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
83 if (calc_cond)
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
84 anorm = chol_mat.abs().sum().row(static_cast<octave_idx_type>(0)).max();
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
85
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
86 F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1), n, h, n, info
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
87 F77_CHAR_ARG_LEN (1)));
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
88
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
89 xrcond = 0.0;
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
90 if (info != 0)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
91 info = -1;
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
92 else if (calc_cond)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
93 {
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
94 octave_idx_type zpocon_info = 0;
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
95
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
96 // Now calculate the condition number for non-singular matrix.
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
97 Array<Complex> z (2*n);
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
98 Complex *pz = z.fortran_vec ();
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
99 Array<double> rz (n);
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
100 double *prz = rz.fortran_vec ();
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
101 F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 ("U", 1), n, h,
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
102 n, anorm, xrcond, pz, prz, zpocon_info
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
103 F77_CHAR_ARG_LEN (1)));
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
104
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
105 if (zpocon_info != 0)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
106 info = -1;
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
107 }
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
108 else
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
109 {
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
110 // If someone thinks of a more graceful way of doing this (or
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
111 // faster for that matter :-)), please let me know!
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
112
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
113 if (n > 1)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
114 for (octave_idx_type j = 0; j < a_nc; j++)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
115 for (octave_idx_type i = j+1; i < a_nr; i++)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
116 chol_mat.xelem (i, j) = 0.0;
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
117 }
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
118
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
119 return info;
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
120 }
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
121
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
122 static ComplexMatrix
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
123 chol2inv_internal (const ComplexMatrix& r)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
124 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
125 ComplexMatrix retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
126
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
127 octave_idx_type r_nr = r.rows ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
128 octave_idx_type r_nc = r.cols ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
129
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
130 if (r_nr == r_nc)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
131 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
132 octave_idx_type n = r_nc;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
133 octave_idx_type info;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
134
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
135 ComplexMatrix tmp = r;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
136
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
137 F77_XFCN (zpotri, ZPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
138 tmp.fortran_vec (), n, info
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
139 F77_CHAR_ARG_LEN (1)));
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
140
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
141 // If someone thinks of a more graceful way of doing this (or
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
142 // faster for that matter :-)), please let me know!
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
143
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
144 if (n > 1)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
145 for (octave_idx_type j = 0; j < r_nc; j++)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
146 for (octave_idx_type i = j+1; i < r_nr; i++)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
147 tmp.xelem (i, j) = std::conj (tmp.xelem (j, i));
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
148
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
149 retval = tmp;
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
150 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
151 else
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
152 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
153
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
154 return retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
155 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
156
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
157 // Compute the inverse of a matrix using the Cholesky factorization.
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
158 ComplexMatrix
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
159 ComplexCHOL::inverse (void) const
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
160 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
161 return chol2inv_internal (chol_mat);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
162 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
163
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
164 void
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
165 ComplexCHOL::set (const ComplexMatrix& R)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
166 {
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
167 if (R.is_square ())
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
168 chol_mat = R;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
169 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
170 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
171 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
172
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
173 void
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
174 ComplexCHOL::update (const ComplexMatrix& u)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
175 {
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
176 octave_idx_type n = chol_mat.rows ();
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
177
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
178 if (u.length () == n)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
179 {
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
180 ComplexMatrix tmp = u;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
181
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
182 OCTAVE_LOCAL_BUFFER (double, w, n);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
183
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
184 F77_XFCN (zch1up, ZCH1UP, (n, chol_mat.fortran_vec (),
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
185 tmp.fortran_vec (), w));
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
186 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
187 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
188 (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
189 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
190
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
191 octave_idx_type
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
192 ComplexCHOL::downdate (const ComplexMatrix& u)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
193 {
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
194 octave_idx_type info = -1;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
195
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
196 octave_idx_type n = chol_mat.rows ();
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
197
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
198 if (u.length () == n)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
199 {
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
200 ComplexMatrix tmp = u;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
201
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
202 OCTAVE_LOCAL_BUFFER (double, w, n);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
203
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
204 F77_XFCN (zch1dn, ZCH1DN, (n, chol_mat.fortran_vec (),
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
205 tmp.fortran_vec (), w, info));
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
206 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
207 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
208 (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
209
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
210 return info;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
211 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
212
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
213 ComplexMatrix
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
214 chol2inv (const ComplexMatrix& r)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
215 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
216 return chol2inv_internal (r);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
217 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
218
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
219 /*
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
220 ;;; Local Variables: ***
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
221 ;;; mode: C++ ***
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
222 ;;; End: ***
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
223 */