Mercurial > octave

annotate liboctave/dbleCHOL.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 "dRowVector.h"
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
33 #include "dbleCHOL.h"
1847
2ffe49eb95a5 [project @ 1996-02-03 12:47:55 by jwe]
jwe
parents: 1368
diff changeset
34 #include "f77-fcn.h"
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
35 #include "lo-error.h"
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
36
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
37 extern "C"
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
38 {
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
39 F77_RET_T
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
40 F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
41 double*, const octave_idx_type&, octave_idx_type&
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
42 F77_CHAR_ARG_LEN_DECL);
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
43
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 (dpotri, DPOTRI) (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 double*, const octave_idx_type&, octave_idx_type&
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
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 (dpocon, DPOCON) (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 double*, const octave_idx_type&, const double&,
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
52 double&, double*, octave_idx_type*,
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 (dch1up, DCH1UP) (const octave_idx_type&, double*, double*, 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 (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*,
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
59 int &);
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 CHOL::init (const Matrix& 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 ();
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
67
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
68 if (a_nr != a_nc)
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
69 {
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
70 (*current_liboctave_error_handler) ("CHOL requires square matrix");
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
71 return -1;
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
72 }
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
73
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
74 octave_idx_type n = a_nc;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
75 octave_idx_type info;
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
76
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
77 chol_mat = a;
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
78 double *h = chol_mat.fortran_vec ();
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
79
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
80 // Calculate the norm of the matrix, for later use.
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
81 double anorm = 0;
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
82 if (calc_cond)
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
83 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
84
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
85 F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 ("U", 1),
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
86 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 dpocon_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<double> z (3*n);
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
98 double *pz = z.fortran_vec ();
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
99 Array<octave_idx_type> iz (n);
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
100 octave_idx_type *piz = iz.fortran_vec ();
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
101 F77_XFCN (dpocon, DPOCON, (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, piz, dpocon_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 (dpocon_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 }
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
118
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
119 return info;
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 Matrix
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
123 chol2inv_internal (const Matrix& 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 Matrix 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;
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
133 octave_idx_type info = 0;
5340
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 Matrix tmp = r;
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
136 double *v = tmp.fortran_vec();
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
137
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
138 if (info == 0)
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
139 {
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
140 F77_XFCN (dpotri, DPOTRI, (F77_CONST_CHAR_ARG2 ("U", 1), n,
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
141 v, n, info
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
142 F77_CHAR_ARG_LEN (1)));
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 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
145 // faster for that matter :-)), please let me know!
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
146
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
147 if (n > 1)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
148 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
149 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
150 tmp.xelem (i, j) = tmp.xelem (j, i);
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
151
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
152 retval = tmp;
5340
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 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
155 else
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
156 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
157
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
158 return retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
159 }
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 // 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
162 Matrix
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
163 CHOL::inverse (void) const
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
164 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
165 return chol2inv_internal (chol_mat);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
166 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
167
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
168 void
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
169 CHOL::set (const Matrix& R)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
170 {
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
171 if (R.is_square ())
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
172 chol_mat = R;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
173 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
174 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
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
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
177 void
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
178 CHOL::update (const Matrix& u)
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 int n = chol_mat.rows ();
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 if (u.length () == 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 Matrix tmp = u;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
185
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
186 OCTAVE_LOCAL_BUFFER (double, w, n);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
187
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
188 F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (),
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
189 tmp.fortran_vec (), w));
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 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
192 (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
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
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
195 octave_idx_type
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
196 CHOL::downdate (const Matrix& u)
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 octave_idx_type info = -1;
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 octave_idx_type n = chol_mat.rows ();
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 if (u.length () == 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 Matrix tmp = u;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
205
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
206 OCTAVE_LOCAL_BUFFER (double, w, n);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
207
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
208 F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (),
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
209 tmp.fortran_vec (), w, info));
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
210 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
211 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
212 (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
213
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
214 return info;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
215 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
216
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
217 Matrix
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
218 chol2inv (const Matrix& r)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
219 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
220 return chol2inv_internal (r);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
221 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
222
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
223 /*
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
224 ;;; Local Variables: ***
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
225 ;;; mode: C++ ***
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
226 ;;; End: ***
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
227 */