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

annotate src/DLD-FUNCTIONS/chol.cc @ 13953:642e43164af6

fix behavior of chol (..., 'lower') to be compatible with matlab * chol.cc: transpose input matrix prior to factorization when chol (..., 'lower') is invoked so that only the lower triangular part is used.

author	Carlo de Falco <kingcrimson@tiscali.it>
date	Mon, 28 Nov 2011 12:39:39 +0100
parents	5fa482628bf6
children	b0cdd60db5e5

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

Mercurial > octave-nkf

annotate src/DLD-FUNCTIONS/chol.cc @ 13953:642e43164af6