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

annotate src/DLD-FUNCTIONS/chol.cc @ 7789:82be108cc558

First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code

author	David Bateman <dbateman@free.fr>
date	Sun, 27 Apr 2008 22:34:17 +0200
parents	efccca5f2ad7
children	87865ed7405f

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
7700 efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	24 // The cholupdate, cholinsert, choldelete and cholshift functions were
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	25 // written by Jaroslav Hajek <highegg@gmail.com>, Copyright (C) 2008
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	26 // VZLU Prague, a.s., Czech Republic.
7554 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"
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	34 #include "fCmplxCHOL.h"
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	35 #include "floatCHOL.h"
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	36 #include "SparseCmplxCHOL.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	37 #include "SparsedbleCHOL.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	38 #include "oct-spparms.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	39 #include "sparse-util.h"
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	40
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	41 #include "ov-re-sparse.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	42 #include "ov-cx-sparse.h"
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	43 #include "defun-dld.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	44 #include "error.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	45 #include "gripes.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	46 #include "oct-obj.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	47 #include "utils.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	48
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	49 DEFUN_DLD (chol, args, nargout,
3548 ab7fa5a8f23f [project @ 2000-02-03 01:17:15 by jwe] jwe parents: 3372 diff changeset	50 "-- texinfo --\n\
7650 eb7bdde776f2 Texinfo fixes John W. Eaton <jwe@octave.org> parents: 7559 diff changeset	51 @deftypefn {Loadable Function} {@var{r} =} chol (@var{a})\n\
eb7bdde776f2 Texinfo fixes John W. Eaton <jwe@octave.org> parents: 7559 diff changeset	52 @deftypefnx {Loadable Function} {[@var{r}, @var{p}] =} chol (@var{a})\n\
eb7bdde776f2 Texinfo fixes John W. Eaton <jwe@octave.org> parents: 7559 diff changeset	53 @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	54 @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	55 @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	56 @cindex Cholesky factorization\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	57 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	58 matrix @var{a}, where\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	59 @iftex\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	60 @tex\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	61 $ R^T R = A $.\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	62 @end tex\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	63 @end iftex\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	64 @ifinfo\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	65 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	66 @example\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	67 @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	68 @end example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	69 @end ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	70 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	71 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	72 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	73 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	74 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	75 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	76 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	77 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	78 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	79 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	80 @iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	81 @tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	82 $ 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	83 @end tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	84 @end iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	85 @ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	86 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	87 @example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	88 @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	89 @end example\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	90 @end ifinfo\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	91 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	92 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	93 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	94 such that\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	95 @iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	96 @tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	97 $ 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	98 @end tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	99 @end iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	100 @ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	101 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	102 @example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	103 @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	104 @end example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	105 @end ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	106 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	107 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	108 @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	109 @iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	110 @tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	111 $ L L^T = A $.\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	112 @end tex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	113 @end iftex\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	114 @ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	115 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	116 @example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	117 @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	118 @end example\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	119 @end ifinfo\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	120 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	121 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	122 sparse matrices.\n\
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	123 @seealso{cholinv, chol2inv}\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3244 diff changeset	124 @end deftypefn")
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	125 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	126 octave_value_list retval;
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	127 int nargin = args.length ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	128 bool LLt = false;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	129 bool vecout = false;
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	130
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	131 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	132 \|\| (! args(0).is_sparse_type () && nargout > 2))
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	133 {
5823 080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe] jwe parents: 5760 diff changeset	134 print_usage ();
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	135 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	136 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	137
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	138 int n = 1;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	139 while (n < nargin && ! error_state)
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 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	142
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	143 if (! error_state )
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	144 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	145 if (tmp.compare ("vector") == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	146 vecout = true;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	147 else if (tmp.compare ("lower") == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	148 LLt = true;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	149 else if (tmp.compare ("upper") == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	150 LLt = false;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	151 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	152 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	153 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	154 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	155 error ("chol: expecting trailing string arguments");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	156 }
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	157
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	158 if (! error_state)
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_value arg = args(0);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	161
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	162 octave_idx_type nr = arg.rows ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	163 octave_idx_type nc = arg.columns ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	164 bool natural = (nargout != 3);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	165
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	166 int arg_is_empty = empty_arg ("chol", nr, nc);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	167
7515 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 retval;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	170 if (arg_is_empty > 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	171 return octave_value (Matrix ());
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_sparse_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 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	176 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	177 SparseMatrix m = arg.sparse_matrix_value ();
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	178
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	179 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	180 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	181 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	182 SparseCHOL fact (m, info, natural);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	183 if (nargout == 3)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	184 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	185 retval(2) = fact.perm ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	186 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	187 retval(2) = fact.Q();
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	188
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	189 if (nargout > 1 \|\| info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	190 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	191 retval(1) = fact.P();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	192 if (LLt)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	193 retval(0) = fact.L();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	194 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	195 retval(0) = fact.R();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	196 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	197 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	198 error ("chol: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	199 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	200 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	201 else if (arg.is_complex_type ())
3243 dd00769643ae [project @ 1999-05-28 04:19:00 by jwe] jwe parents: 3014 diff changeset	202 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	203 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	204
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	205 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	206 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	207 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	208 SparseComplexCHOL fact (m, info, natural);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	209
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	210 if (nargout == 3)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	211 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	212 retval(2) = fact.perm ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	213 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	214 retval(2) = fact.Q();
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 if (nargout > 1 \|\| info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	217 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	218 retval(1) = fact.P();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	219 if (LLt)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	220 retval(0) = fact.L();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	221 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	222 retval(0) = fact.R();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	223 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	224 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	225 error ("chol: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	226 }
3243 dd00769643ae [project @ 1999-05-28 04:19:00 by jwe] jwe parents: 3014 diff changeset	227 }
3244 7f8827f4040a [project @ 1999-05-29 14:01:55 by jwe] jwe parents: 3243 diff changeset	228 else
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	229 gripe_wrong_type_arg ("chol", arg);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	230 }
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	231 else if (arg.is_single_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	232 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	233 if (arg.is_real_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	234 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	235 FloatMatrix m = arg.float_matrix_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	236
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	237 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	238 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	239 octave_idx_type info;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	240 FloatCHOL fact (m, info);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	241 if (nargout == 2 \|\| info == 0)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	242 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	243 retval(1) = static_cast<float> (info);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	244 if (LLt)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	245 retval(0) = fact.chol_matrix ().transpose ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	246 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	247 retval(0) = fact.chol_matrix ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	248 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	249 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	250 error ("chol: matrix not positive definite");
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	251 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	252 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	253 else if (arg.is_complex_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	254 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	255 FloatComplexMatrix m = arg.float_complex_matrix_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	256
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	257 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	258 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	259 octave_idx_type info;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	260 FloatComplexCHOL fact (m, info);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	261 if (nargout == 2 \|\| info == 0)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	262 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	263 retval(1) = static_cast<float> (info);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	264 if (LLt)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	265 retval(0) = fact.chol_matrix ().hermitian ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	266 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	267 retval(0) = fact.chol_matrix ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	268 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	269 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	270 error ("chol: matrix not positive definite");
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	271 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	272 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	273 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	274 gripe_wrong_type_arg ("chol", arg);
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7700 diff changeset	275 }
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	276 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	277 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	278 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	279 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	280 Matrix m = arg.matrix_value ();
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	281
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	282 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	283 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	284 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	285 CHOL fact (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	286 if (nargout == 2 \|\| info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	287 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	288 retval(1) = static_cast<double> (info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	289 if (LLt)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	290 retval(0) = fact.chol_matrix ().transpose ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	291 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	292 retval(0) = fact.chol_matrix ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	293 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	294 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	295 error ("chol: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	296 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	297 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	298 else if (arg.is_complex_type ())
3243 dd00769643ae [project @ 1999-05-28 04:19:00 by jwe] jwe parents: 3014 diff changeset	299 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	300 ComplexMatrix m = arg.complex_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	301
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	302 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	303 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	304 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	305 ComplexCHOL fact (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	306 if (nargout == 2 \|\| info == 0)
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 retval(1) = static_cast<double> (info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	309 if (LLt)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	310 retval(0) = fact.chol_matrix ().hermitian ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	311 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	312 retval(0) = fact.chol_matrix ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	313 }
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 ("chol: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	316 }
3243 dd00769643ae [project @ 1999-05-28 04:19:00 by jwe] jwe parents: 3014 diff changeset	317 }
3244 7f8827f4040a [project @ 1999-05-29 14:01:55 by jwe] jwe parents: 3243 diff changeset	318 else
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	319 gripe_wrong_type_arg ("chol", arg);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	320 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	321 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	322
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	323 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	324 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	325
5760 8d7162924bd3 [project @ 2006-04-14 04:01:37 by jwe] jwe parents: 5448 diff changeset	326 DEFUN_DLD (cholinv, args, ,
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	327 "-- texinfo --\n\
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	328 @deftypefn {Loadable Function} {} cholinv (@var{a})\n\
5448 ebe5d7d15522 [project @ 2005-09-14 18:55:04 by jwe] jwe parents: 5343 diff changeset	329 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	330 symmetric positive definite matrix @var{a}.\n\
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	331 @seealso{chol, chol2inv}\n\
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	332 @end deftypefn")
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	333 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	334 octave_value retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	335
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	336 int nargin = args.length ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	337
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	338 if (nargin == 1)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	339 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	340 octave_value arg = args(0);
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	341
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	342 octave_idx_type nr = arg.rows ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	343 octave_idx_type nc = arg.columns ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	344
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	345 if (nr == 0 \|\| nc == 0)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	346 retval = Matrix ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	347 else
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	348 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	349 if (arg.is_sparse_type ())
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	350 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	351 if (arg.is_real_type ())
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	352 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	353 SparseMatrix m = arg.sparse_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 SparseCHOL 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 }
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	364 }
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	365 else if (arg.is_complex_type ())
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	366 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	367 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	368
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	369 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	370 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	371 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	372 SparseComplexCHOL chol (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	373 if (info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	374 retval = chol.inverse ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	375 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	376 error ("cholinv: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	377 }
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	378 }
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	379 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	380 gripe_wrong_type_arg ("cholinv", arg);
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	381 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	382 else
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	383 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	384 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	385 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	386 Matrix m = arg.matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	387
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	388 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	389 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	390 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	391 CHOL chol (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	392 if (info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	393 retval = chol.inverse ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	394 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	395 error ("cholinv: matrix not positive definite");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	396 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	397 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	398 else if (arg.is_complex_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	399 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	400 ComplexMatrix m = arg.complex_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	401
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	402 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	403 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	404 octave_idx_type info;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	405 ComplexCHOL chol (m, info);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	406 if (info == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	407 retval = chol.inverse ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	408 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	409 error ("cholinv: matrix not positive definite");
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 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	412 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	413 gripe_wrong_type_arg ("chol", arg);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	414 }
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	415 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	416 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	417 else
5823 080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe] jwe parents: 5760 diff changeset	418 print_usage ();
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	419
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	420 return retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	421 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	422
5760 8d7162924bd3 [project @ 2006-04-14 04:01:37 by jwe] jwe parents: 5448 diff changeset	423 DEFUN_DLD (chol2inv, args, ,
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	424 "-- texinfo --\n\
5343 be0bd7a760c9 [project @ 2005-05-06 18:04:16 by jwe] jwe parents: 5340 diff changeset	425 @deftypefn {Loadable Function} {} chol2inv (@var{u})\n\
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	426 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	427 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	428 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	429 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	430 using @code{inv}.\n\
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	431 @seealso{chol, cholinv}\n\
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	432 @end deftypefn")
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	433 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	434 octave_value retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	435
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	436 int nargin = args.length ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	437
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	438 if (nargin == 1)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	439 {
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	440 octave_value arg = args(0);
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 octave_idx_type nr = arg.rows ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	443 octave_idx_type nc = arg.columns ();
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 if (nr == 0 \|\| nc == 0)
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	446 retval = Matrix ();
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	447 else
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	448 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	449 if (arg.is_sparse_type ())
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	450 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	451 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	452 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	453 SparseMatrix r = arg.sparse_matrix_value ();
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	454
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	455 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	456 retval = chol2inv (r);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	457 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	458 else if (arg.is_complex_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	459 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	460 SparseComplexMatrix r = arg.sparse_complex_matrix_value ();
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	461
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	462 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	463 retval = chol2inv (r);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	464 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	465 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	466 gripe_wrong_type_arg ("chol2inv", arg);
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	467 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	468 else
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	469 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	470 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	471 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	472 Matrix r = arg.matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	473
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	474 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	475 retval = chol2inv (r);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	476 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	477 else if (arg.is_complex_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	478 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	479 ComplexMatrix r = arg.complex_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	480
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	481 if (! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	482 retval = chol2inv (r);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	483 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	484 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	485 gripe_wrong_type_arg ("chol2inv", arg);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	486 }
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	487 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	488 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	489 else
5823 080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe] jwe parents: 5760 diff changeset	490 print_usage ();
5340 15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	491
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	492 return retval;
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	493 }
15843d76156d [project @ 2005-05-06 16:26:58 by jwe] jwe parents: 5307 diff changeset	494
7554 40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	495 DEFUN_DLD (cholupdate, args, nargout,
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	496 "-- texinfo --\n\
7650 eb7bdde776f2 Texinfo fixes John W. Eaton <jwe@octave.org> parents: 7559 diff changeset	497 @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	498 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	499 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	500 upper triangular matrix @var{R1} such that\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	501 @itemize @bullet\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	502 @item\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	503 @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	504 if @var{op} is \"+\"\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	505 @item\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	506 @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	507 if @var{op} is \"-\"\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	508 @end itemize\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	509 \n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	510 If @var{op} is \"-\", @var{info} is set to\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	511 @itemize\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	512 @item 0 if the downdate was successful,\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	513 @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	514 @item 2 if @var{R} is singular.\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	515 @end itemize\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	516 \n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	517 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	518 @seealso{chol, qrupdate}\n\
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	519 @end deftypefn")
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	520 {
7559 07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	521 octave_idx_type nargin = args.length ();
7554 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 octave_value_list retval;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	524
7559 07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	525 if (nargin > 3 \|\| nargin < 2)
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	526 {
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	527 print_usage ();
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	528 return retval;
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	529 }
7554 40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	530
7559 07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	531 octave_value argr = args(0);
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	532 octave_value argu = args(1);
7554 40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	533
7559 07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	534 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	535 && (nargin < 3 \|\| args(2).is_string ()))
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	536 {
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	537 octave_idx_type n = argr.rows ();
7554 40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	538
7559 07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	539 std::string op = (nargin < 3) ? "+" : args(2).string_value ();
7554 40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	540
7559 07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	541 bool down = op == "-";
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	542
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	543 if (down \|\| op == "+")
07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	544 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	545 {
7559 07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	546 if (argr.is_real_matrix () && argu.is_real_matrix ())
7554 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 // real case
7559 07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	549 Matrix R = argr.matrix_value ();
7554 40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	550 Matrix u = argu.matrix_value ();
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 CHOL fact;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	553 fact.set (R);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	554 int err = 0;
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 if (down)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	557 err = fact.downdate (u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	558 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	559 fact.update (u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	560
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	561 if (nargout > 1)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	562 retval(1) = err;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	563 else if (err)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	564 error ("cholupdate: downdate violates positiveness");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	565
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	566 retval(0) = fact.chol_matrix ();
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	567 }
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	568 else
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 // complex case
7559 07522d7dcdf8 fixes to QR and Cholesky updating code Jaroslav Hajek <highegg@gmail.com> parents: 7554 diff changeset	571 ComplexMatrix R = argr.complex_matrix_value ();
7554 40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	572 ComplexMatrix u = argu.complex_matrix_value ();
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 ComplexCHOL fact;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	575 fact.set (R);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	576 int err = 0;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	577
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	578 if (down)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	579 err = fact.downdate (u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	580 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	581 fact.update (u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	582
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	583 if (nargout > 1)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	584 retval(1) = err;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	585 else if (err)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	586 error ("cholupdate: downdate violates positiveness");
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 retval(0) = fact.chol_matrix ();
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 }
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	591 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	592 error ("cholupdate: dimension mismatch");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	593 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	594 error ("cholupdate: op must be \"+\" or \"-\"");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	595 }
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	596 else
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	597 print_usage ();
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 return retval;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	600 }
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 /*
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	603 %!test
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	604 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	605 %! -0.131721 0.738529 0.019851 -0.140295 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	606 %! 0.124120 0.019851 0.354879 -0.059472 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	607 %! -0.061673 -0.140295 -0.059472 0.600939 ];
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	608 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	609 %! u = [ 0.98950 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	610 %! 0.39844 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	611 %! 0.63484 ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	612 %! 0.13351 ];
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	613 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	614 %! R = chol(A);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	615 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	616 %! R1 = cholupdate(R,u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	617 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	618 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	619 %! assert(norm(R1'R1 - R'R - uu',Inf) < 1e1eps)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	620 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	621 %! R1 = cholupdate(R1,u,"-");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	622 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	623 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	624 %! assert(norm(R1 - R,Inf) < 1e1*eps)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	625 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	626 %!test
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	627 %! 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	628 %! -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	629 %! 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	630 %! -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	631 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	632 %! u = [ 0.54267 + 0.91519i ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	633 %! 0.99647 + 0.43141i ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	634 %! 0.83760 + 0.68977i ;
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	635 %! 0.39160 + 0.90378i ];
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	636 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	637 %! R = chol(A);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	638 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	639 %! R1 = cholupdate(R,u);
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	640 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	641 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	642 %! assert(norm(R1'R1 - R'R - uu',Inf) < 1e1eps)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	643 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	644 %! R1 = cholupdate(R1,u,"-");
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	645 %!
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	646 %! assert(norm(triu(R1)-R1,Inf) == 0)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	647 %! assert(norm(R1 - R,Inf) < 1e1*eps)
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	648 */
40574114c514 implement Cholesky factorization updating Jaroslav Hajek <highegg@gmail.com> parents: 7515 diff changeset	649
7700 efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	650 DEFUN_DLD (cholinsert, args, nargout,
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	651 "-- texinfo --\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	652 @deftypefn {Loadable Function} {[@var{R1}, @var{info}] =} cholinsert (@var{R}, @var{j}, @var{u})\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	653 Given a Cholesky@tie{}factorization of a real symmetric or complex hermitian\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	654 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper triangular,\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	655 return the QR@tie{}factorization of\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	656 @var{A1}, where @w{A1(p,p) = A}, @w{A1(:,j) = A1(j,:)' = u} and\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	657 @w{p = [1:j-1,j+1:n+1]}. @w{u(j)} should be positive.\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	658 On return, @var{info} is set to\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	659 @itemize\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	660 @item 0 if the insertion was successful,\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	661 @item 1 if @var{A1} is not positive definite,\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	662 @item 2 if @var{R} is singular.\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	663 @end itemize\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	664 \n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	665 If @var{info} is not present, an error message is printed in cases 1 and 2.\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	666 @seealso{chol, cholupdate, choldelete}\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	667 @end deftypefn")
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	668 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	669 octave_idx_type nargin = args.length ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	670
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	671 octave_value_list retval;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	672
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	673 if (nargin != 3)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	674 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	675 print_usage ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	676 return retval;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	677 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	678
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	679 octave_value argr = args(0);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	680 octave_value argj = args(1);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	681 octave_value argu = args(2);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	682
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	683 if (argr.is_matrix_type () && argu.is_matrix_type ()
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	684 && argj.is_real_scalar ())
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	685 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	686 octave_idx_type n = argr.rows ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	687 octave_idx_type j = argj.scalar_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	688
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	689 if (argr.columns () == n && argu.rows () == n+1 && argu.columns () == 1)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	690 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	691 if (j > 0 && j <= n+1)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	692 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	693 if (argr.is_real_matrix () && argu.is_real_matrix ())
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	694 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	695 // real case
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	696 Matrix R = argr.matrix_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	697 Matrix u = argu.matrix_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	698
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	699 CHOL fact;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	700 fact.set (R);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	701 int err = fact.insert_sym (u, j-1);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	702
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	703 if (nargout > 1)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	704 retval(1) = err;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	705 else if (err)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	706 error ("cholinsert: insertion violates positiveness");
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	707
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	708 retval(0) = fact.chol_matrix ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	709 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	710 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	711 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	712 // complex case
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	713 ComplexMatrix R = argr.complex_matrix_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	714 ComplexMatrix u = argu.complex_matrix_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	715
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	716 ComplexCHOL fact;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	717 fact.set (R);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	718 int err = fact.insert_sym (u, j-1);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	719
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	720 if (nargout > 1)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	721 retval(1) = err;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	722 else if (err)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	723 error ("cholinsert: insertion violates positiveness");
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	724
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	725 retval(0) = fact.chol_matrix ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	726 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	727 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	728 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	729 error ("cholinsert: index out of range");
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	730 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	731 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	732 error ("cholinsert: dimension mismatch");
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	733 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	734 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	735 print_usage ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	736
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	737 return retval;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	738 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	739
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	740 /*
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	741 %!test
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	742 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	743 %! -0.131721 0.738529 0.019851 -0.140295 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	744 %! 0.124120 0.019851 0.354879 -0.059472 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	745 %! -0.061673 -0.140295 -0.059472 0.600939 ];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	746 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	747 %! u = [ 0.35080 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	748 %! 0.63930 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	749 %! 3.31057 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	750 %! -0.13825 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	751 %! 0.45266 ];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	752 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	753 %! R = chol(A);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	754 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	755 %! j = 3; p = [1:j-1, j+1:5];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	756 %! R1 = cholinsert(R,j,u); A1 = R1'*R1;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	757 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	758 %! assert(norm(triu(R1)-R1,Inf) == 0)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	759 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	760 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	761 %!test
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	762 %! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	763 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	764 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	765 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	766 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	767 %! u = [ 0.35080 + 0.04298i;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	768 %! 0.63930 + 0.23778i;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	769 %! 3.31057 + 0.00000i;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	770 %! -0.13825 + 0.19879i;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	771 %! 0.45266 + 0.50020i];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	772 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	773 %! R = chol(A);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	774 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	775 %! j = 3; p = [1:j-1, j+1:5];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	776 %! R1 = cholinsert(R,j,u); A1 = R1'*R1;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	777 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	778 %! assert(norm(triu(R1)-R1,Inf) == 0)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	779 %! assert(norm(A1(p,p) - A,Inf) < 1e1*eps)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	780 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	781 */
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	782
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	783 DEFUN_DLD (choldelete, args, nargout,
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	784 "-- texinfo --\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	785 @deftypefn {Loadable Function} {@var{R1} =} choldelete (@var{R}, @var{j})\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	786 Given a Cholesky@tie{}factorization of a real symmetric or complex hermitian\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	787 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper triangular,\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	788 return the QR@tie{}factorization of @w{A(p,p)}, where @w{p = [1:j-1,j+1:n+1]}.\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	789 @seealso{chol, cholupdate, cholinsert}\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	790 @end deftypefn")
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	791 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	792 octave_idx_type nargin = args.length ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	793
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	794 octave_value_list retval;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	795
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	796 if (nargin != 2)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	797 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	798 print_usage ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	799 return retval;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	800 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	801
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	802 octave_value argr = args(0);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	803 octave_value argj = args(1);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	804
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	805 if (argr.is_matrix_type () && argj.is_real_scalar ())
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	806 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	807 octave_idx_type n = argr.rows ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	808 octave_idx_type j = argj.scalar_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	809
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	810 if (argr.columns () == n)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	811 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	812 if (j > 0 && j <= n)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	813 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	814 if (argr.is_real_matrix ())
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	815 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	816 // real case
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	817 Matrix R = argr.matrix_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	818
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	819 CHOL fact;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	820 fact.set (R);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	821 fact.delete_sym (j-1);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	822
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	823 retval(0) = fact.chol_matrix ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	824 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	825 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	826 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	827 // complex case
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	828 ComplexMatrix R = argr.complex_matrix_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	829
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	830 ComplexCHOL fact;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	831 fact.set (R);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	832 fact.delete_sym (j-1);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	833
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	834 retval(0) = fact.chol_matrix ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	835 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	836 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	837 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	838 error ("choldelete: index out of range");
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	839 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	840 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	841 error ("choldelete: dimension mismatch");
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	842 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	843 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	844 print_usage ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	845
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	846 return retval;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	847 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	848
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	849 /*
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	850 %!test
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	851 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	852 %! -0.131721 0.738529 0.019851 -0.140295 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	853 %! 0.124120 0.019851 0.354879 -0.059472 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	854 %! -0.061673 -0.140295 -0.059472 0.600939 ];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	855 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	856 %! R = chol(A);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	857 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	858 %! j = 3; p = [1:j-1,j+1:4];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	859 %! R1 = choldelete(R,j);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	860 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	861 %! assert(norm(triu(R1)-R1,Inf) == 0)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	862 %! assert(norm(R1'R1 - A(p,p),Inf) < 1e1eps)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	863 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	864 %!test
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	865 %! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	866 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	867 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	868 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	869 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	870 %! R = chol(A);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	871 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	872 %! j = 3; p = [1:j-1,j+1:4];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	873 %! R1 = choldelete(R,j);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	874 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	875 %! assert(norm(triu(R1)-R1,Inf) == 0)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	876 %! assert(norm(R1'R1 - A(p,p),Inf) < 1e1eps)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	877 */
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	878
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	879 DEFUN_DLD (cholshift, args, nargout,
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	880 "-- texinfo --\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	881 @deftypefn {Loadable Function} {@var{R1} =} cholshift (@var{R}, @var{i}, @var{j})\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	882 Given a Cholesky@tie{}factorization of a real symmetric or complex hermitian\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	883 positive definite matrix @w{@var{A} = @var{R}'*@var{R}}, @var{R}@tie{}upper triangular,\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	884 return the QR@tie{}factorization of\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	885 @w{@var{A}(p,p)}, where @w{p} is the permutation @*\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	886 @code{p = [1:i-1, shift(i:j, 1), j+1:n]} if @w{@var{i} < @var{j}} @*\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	887 or @*\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	888 @code{p = [1:j-1, shift(j:i,-1), i+1:n]} if @w{@var{j} < @var{i}}. @*\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	889 \n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	890 @seealso{chol, cholinsert, choldelete}\n\
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	891 @end deftypefn")
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	892 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	893 octave_idx_type nargin = args.length ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	894
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	895 octave_value_list retval;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	896
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	897 if (nargin != 3)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	898 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	899 print_usage ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	900 return retval;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	901 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	902
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	903 octave_value argr = args(0);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	904 octave_value argi = args(1);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	905 octave_value argj = args(2);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	906
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	907 if (argr.is_matrix_type () && argi.is_real_scalar () && argj.is_real_scalar ())
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	908 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	909 octave_idx_type n = argr.rows ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	910 octave_idx_type i = argi.scalar_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	911 octave_idx_type j = argj.scalar_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	912
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	913 if (argr.columns () == n)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	914 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	915 if (j > 0 && j <= n+1 && i > 0 && i <= n+1)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	916 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	917 if (argr.is_real_matrix ())
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	918 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	919 // real case
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	920 Matrix R = argr.matrix_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	921
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	922 CHOL fact;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	923 fact.set (R);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	924 fact.shift_sym (i-1, j-1);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	925
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	926 retval(0) = fact.chol_matrix ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	927 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	928 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	929 {
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	930 // complex case
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	931 ComplexMatrix R = argr.complex_matrix_value ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	932
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	933 ComplexCHOL fact;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	934 fact.set (R);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	935 fact.shift_sym (i-1, j-1);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	936
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	937 retval(0) = fact.chol_matrix ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	938 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	939 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	940 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	941 error ("cholshift: index out of range");
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	942 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	943 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	944 error ("cholshift: dimension mismatch");
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	945 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	946 else
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	947 print_usage ();
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	948
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	949 return retval;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	950 }
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	951
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	952 /*
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	953 %!test
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	954 %! A = [ 0.436997 -0.131721 0.124120 -0.061673 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	955 %! -0.131721 0.738529 0.019851 -0.140295 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	956 %! 0.124120 0.019851 0.354879 -0.059472 ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	957 %! -0.061673 -0.140295 -0.059472 0.600939 ];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	958 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	959 %! R = chol(A);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	960 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	961 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	962 %! R1 = cholshift(R,i,j);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	963 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	964 %! assert(norm(triu(R1)-R1,Inf) == 0)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	965 %! assert(norm(R1'R1 - A(p,p),Inf) < 1e1eps)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	966 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	967 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	968 %! R1 = cholshift(R,i,j);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	969 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	970 %! assert(norm(triu(R1)-R1,Inf) == 0)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	971 %! assert(norm(R1'R1 - A(p,p),Inf) < 1e1eps)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	972 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	973 %!test
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	974 %! A = [ 0.5585528 + 0.0000000i -0.1662088 - 0.0315341i 0.0107873 + 0.0236411i -0.0276775 - 0.0186073i ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	975 %! -0.1662088 + 0.0315341i 0.6760061 + 0.0000000i 0.0011452 - 0.0475528i 0.0145967 + 0.0247641i ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	976 %! 0.0107873 - 0.0236411i 0.0011452 + 0.0475528i 0.6263149 - 0.0000000i -0.1585837 - 0.0719763i ;
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	977 %! -0.0276775 + 0.0186073i 0.0145967 - 0.0247641i -0.1585837 + 0.0719763i 0.6034234 - 0.0000000i ];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	978 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	979 %! R = chol(A);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	980 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	981 %! i = 1; j = 3; p = [1:i-1, shift(i:j,-1), j+1:4];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	982 %! R1 = cholshift(R,i,j);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	983 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	984 %! assert(norm(triu(R1)-R1,Inf) == 0)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	985 %! assert(norm(R1'R1 - A(p,p),Inf) < 1e1eps)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	986 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	987 %! j = 1; i = 3; p = [1:j-1, shift(j:i,+1), i+1:4];
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	988 %! R1 = cholshift(R,i,j);
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	989 %!
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	990 %! assert(norm(triu(R1)-R1,Inf) == 0)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	991 %! assert(norm(R1'R1 - A(p,p),Inf) < 1e1eps)
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	992 */
efccca5f2ad7 more QR & Cholesky updating functions Jaroslav Hajek <highegg@gmail.com> parents: 7650 diff changeset	993
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	994 /*
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	995 ;;; Local Variables: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	996 ;;; mode: C++ ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	997 ;;; End: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	998 */
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	999

Mercurial > octave-nkf

annotate src/DLD-FUNCTIONS/chol.cc @ 7789:82be108cc558