Mercurial > octave-nkf

annotate src/DLD-FUNCTIONS/chol.cc @ 7650:eb7bdde776f2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Texinfo fixes
author John W. Eaton <jwe@octave.org>
date Wed, 26 Mar 2008 16:27:14 -0400
parents 07522d7dcdf8
children efccca5f2ad7
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
1 /*
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
2
7017
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
3 Copyright (C) 1996, 1997, 1999, 2000, 2002, 2005, 2006, 2007
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
4 John W. Eaton
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
5
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
6 This file is part of Octave.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
7
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
8 Octave is free software; you can redistribute it and/or modify it
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
9 under the terms of the GNU General Public License as published by the
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5823
diff changeset
10 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5823
diff changeset
11 option) any later version.
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
12
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
13 Octave is distributed in the hope that it will be useful, but WITHOUT
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
16 for more details.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
17
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
18 You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5823
diff changeset
19 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 5823
diff changeset
20 <http://www.gnu.org/licenses/>.
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
21
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
22 */
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
23
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
24 // The cholupdate function was written by Jaroslav Hajek
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
25 // <highegg@gmail.com>, Copyright (C) 2008 VZLU Prague, a.s., Czech
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
26 // Republic.
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
27
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
28 #ifdef HAVE_CONFIG_H
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
29 #include <config.h>
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
30 #endif
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
31
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
32 #include "CmplxCHOL.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
33 #include "dbleCHOL.h"
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
34 #include "SparseCmplxCHOL.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
35 #include "SparsedbleCHOL.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
36 #include "oct-spparms.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
37 #include "sparse-util.h"
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
38
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
39 #include "ov-re-sparse.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
40 #include "ov-cx-sparse.h"
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
41 #include "defun-dld.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
42 #include "error.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
43 #include "gripes.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
44 #include "oct-obj.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
45 #include "utils.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
46
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
47 DEFUN_DLD (chol, args, nargout,
3548
ab7fa5a8f23f [project @ 2000-02-03 01:17:15 by jwe]
jwe
parents: 3372
diff changeset
48 "-*- texinfo -*-\n\
7650
eb7bdde776f2 Texinfo fixes
John W. Eaton <jwe@octave.org>
parents: 7559
diff changeset
49 @deftypefn {Loadable Function} {@var{r} =} chol (@var{a})\n\
eb7bdde776f2 Texinfo fixes
John W. Eaton <jwe@octave.org>
parents: 7559
diff changeset
50 @deftypefnx {Loadable Function} {[@var{r}, @var{p}] =} chol (@var{a})\n\
eb7bdde776f2 Texinfo fixes
John W. Eaton <jwe@octave.org>
parents: 7559
diff changeset
51 @deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}] =} chol (@var{s})\n\
eb7bdde776f2 Texinfo fixes
John W. Eaton <jwe@octave.org>
parents: 7559
diff changeset
52 @deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}] =} chol (@var{s}, 'vector')\n\
eb7bdde776f2 Texinfo fixes
John W. Eaton <jwe@octave.org>
parents: 7559
diff changeset
53 @deftypefnx {Loadable Function} {[@var{l}, @dots{}] =} chol (@dots{}, 'lower')\n\
3372
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
54 @cindex Cholesky factorization\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
55 Compute the Cholesky factor, @var{r}, of the symmetric positive definite\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
56 matrix @var{a}, where\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
57 @iftex\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
58 @tex\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
59 $ R^T R = A $.\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
60 @end tex\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
61 @end iftex\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
62 @ifinfo\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
63 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
64 @example\n\
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
65 @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
66 @end example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
67 @end ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
68 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
69 Called with one output argument @code{chol} fails if @var{a} or @var{s} is\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
70 not positive definite. With two or more output arguments @var{p} flags\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
71 whether the matrix was positive definite and @code{chol} does not fail. A\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
72 zero value indicated that the matrix was positive definite and the @var{r}\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
73 gives the factorization, annd @var{p} will have a positive value otherwise.\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
74 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
75 If called with 3 outputs then a sparsity preserving row/column permutation\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
76 is applied to @var{a} prior to the factorization. That is @var{r}\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
77 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
78 @iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
79 @tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
80 $ 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
81 @end tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
82 @end iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
83 @ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
84 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
85 @example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
86 @var{r}' * @var{r} = @var{q}' * @var{a} * @var{q}.\n\
3372
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
87 @end example\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
88 @end ifinfo\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 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
91 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
92 such that\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
93 @iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
94 @tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
95 $ 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
96 @end tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
97 @end iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
98 @ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
99 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
100 @example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
101 @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
102 @end example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
103 @end ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
104 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
105 Called with either a sparse or full matrix and uing the 'lower' flag,\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
106 @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
107 @iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
108 @tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
109 $ L L^T = A $.\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
110 @end tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
111 @end iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
112 @ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
113 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
114 @example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
115 @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
116 @end example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
117 @end ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
118 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
119 In general the lower trinagular factorization is significantly faster for\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
120 sparse matrices.\n\
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
121 @seealso{cholinv, chol2inv}\n\
3372
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe]
jwe
parents: 3244
diff changeset
122 @end deftypefn")
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
123 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
124 octave_value_list retval;
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
125 int nargin = args.length ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
126 bool LLt = false;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
127 bool vecout = false;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
128
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
129 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
130 || (! args(0).is_sparse_type () && nargout > 2))
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
131 {
5823
080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe]
jwe
parents: 5760
diff changeset
132 print_usage ();
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
133 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
134 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
135
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
136 int n = 1;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
137 while (n < nargin && ! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
138 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
139 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
140
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
141 if (! error_state )
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
142 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
143 if (tmp.compare ("vector") == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
144 vecout = true;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
145 else if (tmp.compare ("lower") == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
146 LLt = true;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
147 else if (tmp.compare ("upper") == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
148 LLt = false;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
149 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
150 error ("chol: unexpected second or third input");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
151 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
152 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
153 error ("chol: expecting trailing string arguments");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
154 }
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
155
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
156 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
157 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
158 octave_value arg = args(0);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
159
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
160 octave_idx_type nr = arg.rows ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
161 octave_idx_type nc = arg.columns ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
162 bool natural = (nargout != 3);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
163
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
164 int arg_is_empty = empty_arg ("chol", nr, nc);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
165
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
166 if (arg_is_empty < 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
167 return retval;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
168 if (arg_is_empty > 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
169 return octave_value (Matrix ());
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
170
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
171 if (arg.is_sparse_type ())
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 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
174 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
175 SparseMatrix m = arg.sparse_matrix_value ();
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
176
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
177 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
178 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
179 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
180 SparseCHOL fact (m, info, natural);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
181 if (nargout == 3)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
182 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
183 retval(2) = fact.perm ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
184 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
185 retval(2) = fact.Q();
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
186
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
187 if (nargout > 1 || info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
188 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
189 retval(1) = fact.P();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
190 if (LLt)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
191 retval(0) = fact.L();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
192 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
193 retval(0) = fact.R();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
194 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
195 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
196 error ("chol: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
197 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
198 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
199 else if (arg.is_complex_type ())
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3014
diff changeset
200 {
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
201 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
202
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
203 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
204 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
205 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
206 SparseComplexCHOL fact (m, info, natural);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
207
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
208 if (nargout == 3)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
209 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
210 retval(2) = fact.perm ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
211 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
212 retval(2) = fact.Q();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
213
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
214 if (nargout > 1 || info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
215 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
216 retval(1) = fact.P();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
217 if (LLt)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
218 retval(0) = fact.L();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
219 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
220 retval(0) = fact.R();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
221 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
222 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
223 error ("chol: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
224 }
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3014
diff changeset
225 }
3244
7f8827f4040a [project @ 1999-05-29 14:01:55 by jwe]
jwe
parents: 3243
diff changeset
226 else
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
227 gripe_wrong_type_arg ("chol", arg);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
228 }
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
229 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
230 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
231 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
232 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
233 Matrix m = arg.matrix_value ();
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
234
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
235 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
236 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
237 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
238 CHOL fact (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
239 if (nargout == 2 || info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
240 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
241 retval(1) = static_cast<double> (info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
242 if (LLt)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
243 retval(0) = fact.chol_matrix ().transpose ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
244 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
245 retval(0) = fact.chol_matrix ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
246 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
247 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
248 error ("chol: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
249 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
250 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
251 else if (arg.is_complex_type ())
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3014
diff changeset
252 {
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
253 ComplexMatrix m = arg.complex_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
254
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
255 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
256 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
257 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
258 ComplexCHOL fact (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
259 if (nargout == 2 || info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
260 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
261 retval(1) = static_cast<double> (info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
262 if (LLt)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
263 retval(0) = fact.chol_matrix ().hermitian ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
264 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
265 retval(0) = fact.chol_matrix ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
266 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
267 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
268 error ("chol: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
269 }
3243
dd00769643ae [project @ 1999-05-28 04:19:00 by jwe]
jwe
parents: 3014
diff changeset
270 }
3244
7f8827f4040a [project @ 1999-05-29 14:01:55 by jwe]
jwe
parents: 3243
diff changeset
271 else
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
272 gripe_wrong_type_arg ("chol", arg);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
273 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
274 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
275
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
276 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
277 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
278
5760
8d7162924bd3 [project @ 2006-04-14 04:01:37 by jwe]
jwe
parents: 5448
diff changeset
279 DEFUN_DLD (cholinv, args, ,
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
280 "-*- texinfo -*-\n\
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
281 @deftypefn {Loadable Function} {} cholinv (@var{a})\n\
5448
ebe5d7d15522 [project @ 2005-09-14 18:55:04 by jwe]
jwe
parents: 5343
diff changeset
282 Use the Cholesky factorization to compute the inverse of the\n\
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
283 symmetric positive definite matrix @var{a}.\n\
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
284 @seealso{chol, chol2inv}\n\
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
285 @end deftypefn")
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
286 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
287 octave_value retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
288
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
289 int nargin = args.length ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
290
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
291 if (nargin == 1)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
292 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
293 octave_value arg = args(0);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
294
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
295 octave_idx_type nr = arg.rows ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
296 octave_idx_type nc = arg.columns ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
297
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
298 if (nr == 0 || nc == 0)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
299 retval = Matrix ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
300 else
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
301 {
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
302 if (arg.is_sparse_type ())
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
303 {
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
304 if (arg.is_real_type ())
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
305 {
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
306 SparseMatrix m = arg.sparse_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
307
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
308 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
309 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
310 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
311 SparseCHOL chol (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
312 if (info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
313 retval = chol.inverse ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
314 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
315 error ("cholinv: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
316 }
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
317 }
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
318 else if (arg.is_complex_type ())
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
319 {
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
320 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
321
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
322 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
323 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
324 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
325 SparseComplexCHOL chol (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
326 if (info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
327 retval = chol.inverse ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
328 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
329 error ("cholinv: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
330 }
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
331 }
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
332 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
333 gripe_wrong_type_arg ("cholinv", arg);
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
334 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
335 else
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
336 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
337 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
338 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
339 Matrix m = arg.matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
340
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
341 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
342 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
343 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
344 CHOL chol (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
345 if (info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
346 retval = chol.inverse ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
347 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
348 error ("cholinv: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
349 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
350 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
351 else if (arg.is_complex_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
352 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
353 ComplexMatrix m = arg.complex_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
354
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
355 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
356 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
357 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
358 ComplexCHOL chol (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
359 if (info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
360 retval = chol.inverse ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
361 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
362 error ("cholinv: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
363 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
364 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
365 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
366 gripe_wrong_type_arg ("chol", arg);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
367 }
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
368 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
369 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
370 else
5823
080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe]
jwe
parents: 5760
diff changeset
371 print_usage ();
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
372
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
373 return retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
374 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
375
5760
8d7162924bd3 [project @ 2006-04-14 04:01:37 by jwe]
jwe
parents: 5448
diff changeset
376 DEFUN_DLD (chol2inv, args, ,
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
377 "-*- texinfo -*-\n\
5343
be0bd7a760c9 [project @ 2005-05-06 18:04:16 by jwe]
jwe
parents: 5340
diff changeset
378 @deftypefn {Loadable Function} {} chol2inv (@var{u})\n\
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
379 Invert a symmetric, positive definite square matrix from its Cholesky\n\
5343
be0bd7a760c9 [project @ 2005-05-06 18:04:16 by jwe]
jwe
parents: 5340
diff changeset
380 decomposition, @var{u}. Note that @var{u} should be an upper-triangular\n\
be0bd7a760c9 [project @ 2005-05-06 18:04:16 by jwe]
jwe
parents: 5340
diff changeset
381 matrix with positive diagonal elements. @code{chol2inv (@var{u})}\n\
be0bd7a760c9 [project @ 2005-05-06 18:04:16 by jwe]
jwe
parents: 5340
diff changeset
382 provides @code{inv (@var{u}'*@var{u})} but it is much faster than\n\
be0bd7a760c9 [project @ 2005-05-06 18:04:16 by jwe]
jwe
parents: 5340
diff changeset
383 using @code{inv}.\n\
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
384 @seealso{chol, cholinv}\n\
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
385 @end deftypefn")
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
386 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
387 octave_value retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
388
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
389 int nargin = args.length ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
390
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
391 if (nargin == 1)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
392 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
393 octave_value arg = args(0);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
394
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
395 octave_idx_type nr = arg.rows ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
396 octave_idx_type nc = arg.columns ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
397
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
398 if (nr == 0 || nc == 0)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
399 retval = Matrix ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
400 else
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
401 {
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
402 if (arg.is_sparse_type ())
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
403 {
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
404 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
405 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
406 SparseMatrix r = arg.sparse_matrix_value ();
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
407
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
408 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
409 retval = chol2inv (r);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
410 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
411 else if (arg.is_complex_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
412 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
413 SparseComplexMatrix r = arg.sparse_complex_matrix_value ();
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
414
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
415 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
416 retval = chol2inv (r);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
417 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
418 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
419 gripe_wrong_type_arg ("chol2inv", arg);
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
420 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
421 else
7515
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
422 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
423 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
424 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
425 Matrix r = arg.matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
426
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
427 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
428 retval = chol2inv (r);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
429 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
430 else if (arg.is_complex_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
431 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
432 ComplexMatrix r = arg.complex_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
433
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
434 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
435 retval = chol2inv (r);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
436 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
437 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
438 gripe_wrong_type_arg ("chol2inv", arg);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions
David Bateman <dbateman@free.fr>
parents: 7017
diff changeset
439 }
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
440 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
441 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
442 else
5823
080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe]
jwe
parents: 5760
diff changeset
443 print_usage ();
5340
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
444
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
445 return retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
446 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe]
jwe
parents: 5307
diff changeset
447
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
448 DEFUN_DLD (cholupdate, args, nargout,
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
449 "-*- texinfo -*-\n\
7650
eb7bdde776f2 Texinfo fixes
John W. Eaton <jwe@octave.org>
parents: 7559
diff changeset
450 @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
451 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
452 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
453 upper triangular matrix @var{R1} such that\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
454 @itemize @bullet\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
455 @item\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
456 @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
457 if @var{op} is \"+\"\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
458 @item\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
459 @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
460 if @var{op} is \"-\"\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
461 @end itemize\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
462 \n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
463 If @var{op} is \"-\", @var{info} is set to\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
464 @itemize\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
465 @item 0 if the downdate was successful,\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
466 @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
467 @item 2 if @var{R} is singular.\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
468 @end itemize\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
469 \n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
470 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
471 @seealso{chol, qrupdate}\n\
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
472 @end deftypefn")
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
473 {
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
474 octave_idx_type nargin = args.length ();
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
475
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
476 octave_value_list retval;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
477
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
478 if (nargin > 3 || nargin < 2)
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
479 {
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
480 print_usage ();
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
481 return retval;
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
482 }
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
483
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
484 octave_value argr = args(0);
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
485 octave_value argu = args(1);
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
486
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
487 if (argr.is_matrix_type () && argu.is_matrix_type ()
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
488 && (nargin < 3 || args(2).is_string ()))
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
489 {
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
490 octave_idx_type n = argr.rows ();
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
491
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
492 std::string op = (nargin < 3) ? "+" : args(2).string_value ();
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
493
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
494 bool down = op == "-";
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
495
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
496 if (down || op == "+")
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
497 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
498 {
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
499 if (argr.is_real_matrix () && argu.is_real_matrix ())
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
500 {
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
501 // real case
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
502 Matrix R = argr.matrix_value ();
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
503 Matrix u = argu.matrix_value ();
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
504
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
505 CHOL fact;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
506 fact.set (R);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
507 int err = 0;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
508
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
509 if (down)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
510 err = fact.downdate (u);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
511 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
512 fact.update (u);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
513
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
514 if (nargout > 1)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
515 retval(1) = err;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
516 else if (err)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
517 error ("cholupdate: downdate violates positiveness");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
518
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
519 retval(0) = fact.chol_matrix ();
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
520 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
521 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
522 {
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
523 // complex case
7559
07522d7dcdf8 fixes to QR and Cholesky updating code
Jaroslav Hajek <highegg@gmail.com>
parents: 7554
diff changeset
524 ComplexMatrix R = argr.complex_matrix_value ();
7554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
525 ComplexMatrix u = argu.complex_matrix_value ();
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
526
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
527 ComplexCHOL fact;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
528 fact.set (R);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
529 int err = 0;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
530
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
531 if (down)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
532 err = fact.downdate (u);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
533 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
534 fact.update (u);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
535
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
536 if (nargout > 1)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
537 retval(1) = err;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
538 else if (err)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
539 error ("cholupdate: downdate violates positiveness");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
540
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
541 retval(0) = fact.chol_matrix ();
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
542 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
543 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
544 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
545 error ("cholupdate: dimension mismatch");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
546 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
547 error ("cholupdate: op must be \"+\" or \"-\"");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
548 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
549 else
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
550 print_usage ();
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
551
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
552 return retval;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
553 }
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
554
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
555 /*
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
556 %!test
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
557 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
558 %! -0.131721 0.738529 0.019851 -0.140295 ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
559 %! 0.124120 0.019851 0.354879 -0.059472 ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
560 %! -0.061673 -0.140295 -0.059472 0.600939 ];
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
561 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
562 %! u = [ 0.98950 ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
563 %! 0.39844 ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
564 %! 0.63484 ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
565 %! 0.13351 ];
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
566 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
567 %! R = chol(A);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
568 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
569 %! R1 = cholupdate(R,u);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
570 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
571 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
572 %! 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
573 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
574 %! R1 = cholupdate(R1,u,"-");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
575 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
576 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
577 %! assert(norm(R1 - R,Inf) < 1e1*eps)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
578 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
579 %!test
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
580 %! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
581 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
582 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
583 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ];
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
584 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
585 %! u = [ 0.54267 + 0.91519i ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
586 %! 0.99647 + 0.43141i ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
587 %! 0.83760 + 0.68977i ;
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
588 %! 0.39160 + 0.90378i ];
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
589 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
590 %! R = chol(A);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
591 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
592 %! R1 = cholupdate(R,u);
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
593 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
594 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
595 %! 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
596 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
597 %! R1 = cholupdate(R1,u,"-");
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
598 %!
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
599 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
600 %! assert(norm(R1 - R,Inf) < 1e1*eps)
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
601 */
40574114c514 implement Cholesky factorization updating
Jaroslav Hajek <highegg@gmail.com>
parents: 7515
diff changeset
602
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
603 /*
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
604 ;;; Local Variables: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
605 ;;; mode: C++ ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
606 ;;; End: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
607 */
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
608