Mercurial > octave

annotate liboctave/dbleCHOL.cc @ 8377:25bc2d31e1bf

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
improve OCTAVE_LOCAL_BUFFER
author Jaroslav Hajek <highegg@gmail.com>
date Wed, 29 Oct 2008 16:52:10 +0100
parents 7c9ba697a479
children d66c9b6e506a
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"
8377
25bc2d31e1bf improve OCTAVE_LOCAL_BUFFER
Jaroslav Hajek <highegg@gmail.com>
parents: 7725
diff changeset
36 #include "oct-locbuf.h"
457
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
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
41 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
42 double*, const octave_idx_type&, octave_idx_type&
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
43 F77_CHAR_ARG_LEN_DECL);
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
44
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
45 F77_RET_T
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
46 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
47 double*, const octave_idx_type&, octave_idx_type&
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
48 F77_CHAR_ARG_LEN_DECL);
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
49
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
50 F77_RET_T
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
51 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
52 double*, const octave_idx_type&, const double&,
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
53 double&, double*, octave_idx_type*,
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
54 octave_idx_type& F77_CHAR_ARG_LEN_DECL);
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
55 F77_RET_T
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
56 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
57
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
58 F77_RET_T
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
59 F77_FUNC (dch1dn, DCH1DN) (const octave_idx_type&, double*, double*, double*,
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
60 octave_idx_type&);
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
61
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
62 F77_RET_T
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
63 F77_FUNC (dqrshc, DQRSHC) (const octave_idx_type&, const octave_idx_type&, const octave_idx_type&,
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
64 double*, double*, const octave_idx_type&, const octave_idx_type&);
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
65
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
66 F77_RET_T
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
67 F77_FUNC (dchinx, DCHINX) (const octave_idx_type&, const double*, double*, const octave_idx_type&,
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
68 const double*, octave_idx_type&);
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
69
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
70 F77_RET_T
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
71 F77_FUNC (dchdex, DCHDEX) (const octave_idx_type&, const double*, double*, const octave_idx_type&);
457
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
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
75 CHOL::init (const Matrix& a, bool calc_cond)
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
76 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
77 octave_idx_type a_nr = a.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
78 octave_idx_type a_nc = a.cols ();
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
79
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
80 if (a_nr != a_nc)
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
81 {
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
82 (*current_liboctave_error_handler) ("CHOL requires square matrix");
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
83 return -1;
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
84 }
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
85
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
86 octave_idx_type n = a_nc;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5105
diff changeset
87 octave_idx_type info;
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
88
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
89 chol_mat = a;
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
90 double *h = chol_mat.fortran_vec ();
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
91
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
92 // Calculate the norm of the matrix, for later use.
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
93 double anorm = 0;
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
94 if (calc_cond)
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
95 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
96
4552
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
97 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
98 n, h, n, info
6f3382e08a52 [project @ 2003-10-27 20:38:02 by jwe]
jwe
parents: 4192
diff changeset
99 F77_CHAR_ARG_LEN (1)));
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
100
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
101 xrcond = 0.0;
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
102 if (info != 0)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
103 info = -1;
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
104 else if (calc_cond)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
105 {
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
106 octave_idx_type dpocon_info = 0;
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
107
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
108 // 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
109 Array<double> z (3*n);
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
110 double *pz = z.fortran_vec ();
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
111 Array<octave_idx_type> iz (n);
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
112 octave_idx_type *piz = iz.fortran_vec ();
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
113 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
114 n, anorm, xrcond, pz, piz, dpocon_info
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
115 F77_CHAR_ARG_LEN (1)));
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
116
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
117 if (dpocon_info != 0)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
118 info = -1;
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
119 }
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
120 else
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
121 {
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
122 // 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
123 // faster for that matter :-)), please let me know!
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
124
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
125 if (n > 1)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
126 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
127 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
128 chol_mat.xelem (i, j) = 0.0;
1944
8cb4d3008c76 [project @ 1996-02-14 00:45:04 by jwe]
jwe
parents: 1882
diff changeset
129 }
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
130
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
131 return info;
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
132 }
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
133
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
134 static Matrix
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
135 chol2inv_internal (const Matrix& 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 Matrix retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
138
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
139 octave_idx_type r_nr = r.rows ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
140 octave_idx_type r_nc = r.cols ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
141
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
142 if (r_nr == r_nc)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
143 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
144 octave_idx_type n = r_nc;
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
145 octave_idx_type info = 0;
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
146
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
147 Matrix tmp = r;
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
148 double *v = tmp.fortran_vec();
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
149
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
150 if (info == 0)
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
151 {
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
152 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
153 v, n, info
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
154 F77_CHAR_ARG_LEN (1)));
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
155
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
156 // 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
157 // faster for that matter :-)), please let me know!
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
158
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
159 if (n > 1)
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
160 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
161 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
162 tmp.xelem (i, j) = tmp.xelem (j, i);
6486
e978a9233cf6 [project @ 2007-04-04 15:16:46 by jwe]
jwe
parents: 5340
diff changeset
163
7482
29980c6b8604 don't check f77_exception_encountered
John W. Eaton <jwe@octave.org>
parents: 7017
diff changeset
164 retval = tmp;
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
165 }
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 else
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
168 (*current_liboctave_error_handler) ("chol2inv requires square matrix");
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
169
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
170 return retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
171 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
172
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
173 // 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
174 Matrix
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
175 CHOL::inverse (void) const
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
176 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
177 return chol2inv_internal (chol_mat);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
178 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
179
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
180 void
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
181 CHOL::set (const Matrix& R)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
182 {
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
183 if (R.is_square ())
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
184 chol_mat = R;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
185 else
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
186 (*current_liboctave_error_handler) ("CHOL requires square matrix");
7554
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
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
189 void
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
190 CHOL::update (const Matrix& u)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
191 {
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
192 octave_idx_type n = chol_mat.rows ();
7554
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 if (u.length () == n)
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 Matrix tmp = 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_LOCAL_BUFFER (double, w, 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 F77_XFCN (dch1up, DCH1UP, (n, chol_mat.fortran_vec (),
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
201 tmp.fortran_vec (), w));
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
202 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
203 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
204 (*current_liboctave_error_handler) ("CHOL update dimension mismatch");
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
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
207 octave_idx_type
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
208 CHOL::downdate (const Matrix& u)
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 octave_idx_type info = -1;
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 octave_idx_type n = chol_mat.rows ();
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 if (u.length () == n)
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 Matrix tmp = u;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
217
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
218 OCTAVE_LOCAL_BUFFER (double, w, n);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
219
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
220 F77_XFCN (dch1dn, DCH1DN, (n, chol_mat.fortran_vec (),
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
221 tmp.fortran_vec (), w, info));
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
222 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
223 else
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
224 (*current_liboctave_error_handler) ("CHOL downdate dimension mismatch");
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
225
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
226 return info;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
227 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7482
diff changeset
228
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
229 octave_idx_type
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
230 CHOL::insert_sym (const Matrix& u, octave_idx_type j)
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
231 {
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
232 octave_idx_type info = -1;
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
233
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
234 octave_idx_type n = chol_mat.rows ();
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
235
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
236 if (u.length () != n+1)
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
237 (*current_liboctave_error_handler) ("CHOL insert dimension mismatch");
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
238 else if (j < 0 || j > n)
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
239 (*current_liboctave_error_handler) ("CHOL insert index out of range");
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
240 else
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
241 {
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
242 Matrix chol_mat1 (n+1, n+1);
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
243
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
244 F77_XFCN (dchinx, DCHINX, (n, chol_mat.data (), chol_mat1.fortran_vec (),
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
245 j+1, u.data (), info));
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
246
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
247 chol_mat = chol_mat1;
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
248 }
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
249
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
250 return info;
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
251 }
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
252
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
253 void
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
254 CHOL::delete_sym (octave_idx_type j)
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
255 {
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
256 octave_idx_type n = chol_mat.rows ();
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
257
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
258 if (j < 0 || j > n-1)
7725
7c9ba697a479 cosmetic fixes in QR & Cholesky updating codes
Jaroslav Hajek <highegg@gmail.com>
parents: 7700
diff changeset
259 (*current_liboctave_error_handler) ("CHOL delete index out of range");
7700
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
260 else
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
261 {
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
262 Matrix chol_mat1 (n-1, n-1);
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
263
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
264 F77_XFCN (dchdex, DCHDEX, (n, chol_mat.data (), chol_mat1.fortran_vec (), j+1));
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
265
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
266 chol_mat = chol_mat1;
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
267 }
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
268 }
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
269
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
270 void
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
271 CHOL::shift_sym (octave_idx_type i, octave_idx_type j)
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
272 {
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
273 octave_idx_type n = chol_mat.rows ();
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
274 double dummy;
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
275
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
276 if (i < 0 || i > n-1 || j < 0 || j > n-1)
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
277 (*current_liboctave_error_handler) ("CHOL shift index out of range");
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
278 else
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
279 F77_XFCN (dqrshc, DQRSHC, (0, n, n, &dummy, chol_mat.fortran_vec (), i+1, j+1));
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
280 }
efccca5f2ad7 more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents: 7559
diff changeset
281
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
282 Matrix
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
283 chol2inv (const Matrix& r)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
284 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
285 return chol2inv_internal (r);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
286 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
287
457
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
288 /*
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
289 ;;; Local Variables: ***
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
290 ;;; mode: C++ ***
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
291 ;;; End: ***
3d4b4f0fa5ba [project @ 1994-06-06 00:33:33 by jwe]
jwe
parents:
diff changeset
292 */