Mercurial > octave-nkf
annotate src/DLD-FUNCTIONS/chol.cc @ 10154:40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 20 Jan 2010 17:33:41 -0500 |
parents | c0aeedd8fb86 |
children | d0ce5e973937 |
rev | line source |
---|---|
2928 | 1 /* |
2 | |
7017 | 3 Copyright (C) 1996, 1997, 1999, 2000, 2002, 2005, 2006, 2007 |
4 John W. Eaton | |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
5 Copyright (C) 2008, 2009 Jaroslav Hajek |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
6 Copyright (C) 2008, 2009 VZLU Prague |
2928 | 7 |
8 This file is part of Octave. | |
9 | |
10 Octave is free software; you can redistribute it and/or modify it | |
11 under the terms of the GNU General Public License as published by the | |
7016 | 12 Free Software Foundation; either version 3 of the License, or (at your |
13 option) any later version. | |
2928 | 14 |
15 Octave is distributed in the hope that it will be useful, but WITHOUT | |
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
17 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
18 for more details. | |
19 | |
20 You should have received a copy of the GNU General Public License | |
7016 | 21 along with Octave; see the file COPYING. If not, see |
22 <http://www.gnu.org/licenses/>. | |
2928 | 23 |
24 */ | |
25 | |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
26 |
2928 | 27 #ifdef HAVE_CONFIG_H |
28 #include <config.h> | |
29 #endif | |
30 | |
31 #include "CmplxCHOL.h" | |
32 #include "dbleCHOL.h" | |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7700
diff
changeset
|
33 #include "fCmplxCHOL.h" |
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7700
diff
changeset
|
34 #include "floatCHOL.h" |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
35 #include "SparseCmplxCHOL.h" |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
36 #include "SparsedbleCHOL.h" |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
37 #include "oct-spparms.h" |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
38 #include "sparse-util.h" |
2928 | 39 |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
40 #include "ov-re-sparse.h" |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
41 #include "ov-cx-sparse.h" |
2928 | 42 #include "defun-dld.h" |
43 #include "error.h" | |
44 #include "gripes.h" | |
45 #include "oct-obj.h" | |
46 #include "utils.h" | |
47 | |
9715
9f27172fbd1e
auto-set MatrixType from certain functions
Jaroslav Hajek <highegg@gmail.com>
parents:
9453
diff
changeset
|
48 template <class CHOLT> |
9f27172fbd1e
auto-set MatrixType from certain functions
Jaroslav Hajek <highegg@gmail.com>
parents:
9453
diff
changeset
|
49 static octave_value |
9f27172fbd1e
auto-set MatrixType from certain functions
Jaroslav Hajek <highegg@gmail.com>
parents:
9453
diff
changeset
|
50 get_chol_r (const CHOLT& fact) |
9f27172fbd1e
auto-set MatrixType from certain functions
Jaroslav Hajek <highegg@gmail.com>
parents:
9453
diff
changeset
|
51 { |
9f27172fbd1e
auto-set MatrixType from certain functions
Jaroslav Hajek <highegg@gmail.com>
parents:
9453
diff
changeset
|
52 return octave_value (fact.chol_matrix (), |
9f27172fbd1e
auto-set MatrixType from certain functions
Jaroslav Hajek <highegg@gmail.com>
parents:
9453
diff
changeset
|
53 MatrixType (MatrixType::Upper)); |
9f27172fbd1e
auto-set MatrixType from certain functions
Jaroslav Hajek <highegg@gmail.com>
parents:
9453
diff
changeset
|
54 } |
9f27172fbd1e
auto-set MatrixType from certain functions
Jaroslav Hajek <highegg@gmail.com>
parents:
9453
diff
changeset
|
55 |
9862
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
9758
diff
changeset
|
56 template <class CHOLT> |
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
9758
diff
changeset
|
57 static octave_value |
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
9758
diff
changeset
|
58 get_chol_l (const CHOLT& fact) |
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
9758
diff
changeset
|
59 { |
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
9758
diff
changeset
|
60 return octave_value (fact.chol_matrix ().transpose (), |
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
9758
diff
changeset
|
61 MatrixType (MatrixType::Lower)); |
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
9758
diff
changeset
|
62 } |
c0aeedd8fb86
improve chol Matlab compatibility
Jaroslav Hajek <highegg@gmail.com>
parents:
9758
diff
changeset
|
63 |
2928 | 64 DEFUN_DLD (chol, args, nargout, |
3548 | 65 "-*- texinfo -*-\n\ |
7650 | 66 @deftypefn {Loadable Function} {@var{r} =} chol (@var{a})\n\ |
67 @deftypefnx {Loadable Function} {[@var{r}, @var{p}] =} chol (@var{a})\n\ | |
68 @deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}] =} chol (@var{s})\n\ | |
69 @deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}] =} chol (@var{s}, 'vector')\n\ | |
70 @deftypefnx {Loadable Function} {[@var{l}, @dots{}] =} chol (@dots{}, 'lower')\n\ | |
3372 | 71 @cindex Cholesky factorization\n\ |
72 Compute the Cholesky factor, @var{r}, of the symmetric positive definite\n\ | |
73 matrix @var{a}, where\n\ | |
74 @tex\n\ | |
75 $ R^T R = A $.\n\ | |
76 @end tex\n\ | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8421
diff
changeset
|
77 @ifnottex\n\ |
3372 | 78 \n\ |
79 @example\n\ | |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
80 @var{r}' * @var{r} = @var{a}.\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
81 @end example\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8421
diff
changeset
|
82 @end ifnottex\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
83 \n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
84 Called with one output argument @code{chol} fails if @var{a} or @var{s} is\n\ |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8880
diff
changeset
|
85 not positive definite. With two or more output arguments @var{p} flags\n\ |
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8880
diff
changeset
|
86 whether the matrix was positive definite and @code{chol} does not fail. A\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
87 zero value indicated that the matrix was positive definite and the @var{r}\n\ |
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9064
diff
changeset
|
88 gives the factorization, and @var{p} will have a positive value otherwise.\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
89 \n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
90 If called with 3 outputs then a sparsity preserving row/column permutation\n\ |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8880
diff
changeset
|
91 is applied to @var{a} prior to the factorization. That is @var{r}\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
92 is the factorization of @code{@var{a}(@var{q},@var{q})} such that\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
93 @tex\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
94 $ R^T R = Q^T A Q$.\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
95 @end tex\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8421
diff
changeset
|
96 @ifnottex\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
97 \n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
98 @example\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
99 @var{r}' * @var{r} = @var{q}' * @var{a} * @var{q}.\n\ |
3372 | 100 @end example\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8421
diff
changeset
|
101 @end ifnottex\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
102 \n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
103 The sparsity preserving permutation is generally returned as a matrix.\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
104 However, given the flag 'vector', @var{q} will be returned as a vector\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
105 such that\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
106 @tex\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
107 $ R^T R = A (Q, Q)$.\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
108 @end tex\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8421
diff
changeset
|
109 @ifnottex\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
110 \n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
111 @example\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
112 @var{r}' * @var{r} = a (@var{q}, @var{q}).\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
113 @end example\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8421
diff
changeset
|
114 @end ifnottex\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
115 \n\ |
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9064
diff
changeset
|
116 Called with either a sparse or full matrix and using the 'lower' flag,\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
117 @code{chol} returns the lower triangular factorization such that\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
118 @tex\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
119 $ L L^T = A $.\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
120 @end tex\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8421
diff
changeset
|
121 @ifnottex\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
122 \n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
123 @example\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
124 @var{l} * @var{l}' = @var{a}.\n\ |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
125 @end example\n\ |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
8421
diff
changeset
|
126 @end ifnottex\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
127 \n\ |
9065
8207b833557f
Cleanup documentation for arith.texi, linalg.texi, nonlin.texi
Rik <rdrider0-list@yahoo.com>
parents:
9064
diff
changeset
|
128 In general the lower triangular factorization is significantly faster for\n\ |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
129 sparse matrices.\n\ |
5340 | 130 @seealso{cholinv, chol2inv}\n\ |
3372 | 131 @end deftypefn") |
2928 | 132 { |
133 octave_value_list retval; | |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
134 int nargin = args.length (); |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
135 bool LLt = false; |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
136 bool vecout = false; |
2928 | 137 |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
138 if (nargin < 1 || nargin > 3 || nargout > 3 |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
139 || (! args(0).is_sparse_type () && nargout > 2)) |
2928 | 140 { |
5823 | 141 print_usage (); |
2928 | 142 return retval; |
143 } | |
144 | |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
145 int n = 1; |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
146 while (n < nargin && ! error_state) |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
147 { |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
148 std::string tmp = args(n++).string_value (); |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
149 |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
150 if (! error_state ) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
151 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
152 if (tmp.compare ("vector") == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
153 vecout = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
154 else if (tmp.compare ("lower") == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
155 LLt = true; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
156 else if (tmp.compare ("upper") == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
157 LLt = false; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
158 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
159 error ("chol: unexpected second or third input"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
160 } |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
161 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
162 error ("chol: expecting trailing string arguments"); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
163 } |
2928 | 164 |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
165 if (! error_state) |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
166 { |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
167 octave_value arg = args(0); |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
168 |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
169 octave_idx_type nr = arg.rows (); |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
170 octave_idx_type nc = arg.columns (); |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
171 bool natural = (nargout != 3); |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
172 |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
173 int arg_is_empty = empty_arg ("chol", nr, nc); |
2928 | 174 |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
175 if (arg_is_empty < 0) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
176 return retval; |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
177 if (arg_is_empty > 0) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
178 return octave_value (Matrix ()); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
179 |
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
180 if (arg.is_sparse_type ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
181 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
182 if (arg.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
183 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
184 SparseMatrix m = arg.sparse_matrix_value (); |
2928 | 185 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
186 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
187 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
188 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
189 SparseCHOL fact (m, info, natural); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
190 if (nargout == 3) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
191 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
192 if (vecout) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
193 retval(2) = fact.perm (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
194 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
195 retval(2) = fact.Q(); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
196 } |
2928 | 197 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
198 if (nargout > 1 || info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
199 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
200 retval(1) = fact.P(); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
201 if (LLt) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
202 retval(0) = fact.L(); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
203 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
204 retval(0) = fact.R(); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
205 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
206 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
207 error ("chol: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
208 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
209 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
210 else if (arg.is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
211 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
212 SparseComplexMatrix m = arg.sparse_complex_matrix_value (); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
213 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
214 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
215 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
216 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
217 SparseComplexCHOL fact (m, info, natural); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
218 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
219 if (nargout == 3) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
220 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
221 if (vecout) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
222 retval(2) = fact.perm (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
223 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
224 retval(2) = fact.Q(); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
225 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
226 |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
227 if (nargout > 1 || info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
228 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
229 retval(1) = fact.P(); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
230 if (LLt) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
231 retval(0) = fact.L(); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
232 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
233 retval(0) = fact.R(); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
234 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
235 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
236 error ("chol: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
237 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
238 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
239 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
240 gripe_wrong_type_arg ("chol", arg); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
241 } |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7700
diff
changeset
|
242 else if (arg.is_single_type ()) |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
243 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
244 if (arg.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
245 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
246 FloatMatrix m = arg.float_matrix_value (); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7700
diff
changeset
|
247 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
248 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
249 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
250 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
251 FloatCHOL fact (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
252 if (nargout == 2 || info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
253 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
254 retval(1) = info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
255 if (LLt) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
256 retval(0) = get_chol_l (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
257 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
258 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
259 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
260 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
261 error ("chol: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
262 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
263 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
264 else if (arg.is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
265 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
266 FloatComplexMatrix m = arg.float_complex_matrix_value (); |
7789
82be108cc558
First attempt at single precision tyeps
David Bateman <dbateman@free.fr>
parents:
7700
diff
changeset
|
267 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
268 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
269 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
270 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
271 FloatComplexCHOL fact (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
272 if (nargout == 2 || info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
273 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
274 retval(1) = info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
275 if (LLt) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
276 retval(0) = get_chol_l (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
277 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
278 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
279 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
280 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
281 error ("chol: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
282 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
283 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
284 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
285 gripe_wrong_type_arg ("chol", arg); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
286 } |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
287 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
288 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
289 if (arg.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
290 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
291 Matrix m = arg.matrix_value (); |
2928 | 292 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
293 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
294 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
295 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
296 CHOL fact (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
297 if (nargout == 2 || info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
298 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
299 retval(1) = info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
300 if (LLt) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
301 retval(0) = get_chol_l (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
302 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
303 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
304 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
305 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
306 error ("chol: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
307 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
308 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
309 else if (arg.is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
310 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
311 ComplexMatrix m = arg.complex_matrix_value (); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
312 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
313 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
314 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
315 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
316 ComplexCHOL fact (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
317 if (nargout == 2 || info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
318 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
319 retval(1) = info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
320 if (LLt) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
321 retval(0) = get_chol_l (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
322 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
323 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
324 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
325 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
326 error ("chol: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
327 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
328 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
329 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
330 gripe_wrong_type_arg ("chol", arg); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
331 } |
2928 | 332 } |
333 | |
334 return retval; | |
335 } | |
336 | |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
337 /* |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
338 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
339 %!assert(chol ([2, 1; 1, 1]), [sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)], sqrt (eps)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
340 %!assert(chol (single([2, 1; 1, 1])), single([sqrt(2), 1/sqrt(2); 0, 1/sqrt(2)]), sqrt (eps('single'))); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
341 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
342 %!error chol ([1, 2; 3, 4]); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
343 %!error chol ([1, 2; 3, 4; 5, 6]); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
344 %!error <Invalid call to chol.*> chol (); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
345 %!error <unexpected second or third input.*> chol (1, 2); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
346 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
347 */ |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
348 |
5760 | 349 DEFUN_DLD (cholinv, args, , |
5340 | 350 "-*- texinfo -*-\n\ |
351 @deftypefn {Loadable Function} {} cholinv (@var{a})\n\ | |
5448 | 352 Use the Cholesky factorization to compute the inverse of the\n\ |
5340 | 353 symmetric positive definite matrix @var{a}.\n\ |
354 @seealso{chol, chol2inv}\n\ | |
355 @end deftypefn") | |
356 { | |
357 octave_value retval; | |
358 | |
359 int nargin = args.length (); | |
360 | |
361 if (nargin == 1) | |
362 { | |
363 octave_value arg = args(0); | |
364 | |
365 octave_idx_type nr = arg.rows (); | |
366 octave_idx_type nc = arg.columns (); | |
367 | |
368 if (nr == 0 || nc == 0) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
369 retval = Matrix (); |
5340 | 370 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
371 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
372 if (arg.is_sparse_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
373 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
374 if (arg.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
375 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
376 SparseMatrix m = arg.sparse_matrix_value (); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
377 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
378 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
379 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
380 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
381 SparseCHOL chol (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
382 if (info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
383 retval = chol.inverse (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
384 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
385 error ("cholinv: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
386 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
387 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
388 else if (arg.is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
389 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
390 SparseComplexMatrix m = arg.sparse_complex_matrix_value (); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
391 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
392 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
393 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
394 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
395 SparseComplexCHOL chol (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
396 if (info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
397 retval = chol.inverse (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
398 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
399 error ("cholinv: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
400 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
401 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
402 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
403 gripe_wrong_type_arg ("cholinv", arg); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
404 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
405 else if (arg.is_single_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
406 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
407 if (arg.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
408 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
409 FloatMatrix m = arg.float_matrix_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
410 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
411 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
412 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
413 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
414 FloatCHOL chol (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
415 if (info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
416 retval = chol.inverse (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
417 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
418 error ("cholinv: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
419 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
420 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
421 else if (arg.is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
422 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
423 FloatComplexMatrix m = arg.float_complex_matrix_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
424 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
425 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
426 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
427 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
428 FloatComplexCHOL chol (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
429 if (info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
430 retval = chol.inverse (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
431 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
432 error ("cholinv: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
433 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
434 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
435 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
436 gripe_wrong_type_arg ("chol", arg); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
437 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
438 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
439 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
440 if (arg.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
441 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
442 Matrix m = arg.matrix_value (); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
443 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
444 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
445 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
446 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
447 CHOL chol (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
448 if (info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
449 retval = chol.inverse (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
450 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
451 error ("cholinv: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
452 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
453 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
454 else if (arg.is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
455 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
456 ComplexMatrix m = arg.complex_matrix_value (); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
457 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
458 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
459 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
460 octave_idx_type info; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
461 ComplexCHOL chol (m, info); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
462 if (info == 0) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
463 retval = chol.inverse (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
464 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
465 error ("cholinv: matrix not positive definite"); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
466 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
467 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
468 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
469 gripe_wrong_type_arg ("chol", arg); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
470 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
471 } |
5340 | 472 } |
473 else | |
5823 | 474 print_usage (); |
5340 | 475 |
476 return retval; | |
477 } | |
478 | |
8402
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
479 /* |
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
480 |
8871 | 481 %!shared A, Ainv |
8402
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
482 %! A = [2,0.2;0.2,1]; |
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
483 %! Ainv = inv(A); |
8871 | 484 %!test |
8402
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
485 %! Ainv1 = cholinv(A); |
8871 | 486 %! assert (norm(Ainv-Ainv1),0,1e-10) |
487 %!testif HAVE_CHOLMOD | |
8402
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
488 %! Ainv2 = inv(sparse(A)); |
8871 | 489 %! assert (norm(Ainv-Ainv2),0,1e-10) |
8880
078ca05e4ef8
chol.cc: correct spelling of CHOLMOD in tests
John W. Eaton <jwe@octave.org>
parents:
8872
diff
changeset
|
490 %!testif HAVE_CHOLMOD |
8402
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
491 %! Ainv3 = cholinv(sparse(A)); |
8871 | 492 %! assert (norm(Ainv-Ainv3),0,1e-10) |
8880
078ca05e4ef8
chol.cc: correct spelling of CHOLMOD in tests
John W. Eaton <jwe@octave.org>
parents:
8872
diff
changeset
|
493 %!testif HAVE_CHOLMOD |
8402
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
494 %! Ainv4 = spcholinv(sparse(A)); |
8421
b8ce8738a4d1
chol.cc: Trivial fix for chol tests.
Ben Abbott <bpabbott@mac.com>
parents:
8402
diff
changeset
|
495 %! assert (norm(Ainv-Ainv4),0,1e-10) |
8402
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
496 |
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
497 */ |
2176f2b4599e
Fix sparse cholesky inversion
David Bateman <dbateman@free.fr>
parents:
8284
diff
changeset
|
498 |
5760 | 499 DEFUN_DLD (chol2inv, args, , |
5340 | 500 "-*- texinfo -*-\n\ |
5343 | 501 @deftypefn {Loadable Function} {} chol2inv (@var{u})\n\ |
5340 | 502 Invert a symmetric, positive definite square matrix from its Cholesky\n\ |
5343 | 503 decomposition, @var{u}. Note that @var{u} should be an upper-triangular\n\ |
504 matrix with positive diagonal elements. @code{chol2inv (@var{u})}\n\ | |
505 provides @code{inv (@var{u}'*@var{u})} but it is much faster than\n\ | |
506 using @code{inv}.\n\ | |
5340 | 507 @seealso{chol, cholinv}\n\ |
508 @end deftypefn") | |
509 { | |
510 octave_value retval; | |
511 | |
512 int nargin = args.length (); | |
513 | |
514 if (nargin == 1) | |
515 { | |
516 octave_value arg = args(0); | |
517 | |
518 octave_idx_type nr = arg.rows (); | |
519 octave_idx_type nc = arg.columns (); | |
520 | |
521 if (nr == 0 || nc == 0) | |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
522 retval = Matrix (); |
5340 | 523 else |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
524 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
525 if (arg.is_sparse_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
526 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
527 if (arg.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
528 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
529 SparseMatrix r = arg.sparse_matrix_value (); |
5340 | 530 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
531 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
532 retval = chol2inv (r); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
533 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
534 else if (arg.is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
535 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
536 SparseComplexMatrix r = arg.sparse_complex_matrix_value (); |
5340 | 537 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
538 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
539 retval = chol2inv (r); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
540 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
541 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
542 gripe_wrong_type_arg ("chol2inv", arg); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
543 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
544 else if (arg.is_single_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
545 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
546 if (arg.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
547 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
548 FloatMatrix r = arg.float_matrix_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
549 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
550 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
551 retval = chol2inv (r); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
552 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
553 else if (arg.is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
554 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
555 FloatComplexMatrix r = arg.float_complex_matrix_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
556 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
557 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
558 retval = chol2inv (r); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
559 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
560 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
561 gripe_wrong_type_arg ("chol2inv", arg); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
562 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
563 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
564 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
565 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
566 if (arg.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
567 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
568 Matrix r = arg.matrix_value (); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
569 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
570 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
571 retval = chol2inv (r); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
572 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
573 else if (arg.is_complex_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
574 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
575 ComplexMatrix r = arg.complex_matrix_value (); |
7515
f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
576 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
577 if (! error_state) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
578 retval = chol2inv (r); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
579 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
580 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
581 gripe_wrong_type_arg ("chol2inv", arg); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
582 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
583 } |
5340 | 584 } |
585 else | |
5823 | 586 print_usage (); |
5340 | 587 |
588 return retval; | |
589 } | |
590 | |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
591 DEFUN_DLD (cholupdate, args, nargout, |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
592 "-*- texinfo -*-\n\ |
7650 | 593 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholupdate (@var{R}, @var{u}, @var{op})\n\ |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
594 Update or downdate a Cholesky factorization. Given an upper triangular\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
595 matrix @var{R} and a column vector @var{u}, attempt to determine another\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
596 upper triangular matrix @var{R1} such that\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
597 @itemize @bullet\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
598 @item\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
599 @var{R1}'*@var{R1} = @var{R}'*@var{R} + @var{u}*@var{u}'\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
600 if @var{op} is \"+\"\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
601 @item\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
602 @var{R1}'*@var{R1} = @var{R}'*@var{R} - @var{u}*@var{u}'\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
603 if @var{op} is \"-\"\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
604 @end itemize\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
605 \n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
606 If @var{op} is \"-\", @var{info} is set to\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
607 @itemize\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
608 @item 0 if the downdate was successful,\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
609 @item 1 if @var{R}'*@var{R} - @var{u}*@var{u}' is not positive definite,\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
610 @item 2 if @var{R} is singular.\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
611 @end itemize\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
612 \n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
613 If @var{info} is not present, an error message is printed in cases 1 and 2.\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
614 @seealso{chol, qrupdate}\n\ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
615 @end deftypefn") |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
616 { |
7559
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
617 octave_idx_type nargin = args.length (); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
618 |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
619 octave_value_list retval; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
620 |
7559
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
621 if (nargin > 3 || nargin < 2) |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
622 { |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
623 print_usage (); |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
624 return retval; |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
625 } |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
626 |
7559
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
627 octave_value argr = args(0); |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
628 octave_value argu = args(1); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
629 |
9453
5fcf008c7c41
fix Cholesky updating with scalars
Jaroslav Hajek <highegg@gmail.com>
parents:
9065
diff
changeset
|
630 if (argr.is_numeric_type () && argu.is_numeric_type () |
7559
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
631 && (nargin < 3 || args(2).is_string ())) |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
632 { |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
633 octave_idx_type n = argr.rows (); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
634 |
7559
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
635 std::string op = (nargin < 3) ? "+" : args(2).string_value (); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
636 |
7559
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
637 bool down = op == "-"; |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
638 |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
639 if (down || op == "+") |
07522d7dcdf8
fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents:
7554
diff
changeset
|
640 if (argr.columns () == n && argu.rows () == n && argu.columns () == 1) |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
641 { |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
642 int err = 0; |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
643 if (argr.is_single_type () || argu.is_single_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
644 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
645 if (argr.is_real_type () && argu.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
646 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
647 // real case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
648 FloatMatrix R = argr.float_matrix_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
649 FloatColumnVector u = argu.float_column_vector_value (); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
650 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
651 FloatCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
652 fact.set (R); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
653 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
654 if (down) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
655 err = fact.downdate (u); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
656 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
657 fact.update (u); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
658 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
659 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
660 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
661 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
662 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
663 // complex case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
664 FloatComplexMatrix R = argr.float_complex_matrix_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
665 FloatComplexColumnVector u = argu.float_complex_column_vector_value (); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
666 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
667 FloatComplexCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
668 fact.set (R); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
669 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
670 if (down) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
671 err = fact.downdate (u); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
672 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
673 fact.update (u); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
674 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
675 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
676 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
677 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
678 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
679 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
680 if (argr.is_real_type () && argu.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
681 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
682 // real case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
683 Matrix R = argr.matrix_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
684 ColumnVector u = argu.column_vector_value (); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
685 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
686 CHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
687 fact.set (R); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
688 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
689 if (down) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
690 err = fact.downdate (u); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
691 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
692 fact.update (u); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
693 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
694 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
695 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
696 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
697 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
698 // complex case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
699 ComplexMatrix R = argr.complex_matrix_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
700 ComplexColumnVector u = argu.complex_column_vector_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
701 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
702 ComplexCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
703 fact.set (R); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
704 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
705 if (down) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
706 err = fact.downdate (u); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
707 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
708 fact.update (u); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
709 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
710 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
711 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
712 } |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
713 |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
714 if (nargout > 1) |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
715 retval(1) = err; |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
716 else if (err == 1) |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
717 error ("cholupdate: downdate violates positiveness"); |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
718 else if (err == 2) |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
719 error ("cholupdate: singular matrix"); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
720 } |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
721 else |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
722 error ("cholupdate: dimension mismatch"); |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
723 else |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
724 error ("cholupdate: op must be \"+\" or \"-\""); |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
725 } |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
726 else |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
727 print_usage (); |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
728 |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
729 return retval; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
730 } |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
731 |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
732 /* |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
733 %!shared A, u, Ac, uc |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
734 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
735 %! -0.131721 0.738529 0.019851 -0.140295 ; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
736 %! 0.124120 0.019851 0.354879 -0.059472 ; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
737 %! -0.061673 -0.140295 -0.059472 0.600939 ]; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
738 %! |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
739 %! u = [ 0.98950 ; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
740 %! 0.39844 ; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
741 %! 0.63484 ; |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
742 %! 0.13351 ]; |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
743 %! Ac = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
744 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
745 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
746 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ]; |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
747 %! |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
748 %! uc = [ 0.54267 + 0.91519i ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
749 %! 0.99647 + 0.43141i ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
750 %! 0.83760 + 0.68977i ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
751 %! 0.39160 + 0.90378i ]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
752 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
753 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
754 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
755 %!test |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
756 %! R = chol(A); |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
757 %! |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
758 %! R1 = cholupdate(R,u); |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
759 %! |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
760 %! assert(norm(triu(R1)-R1,Inf) == 0) |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
761 %! assert(norm(R1'*R1 - R'*R - u*u',Inf) < 1e1*eps) |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
762 %! |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
763 %! R1 = cholupdate(R1,u,"-"); |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
764 %! |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
765 %! assert(norm(triu(R1)-R1,Inf) == 0) |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
766 %! assert(norm(R1 - R,Inf) < 1e1*eps) |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
767 %! |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
768 %!test |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
769 %! R = chol(Ac); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
770 %! |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
771 %! R1 = cholupdate(R,uc); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
772 %! |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
773 %! assert(norm(triu(R1)-R1,Inf) == 0) |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
774 %! assert(norm(R1'*R1 - R'*R - uc*uc',Inf) < 1e1*eps) |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
775 %! |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
776 %! R1 = cholupdate(R1,uc,"-"); |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
777 %! |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
778 %! assert(norm(triu(R1)-R1,Inf) == 0) |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
779 %! assert(norm(R1 - R,Inf) < 1e1*eps) |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
780 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
781 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
782 %! R = chol(single(A)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
783 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
784 %! R1 = cholupdate(R,single(u)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
785 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
786 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
787 %! assert(norm(R1'*R1 - R'*R - single(u*u'),Inf) < 1e1*eps('single')) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
788 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
789 %! R1 = cholupdate(R1,single(u),"-"); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
790 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
791 %! assert(norm(triu(R1)-R1,Inf) == 0) |
8066
366821c0c01c
Adjust tolerance to single precision test.
Ben Abbott <bpabbott@mac.com>
parents:
7924
diff
changeset
|
792 %! assert(norm(R1 - R,Inf) < 2e1*eps('single')) |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
793 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
794 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
795 %! R = chol(single(Ac)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
796 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
797 %! R1 = cholupdate(R,single(uc)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
798 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
799 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
800 %! assert(norm(R1'*R1 - R'*R - single(uc*uc'),Inf) < 1e1*eps('single')) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
801 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
802 %! R1 = cholupdate(R1,single(uc),"-"); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
803 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
804 %! assert(norm(triu(R1)-R1,Inf) == 0) |
8066
366821c0c01c
Adjust tolerance to single precision test.
Ben Abbott <bpabbott@mac.com>
parents:
7924
diff
changeset
|
805 %! assert(norm(R1 - R,Inf) < 2e1*eps('single')) |
7554
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
806 */ |
40574114c514
implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents:
7515
diff
changeset
|
807 |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
808 DEFUN_DLD (cholinsert, args, nargout, |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
809 "-*- texinfo -*-\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
810 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholinsert (@var{R}, @var{j}, @var{u})\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
811 Given a Cholesky@tie{}factorization of a real symmetric or complex hermitian\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
812 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper triangular,\n\ |
8284
4ceffd54031a
fix docs for cholinsert, choldelete, cholshift
Jaroslav Hajek <highegg@gmail.com>
parents:
8066
diff
changeset
|
813 return the Cholesky@tie{}factorization of\n\ |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
814 @var{A1}, where @w{A1(p,p) = A}, @w{A1(:,j) = A1(j,:)' = u} and\n\ |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8880
diff
changeset
|
815 @w{p = [1:j-1,j+1:n+1]}. @w{u(j)} should be positive.\n\ |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
816 On return, @var{info} is set to\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
817 @itemize\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
818 @item 0 if the insertion was successful,\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
819 @item 1 if @var{A1} is not positive definite,\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
820 @item 2 if @var{R} is singular.\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
821 @end itemize\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
822 \n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
823 If @var{info} is not present, an error message is printed in cases 1 and 2.\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
824 @seealso{chol, cholupdate, choldelete}\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
825 @end deftypefn") |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
826 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
827 octave_idx_type nargin = args.length (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
828 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
829 octave_value_list retval; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
830 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
831 if (nargin != 3) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
832 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
833 print_usage (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
834 return retval; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
835 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
836 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
837 octave_value argr = args(0); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
838 octave_value argj = args(1); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
839 octave_value argu = args(2); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
840 |
9453
5fcf008c7c41
fix Cholesky updating with scalars
Jaroslav Hajek <highegg@gmail.com>
parents:
9065
diff
changeset
|
841 if (argr.is_numeric_type () && argu.is_numeric_type () |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
842 && argj.is_real_scalar ()) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
843 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
844 octave_idx_type n = argr.rows (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
845 octave_idx_type j = argj.scalar_value (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
846 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
847 if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
848 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
849 if (j > 0 && j <= n+1) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
850 { |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
851 int err = 0; |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
852 if (argr.is_single_type () || argu.is_single_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
853 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
854 if (argr.is_real_type () && argu.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
855 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
856 // real case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
857 FloatMatrix R = argr.float_matrix_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
858 FloatColumnVector u = argu.float_column_vector_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
859 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
860 FloatCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
861 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
862 err = fact.insert_sym (u, j-1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
863 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
864 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
865 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
866 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
867 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
868 // complex case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
869 FloatComplexMatrix R = argr.float_complex_matrix_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
870 FloatComplexColumnVector u = argu.float_complex_column_vector_value (); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
871 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
872 FloatComplexCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
873 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
874 err = fact.insert_sym (u, j-1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
875 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
876 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
877 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
878 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
879 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
880 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
881 if (argr.is_real_type () && argu.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
882 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
883 // real case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
884 Matrix R = argr.matrix_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
885 ColumnVector u = argu.column_vector_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
886 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
887 CHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
888 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
889 err = fact.insert_sym (u, j-1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
890 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
891 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
892 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
893 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
894 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
895 // complex case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
896 ComplexMatrix R = argr.complex_matrix_value (); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
897 ComplexColumnVector u = argu.complex_column_vector_value (); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
898 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
899 ComplexCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
900 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
901 err = fact.insert_sym (u, j-1); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
902 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
903 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
904 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
905 } |
8547
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
906 |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
907 if (nargout > 1) |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
908 retval(1) = err; |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
909 else if (err == 1) |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
910 error ("cholinsert: insertion violates positiveness"); |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
911 else if (err == 2) |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
912 error ("cholinsert: singular matrix"); |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
913 else if (err == 3) |
d66c9b6e506a
imported patch qrupdate.diff
Jaroslav Hajek <highegg@gmail.com>
parents:
8517
diff
changeset
|
914 error ("cholinsert: diagonal element must be real"); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
915 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
916 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
917 error ("cholinsert: index out of range"); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
918 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
919 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
920 error ("cholinsert: dimension mismatch"); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
921 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
922 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
923 print_usage (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
924 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
925 return retval; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
926 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
927 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
928 /* |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
929 %!test |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
930 %! u2 = [ 0.35080 ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
931 %! 0.63930 ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
932 %! 3.31057 ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
933 %! -0.13825 ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
934 %! 0.45266 ]; |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
935 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
936 %! R = chol(A); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
937 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
938 %! j = 3; p = [1:j-1, j+1:5]; |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
939 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1; |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
940 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
941 %! assert(norm(triu(R1)-R1,Inf) == 0) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
942 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
943 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
944 %!test |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
945 %! u2 = [ 0.35080 + 0.04298i; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
946 %! 0.63930 + 0.23778i; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
947 %! 3.31057 + 0.00000i; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
948 %! -0.13825 + 0.19879i; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
949 %! 0.45266 + 0.50020i]; |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
950 %! |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
951 %! R = chol(Ac); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
952 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
953 %! j = 3; p = [1:j-1, j+1:5]; |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
954 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1; |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
955 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
956 %! assert(norm(triu(R1)-R1,Inf) == 0) |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
957 %! assert(norm(A1(p,p) - Ac,Inf) < 1e1*eps) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
958 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
959 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
960 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
961 %! u2 = single ([ 0.35080 ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
962 %! 0.63930 ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
963 %! 3.31057 ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
964 %! -0.13825 ; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
965 %! 0.45266 ]); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
966 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
967 %! R = chol(single(A)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
968 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
969 %! j = 3; p = [1:j-1, j+1:5]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
970 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
971 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
972 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
973 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps('single')) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
974 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
975 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
976 %! u2 = single ([ 0.35080 + 0.04298i; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
977 %! 0.63930 + 0.23778i; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
978 %! 3.31057 + 0.00000i; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
979 %! -0.13825 + 0.19879i; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
980 %! 0.45266 + 0.50020i]); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
981 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
982 %! R = chol(single(Ac)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
983 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
984 %! j = 3; p = [1:j-1, j+1:5]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
985 %! R1 = cholinsert(R,j,u2); A1 = R1'*R1; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
986 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
987 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
988 %! assert(norm(A1(p,p) - single(Ac),Inf) < 1e1*eps('single')) |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
989 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
990 */ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
991 |
7924 | 992 DEFUN_DLD (choldelete, args, , |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
993 "-*- texinfo -*-\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
994 @deftypefn {Loadable Function} {@var{R1} =} choldelete (@var{R}, @var{j})\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
995 Given a Cholesky@tie{}factorization of a real symmetric or complex hermitian\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
996 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper triangular,\n\ |
8284
4ceffd54031a
fix docs for cholinsert, choldelete, cholshift
Jaroslav Hajek <highegg@gmail.com>
parents:
8066
diff
changeset
|
997 return the Cholesky@tie{}factorization of @w{A(p,p)}, where @w{p = [1:j-1,j+1:n+1]}.\n\ |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
998 @seealso{chol, cholupdate, cholinsert}\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
999 @end deftypefn") |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1000 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1001 octave_idx_type nargin = args.length (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1002 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1003 octave_value_list retval; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1004 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1005 if (nargin != 2) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1006 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1007 print_usage (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1008 return retval; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1009 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1010 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1011 octave_value argr = args(0); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1012 octave_value argj = args(1); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1013 |
9453
5fcf008c7c41
fix Cholesky updating with scalars
Jaroslav Hajek <highegg@gmail.com>
parents:
9065
diff
changeset
|
1014 if (argr.is_numeric_type () && argj.is_real_scalar ()) |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1015 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1016 octave_idx_type n = argr.rows (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1017 octave_idx_type j = argj.scalar_value (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1018 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1019 if (argr.columns () == n) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1020 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1021 if (j > 0 && j <= n) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1022 { |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1023 if (argr.is_single_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1024 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1025 if (argr.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1026 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1027 // real case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1028 FloatMatrix R = argr.float_matrix_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1029 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1030 FloatCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1031 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1032 fact.delete_sym (j-1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1033 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1034 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1035 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1036 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1037 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1038 // complex case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1039 FloatComplexMatrix R = argr.float_complex_matrix_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1040 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1041 FloatComplexCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1042 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1043 fact.delete_sym (j-1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1044 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1045 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1046 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1047 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1048 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1049 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1050 if (argr.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1051 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1052 // real case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1053 Matrix R = argr.matrix_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1054 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1055 CHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1056 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1057 fact.delete_sym (j-1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1058 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1059 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1060 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1061 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1062 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1063 // complex case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1064 ComplexMatrix R = argr.complex_matrix_value (); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1065 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1066 ComplexCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1067 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1068 fact.delete_sym (j-1); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1069 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1070 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1071 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1072 } |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1073 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1074 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1075 error ("choldelete: index out of range"); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1076 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1077 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1078 error ("choldelete: dimension mismatch"); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1079 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1080 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1081 print_usage (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1082 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1083 return retval; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1084 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1085 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1086 /* |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1087 %!test |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1088 %! R = chol(A); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1089 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1090 %! j = 3; p = [1:j-1,j+1:4]; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1091 %! R1 = choldelete(R,j); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1092 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1093 %! assert(norm(triu(R1)-R1,Inf) == 0) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1094 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1095 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1096 %!test |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1097 %! R = chol(Ac); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1098 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1099 %! j = 3; p = [1:j-1,j+1:4]; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1100 %! R1 = choldelete(R,j); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1101 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1102 %! assert(norm(triu(R1)-R1,Inf) == 0) |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1103 %! assert(norm(R1'*R1 - Ac(p,p),Inf) < 1e1*eps) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1104 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1105 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1106 %! R = chol(single(A)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1107 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1108 %! j = 3; p = [1:j-1,j+1:4]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1109 %! R1 = choldelete(R,j); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1110 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1111 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1112 %! assert(norm(R1'*R1 - single(A(p,p)),Inf) < 1e1*eps('single')) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1113 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1114 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1115 %! R = chol(single(Ac)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1116 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1117 %! j = 3; p = [1:j-1,j+1:4]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1118 %! R1 = choldelete(R,j); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1119 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1120 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1121 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single')) |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1122 */ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1123 |
7924 | 1124 DEFUN_DLD (cholshift, args, , |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1125 "-*- texinfo -*-\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1126 @deftypefn {Loadable Function} {@var{R1} =} cholshift (@var{R}, @var{i}, @var{j})\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1127 Given a Cholesky@tie{}factorization of a real symmetric or complex hermitian\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1128 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper triangular,\n\ |
8284
4ceffd54031a
fix docs for cholinsert, choldelete, cholshift
Jaroslav Hajek <highegg@gmail.com>
parents:
8066
diff
changeset
|
1129 return the Cholesky@tie{}factorization of\n\ |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1130 @w{@var{A}(p,p)}, where @w{p} is the permutation @*\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1131 @code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @*\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1132 or @*\n\ |
9064
7c02ec148a3c
Check grammar on all .cc files
Rik <rdrider0-list@yahoo.com>
parents:
8880
diff
changeset
|
1133 @code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}. @*\n\ |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1134 \n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1135 @seealso{chol, cholinsert, choldelete}\n\ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1136 @end deftypefn") |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1137 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1138 octave_idx_type nargin = args.length (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1139 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1140 octave_value_list retval; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1141 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1142 if (nargin != 3) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1143 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1144 print_usage (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1145 return retval; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1146 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1147 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1148 octave_value argr = args(0); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1149 octave_value argi = args(1); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1150 octave_value argj = args(2); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1151 |
9453
5fcf008c7c41
fix Cholesky updating with scalars
Jaroslav Hajek <highegg@gmail.com>
parents:
9065
diff
changeset
|
1152 if (argr.is_numeric_type () && argi.is_real_scalar () && argj.is_real_scalar ()) |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1153 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1154 octave_idx_type n = argr.rows (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1155 octave_idx_type i = argi.scalar_value (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1156 octave_idx_type j = argj.scalar_value (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1157 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1158 if (argr.columns () == n) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1159 { |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1160 if (j > 0 && j <= n+1 && i > 0 && i <= n+1) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1161 { |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1162 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1163 if (argr.is_single_type () && argi.is_single_type () && |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1164 argj.is_single_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1165 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1166 if (argr.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1167 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1168 // real case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1169 FloatMatrix R = argr.float_matrix_value (); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1170 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1171 FloatCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1172 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1173 fact.shift_sym (i-1, j-1); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1174 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1175 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1176 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1177 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1178 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1179 // complex case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1180 FloatComplexMatrix R = argr.float_complex_matrix_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1181 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1182 FloatComplexCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1183 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1184 fact.shift_sym (i-1, j-1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1185 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1186 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1187 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1188 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1189 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1190 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1191 if (argr.is_real_type ()) |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1192 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1193 // real case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1194 Matrix R = argr.matrix_value (); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1195 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1196 CHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1197 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1198 fact.shift_sym (i-1, j-1); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1199 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1200 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1201 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1202 else |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1203 { |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1204 // complex case |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1205 ComplexMatrix R = argr.complex_matrix_value (); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1206 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1207 ComplexCHOL fact; |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1208 fact.set (R); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1209 fact.shift_sym (i-1, j-1); |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1210 |
10154
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1211 retval(0) = get_chol_r (fact); |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1212 } |
40dfc0c99116
DLD-FUNCTIONS/*.cc: untabify
John W. Eaton <jwe@octave.org>
parents:
9862
diff
changeset
|
1213 } |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1214 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1215 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1216 error ("cholshift: index out of range"); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1217 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1218 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1219 error ("cholshift: dimension mismatch"); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1220 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1221 else |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1222 print_usage (); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1223 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1224 return retval; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1225 } |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1226 |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1227 /* |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1228 %!test |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1229 %! R = chol(A); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1230 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1231 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4]; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1232 %! R1 = cholshift(R,i,j); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1233 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1234 %! assert(norm(triu(R1)-R1,Inf) == 0) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1235 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1236 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1237 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4]; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1238 %! R1 = cholshift(R,i,j); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1239 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1240 %! assert(norm(triu(R1)-R1,Inf) == 0) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1241 %! assert(norm(R1'*R1 - A(p,p),Inf) < 1e1*eps) |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1242 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1243 %!test |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1244 %! R = chol(Ac); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1245 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1246 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1247 %! R1 = cholshift(R,i,j); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1248 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1249 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1250 %! assert(norm(R1'*R1 - Ac(p,p),Inf) < 1e1*eps) |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1251 %! |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1252 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1253 %! R1 = cholshift(R,i,j); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1254 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1255 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1256 %! assert(norm(R1'*R1 - Ac(p,p),Inf) < 1e1*eps) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1257 |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1258 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1259 %! R = chol(single(A)); |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1260 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1261 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4]; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1262 %! R1 = cholshift(R,i,j); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1263 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1264 %! assert(norm(triu(R1)-R1,Inf) == 0) |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1265 %! assert(norm(R1'*R1 - single(A(p,p)),Inf) < 1e1*eps('single')) |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1266 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1267 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4]; |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1268 %! R1 = cholshift(R,i,j); |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1269 %! |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1270 %! assert(norm(triu(R1)-R1,Inf) == 0) |
7814
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1271 %! assert(norm(R1'*R1 - single(A(p,p)),Inf) < 1e1*eps('single')) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1272 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1273 %!test |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1274 %! R = chol(single(Ac)); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1275 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1276 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1277 %! R1 = cholshift(R,i,j); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1278 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1279 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1280 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single')) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1281 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1282 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4]; |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1283 %! R1 = cholshift(R,i,j); |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1284 %! |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1285 %! assert(norm(triu(R1)-R1,Inf) == 0) |
87865ed7405f
Second set of single precision test code and fix of resulting bugs
David Bateman <dbateman@free.fr>
parents:
7789
diff
changeset
|
1286 %! assert(norm(R1'*R1 - single(Ac(p,p)),Inf) < 1e1*eps('single')) |
7700
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1287 */ |
efccca5f2ad7
more QR & Cholesky updating functions
Jaroslav Hajek <highegg@gmail.com>
parents:
7650
diff
changeset
|
1288 |
2928 | 1289 /* |
1290 ;;; Local Variables: *** | |
1291 ;;; mode: C++ *** | |
1292 ;;; End: *** | |
1293 */ | |
1294 |