octave-nkf: src/DLD-FUNCTIONS/chol.cc annotate

annotate src/DLD-FUNCTIONS/chol.cc @ 7554:40574114c514

implement Cholesky factorization updating

author	Jaroslav Hajek <highegg@gmail.com>
date	Tue, 04 Mar 2008 22:25:50 -0500
parents	f3c00dc0912b
children	07522d7dcdf8

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\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	49 @deftypefn {Loadable Function} {@var{r}} = chol (@var{a})\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	50 @deftypefnx {Loadable Function} {[@var{r}, @var{p}]} = chol (@var{a})\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	51 @deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}]} = chol (@var{s})\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	52 @deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}]} = chol (@var{s}, 'vector')\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 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\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	450 @deftypefn {Loadable Function} {[@var{R1}, @var{info}]} = cholupdate (@var{R}, @var{u}, @var{op})\n\
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 {
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	474 int nargin = args.length ();
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
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	478 octave_value argR,argu,argop;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	479
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	480 if ((nargin == 3 \|\| nargin == 2)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	481 && (argR = args(0), argR.is_matrix_type ())
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	482 && (argu = args(1), argu.is_matrix_type ())
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	483 && (nargin < 3 \|\| (argop = args(2), argop.is_string ())))
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	484 {
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	485 octave_idx_type n = argR.rows ();
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	486
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	487 std::string op = (nargin < 3) ? "+" : argop.string_value();
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	488
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	489 bool down = false;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	490
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	491 if (nargin < 3 \|\| (op == "+") \|\| (down = op == "-"))
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	492 if (argR.columns () == n
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	493 && argu.rows () == n && argu.columns () == 1)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	494 {
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	495 if (argR.is_real_matrix () && argu.is_real_matrix ())
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	496 {
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	497 // real case
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	498 Matrix R = argR.matrix_value ();
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	499 Matrix u = argu.matrix_value ();
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 CHOL fact;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	502 fact.set (R);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	503 int err = 0;
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 if (down)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	506 err = fact.downdate (u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	507 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	508 fact.update (u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	509
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	510 if (nargout > 1)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	511 retval(1) = err;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	512 else if (err)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	513 error ("cholupdate: downdate violates positiveness");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	514
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	515 retval(0) = fact.chol_matrix ();
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	516 }
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	517 else
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 // complex case
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	520 ComplexMatrix R = argR.complex_matrix_value ();
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	521 ComplexMatrix u = argu.complex_matrix_value ();
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 ComplexCHOL fact;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	524 fact.set (R);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	525 int err = 0;
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 if (down)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	528 err = fact.downdate (u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	529 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	530 fact.update (u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	531
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	532 if (nargout > 1)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	533 retval(1) = err;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	534 else if (err)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	535 error ("cholupdate: downdate violates positiveness");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	536
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	537 retval(0) = fact.chol_matrix ();
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	538 }
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	539 }
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	540 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	541 error ("cholupdate: dimension mismatch");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	542 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	543 error ("cholupdate: op must be \"+\" or \"-\"");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	544 }
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	545 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	546 print_usage ();
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	547
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	548 return retval;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	549 }
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	550
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 %!test
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	553 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	554 %! -0.131721 0.738529 0.019851 -0.140295 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	555 %! 0.124120 0.019851 0.354879 -0.059472 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	556 %! -0.061673 -0.140295 -0.059472 0.600939 ];
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	557 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	558 %! u = [ 0.98950 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	559 %! 0.39844 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	560 %! 0.63484 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	561 %! 0.13351 ];
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	562 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	563 %! R = chol(A);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	564 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	565 %! R1 = cholupdate(R,u);
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 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	568 %! assert(norm(R1'R1 - R'R - uu',Inf) < 1e1eps)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	569 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	570 %! R1 = cholupdate(R1,u,"-");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	571 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	572 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	573 %! assert(norm(R1 - R,Inf) < 1e1*eps)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	574 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	575 %!test
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	576 %! 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	577 %! -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	578 %! 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	579 %! -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	580 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	581 %! u = [ 0.54267 + 0.91519i ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	582 %! 0.99647 + 0.43141i ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	583 %! 0.83760 + 0.68977i ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	584 %! 0.39160 + 0.90378i ];
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	585 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	586 %! R = chol(A);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	587 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	588 %! R1 = cholupdate(R,u);
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 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	591 %! assert(norm(R1'R1 - R'R - uu',Inf) < 1e1eps)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	592 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	593 %! R1 = cholupdate(R1,u,"-");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	594 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	595 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	596 %! assert(norm(R1 - R,Inf) < 1e1*eps)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	597 */
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	598
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	599 /*
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	600 ;;; Local Variables: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	601 ;;; mode: C++ ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	602 ;;; End: ***
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

Mercurial > hg > octave-nkf

annotate src/DLD-FUNCTIONS/chol.cc @ 7554:40574114c514