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

annotate src/DLD-FUNCTIONS/chol.cc @ 8562:a6edd5c23cb5

use replacement methods if qrupdate is not available

author	Jaroslav Hajek <highegg@gmail.com>
date	Thu, 22 Jan 2009 11:10:47 +0100
parents	d66c9b6e506a
children	fb1c929dbbb7

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

Mercurial > octave-nkf

annotate src/DLD-FUNCTIONS/chol.cc @ 8562:a6edd5c23cb5