octave-lyh: src/DLD-FUNCTIONS/lu.cc annotate

annotate src/DLD-FUNCTIONS/lu.cc @ 10080:cf70ee43077c

add tests for LU updating

author	Jaroslav Hajek <highegg@gmail.com>
date	Mon, 11 Jan 2010 09:45:35 +0100
parents	f22bbc5d56e9
children	ab1101011a6d

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
8920 eb63fbe60fab update copyright notices John W. Eaton <jwe@octave.org> parents: 8869 diff changeset	3 Copyright (C) 1996, 1997, 1999, 2000, 2003, 2005, 2006, 2007, 2008, 2009
7017 a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe] jwe parents: 7016 diff changeset	4 John W. Eaton
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	5
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	6 This file is part of Octave.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	7
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	8 Octave is free software; you can redistribute it and/or modify it
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	9 under the terms of the GNU General Public License as published by the
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	10 Free Software Foundation; either version 3 of the License, or (at your
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	11 option) any later version.
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	12
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	13 Octave is distributed in the hope that it will be useful, but WITHOUT
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	16 for more details.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	17
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	18 You should have received a copy of the GNU General Public License
7016 93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	19 along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe] jwe parents: 5823 diff changeset	20 <http://www.gnu.org/licenses/>.
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	21
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	22 */
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	23
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	24 #ifdef HAVE_CONFIG_H
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	25 #include <config.h>
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	26 #endif
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	27
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	28 #include "CmplxLU.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	29 #include "dbleLU.h"
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	30 #include "fCmplxLU.h"
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	31 #include "floatLU.h"
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	32 #include "SparseCmplxLU.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	33 #include "SparsedbleLU.h"
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	34
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	35 #include "defun-dld.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	36 #include "error.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	37 #include "gripes.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	38 #include "oct-obj.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	39 #include "utils.h"
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	40 #include "ov-re-sparse.h"
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	41 #include "ov-cx-sparse.h"
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	42
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	43 template <class MT>
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	44 static octave_value
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	45 get_lu_l (const base_lu<MT>& fact)
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	46 {
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	47 MT L = fact.L ();
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	48 if (L.is_square ())
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	49 return octave_value (L, MatrixType (MatrixType::Lower));
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	50 else
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	51 return L;
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	52 }
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	53
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	54 template <class MT>
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	55 static octave_value
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	56 get_lu_u (const base_lu<MT>& fact)
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	57 {
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	58 MT U = fact.U ();
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	59 if (U.is_square () && fact.regular ())
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	60 return octave_value (U, MatrixType (MatrixType::Upper));
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	61 else
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	62 return U;
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	63 }
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	64
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	65 DEFUN_DLD (lu, args, nargout,
3548 ab7fa5a8f23f [project @ 2000-02-03 01:17:15 by jwe] jwe parents: 3372 diff changeset	66 "-- texinfo --\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	67 @deftypefn {Loadable Function} {[@var{l}, @var{u}, @var{p}] =} lu (@var{a})\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	68 @deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} lu (@var{s})\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	69 @deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}] =} lu (@var{s})\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	70 @deftypefnx {Loadable Function} {[@dots{}] =} lu (@var{s}, @var{thres})\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	71 @deftypefnx {Loadable Function} {@var{y} =} lu (@dots{})\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	72 @deftypefnx {Loadable Function} {[@dots{}] =} lu (@dots{}, 'vector')\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	73 @cindex LU decomposition\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	74 Compute the LU decomposition of @var{a}. If @var{a} is full subroutines from\n\
9209 923c7cb7f13f Simplify TeXinfo files by eliminating redundant @iftex followed by @tex construction. Rik <rdrider0-list@yahoo.com> parents: 9065 diff changeset	75 @sc{lapack} are used and if @var{a} is sparse then UMFPACK is used. The\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	76 result is returned in a permuted form, according to the optional return\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	77 value @var{p}. For example, given the matrix @code{a = [1, 2; 3, 4]},\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	78 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	79 @example\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	80 [l, u, p] = lu (a)\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	81 @end example\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	82 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	83 @noindent\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	84 returns\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	85 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	86 @example\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	87 @group\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	88 l =\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	89 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	90 1.00000 0.00000\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	91 0.33333 1.00000\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	92 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	93 u =\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	94 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	95 3.00000 4.00000\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	96 0.00000 0.66667\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	97 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	98 p =\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	99 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	100 0 1\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	101 1 0\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	102 @end group\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	103 @end example\n\
4329 d53c33d93440 [project @ 2003-02-18 20:00:48 by jwe] jwe parents: 3587 diff changeset	104 \n\
5457 c6dc1ccd83a9 [project @ 2005-09-19 19:23:35 by jwe] jwe parents: 5307 diff changeset	105 The matrix is not required to be square.\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	106 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	107 Called with two or three output arguments and a spare input matrix,\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	108 then @dfn{lu} does not attempt to perform sparsity preserving column\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	109 permutations. Called with a fourth output argument, the sparsity\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	110 preserving column transformation @var{Q} is returned, such that\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	111 @code{@var{p} * @var{a} * @var{q} = @var{l} * @var{u}}.\n\
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 Called with a fifth output argument and a sparse input matrix, then\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	114 @dfn{lu} attempts to use a scaling factor @var{r} on the input matrix\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	115 such that @code{@var{p} * (@var{r} \\ @var{a}) * @var{q} = @var{l} * @var{u}}.\n\
9065 8207b833557f Cleanup documentation for arith.texi, linalg.texi, nonlin.texi Rik <rdrider0-list@yahoo.com> parents: 9064 diff changeset	116 This typically leads to a sparser and more stable factorization.\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	117 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	118 An additional input argument @var{thres}, that defines the pivoting\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	119 threshold can be given. @var{thres} can be a scalar, in which case\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	120 it defines UMFPACK pivoting tolerance for both symmetric and unsymmetric\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	121 cases. If @var{thres} is a two element vector, then the first element\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	122 defines the pivoting tolerance for the unsymmetric UMFPACK pivoting\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	123 strategy and the second the symmetric strategy. By default, the values\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	124 defined by @code{spparms} are used and are by default @code{[0.1, 0.001]}.\n\
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 Given the string argument 'vector', @dfn{lu} returns the values of @var{p}\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	127 @var{q} as vector values, such that for full matrix, @code{@var{a}\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	128 (@var{p},:) = @var{l} * @var{u}}, and @code{@var{r}(@var{p},:) * @var{a}\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	129 (:, @var{q}) = @var{l} * @var{u}}.\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	130 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	131 With two output arguments, returns the permuted forms of the upper and\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	132 lower triangular matrices, such that @code{@var{a} = @var{l} * @var{u}}.\n\
9209 923c7cb7f13f Simplify TeXinfo files by eliminating redundant @iftex followed by @tex construction. Rik <rdrider0-list@yahoo.com> parents: 9065 diff changeset	133 With one output argument @var{y}, then the matrix returned by the @sc{lapack}\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	134 routines is returned. If the input matrix is sparse then the matrix @var{l}\n\
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	135 is embedded into @var{u} to give a return value similar to the full case.\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	136 For both full and sparse matrices, @dfn{lu} looses the permutation\n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	137 information.\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 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 issparse = (nargin > 0 && args(0).is_sparse_type ());
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	143 bool scale = (nargout == 5);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	144
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	145 if (nargin < 1 \|\| (issparse && (nargin > 3 \|\| nargout > 5))
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	146 \|\| (!issparse && (nargin > 2 \|\| nargout > 3)))
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: 5457 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 bool vecout = false;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	153 Matrix thres;
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 int n = 1;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	156 while (n < nargin && ! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	157 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	158 if (args (n).is_string ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	159 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	160 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	161
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	162 if (! error_state )
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	163 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	164 if (tmp.compare ("vector") == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	165 vecout = true;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	166 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	167 error ("lu: unrecognized string argument");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	168 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	169 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	170 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	171 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	172 Matrix tmp = args(n++).matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	173
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	174 if (! error_state )
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	175 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	176 if (!issparse)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	177 error ("lu: can not define pivoting threshold for full matrices");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	178 else if (tmp.nelem () == 1)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	179 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	180 thres.resize(1,2);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	181 thres(0) = tmp(0);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	182 thres(1) = tmp(0);
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 else if (tmp.nelem () == 2)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	185 thres = tmp;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	186 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	187 error ("lu: expecting 2 element vector for thres");
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	188 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	189 }
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
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	192 octave_value arg = args(0);
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	193
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4329 diff changeset	194 octave_idx_type nr = arg.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4329 diff changeset	195 octave_idx_type nc = arg.columns ();
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	196
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	197 int arg_is_empty = empty_arg ("lu", nr, nc);
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	198
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	199 if (issparse)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	200 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	201 if (arg_is_empty < 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	202 return retval;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	203 else if (arg_is_empty > 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	204 return octave_value_list (5, SparseMatrix ());
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	205
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	206 ColumnVector Qinit;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	207 if (nargout < 4)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	208 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	209 Qinit.resize (nc);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	210 for (octave_idx_type i = 0; i < nc; i++)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	211 Qinit (i) = i;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	212 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	213
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	214 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	215 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	216 SparseMatrix m = arg.sparse_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	217
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	218 switch (nargout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	219 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	220 case 0:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	221 case 1:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	222 case 2:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	223 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	224 SparseLU fact (m, Qinit, thres, false, true);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	225
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	226 if (nargout < 2)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	227 retval (0) = fact.Y ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	228 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	229 {
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	230 PermMatrix P = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	231 SparseMatrix L = P.transpose () * fact.L ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	232 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	233 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	234
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	235 retval(0) = octave_value (L,
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	236 MatrixType (MatrixType::Permuted_Lower,
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	237 nr, fact.row_perm ()));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	238 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	239 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	240 break;
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	241
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	242 case 3:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	243 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	244 SparseLU fact (m, Qinit, thres, false, true);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	245
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	246 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	247 retval (2) = fact.Pr_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	248 else
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	249 retval(2) = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	250
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	251 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	252 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	253 retval(0) = octave_value (fact.L (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	254 MatrixType (MatrixType::Lower));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	255 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	256 break;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	257
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	258 case 4:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	259 default:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	260 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	261 SparseLU fact (m, thres, scale);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	262
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	263 if (scale)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	264 retval(4) = fact.R ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	265
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	266 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	267 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	268 retval(3) = fact.Pc_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	269 retval(2) = fact.Pr_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	270 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	271 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	272 {
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	273 retval(3) = fact.Pc_mat ();
3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	274 retval(2) = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	275 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	276 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	277 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	278 retval(0) = octave_value (fact.L (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	279 MatrixType (MatrixType::Lower));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	280 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	281 break;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	282 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	283 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	284 else if (arg.is_complex_type ())
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	285 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	286 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	287
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	288 switch (nargout)
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	289 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	290 case 0:
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	291 case 1:
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	292 case 2:
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	293 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	294 SparseComplexLU fact (m, Qinit, thres, false, true);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	295
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	296 if (nargout < 2)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	297 retval (0) = fact.Y ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	298 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	299 {
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	300 PermMatrix P = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	301 SparseComplexMatrix L = P.transpose () * fact.L ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	302 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	303 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	304
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	305 retval(0) = octave_value (L,
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	306 MatrixType (MatrixType::Permuted_Lower,
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	307 nr, fact.row_perm ()));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	308 }
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	309 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	310 break;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	311
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	312 case 3:
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	313 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	314 SparseComplexLU fact (m, Qinit, thres, false, true);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	315
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	316 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	317 retval (2) = fact.Pr_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	318 else
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	319 retval(2) = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	320
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	321 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	322 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	323 retval(0) = octave_value (fact.L (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	324 MatrixType (MatrixType::Lower));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	325 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	326 break;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	327
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	328 case 4:
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	329 default:
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	330 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	331 SparseComplexLU fact (m, thres, scale);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	332
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	333 if (scale)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	334 retval(4) = fact.R ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	335
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	336 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	337 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	338 retval(3) = fact.Pc_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	339 retval(2) = fact.Pr_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	340 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	341 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	342 {
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	343 retval(3) = fact.Pc_mat ();
3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	344 retval(2) = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	345 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	346 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	347 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	348 retval(0) = octave_value (fact.L (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	349 MatrixType (MatrixType::Lower));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	350 }
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	351 break;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	352 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	353 }
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	354 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	355 gripe_wrong_type_arg ("lu", arg);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	356 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	357 else
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	358 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	359 if (arg_is_empty < 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	360 return retval;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	361 else if (arg_is_empty > 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	362 return octave_value_list (3, Matrix ());
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	363
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	364 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	365 {
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	366 if (arg.is_single_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	367 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	368 FloatMatrix m = arg.float_matrix_value ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	369
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	370 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	371 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	372 FloatLU fact (m);
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	373
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	374 switch (nargout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	375 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	376 case 0:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	377 case 1:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	378 retval(0) = fact.Y ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	379 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	380
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	381 case 2:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	382 {
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	383 PermMatrix P = fact.P ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	384 FloatMatrix L = P.transpose () * fact.L ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	385 retval(1) = get_lu_u (fact);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	386 retval(0) = L;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	387 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	388 break;
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	389
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	390 case 3:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	391 default:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	392 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	393 if (vecout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	394 retval(2) = fact.P_vec ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	395 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	396 retval(2) = fact.P ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	397 retval(1) = get_lu_u (fact);
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	398 retval(0) = get_lu_l (fact);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	399 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	400 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	401 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	402 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	403 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	404 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	405 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	406 Matrix m = arg.matrix_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	407
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	408 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	409 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	410 LU fact (m);
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	411
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	412 switch (nargout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	413 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	414 case 0:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	415 case 1:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	416 retval(0) = fact.Y ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	417 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	418
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	419 case 2:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	420 {
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	421 PermMatrix P = fact.P ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	422 Matrix L = P.transpose () * fact.L ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	423 retval(1) = get_lu_u (fact);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	424 retval(0) = L;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	425 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	426 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	427
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	428 case 3:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	429 default:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	430 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	431 if (vecout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	432 retval(2) = fact.P_vec ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	433 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	434 retval(2) = fact.P ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	435 retval(1) = get_lu_u (fact);
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	436 retval(0) = get_lu_l (fact);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	437 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	438 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	439 }
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	440 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	441 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	442 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	443 else if (arg.is_complex_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	444 {
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	445 if (arg.is_single_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	446 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	447 FloatComplexMatrix m = arg.float_complex_matrix_value ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	448
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	449 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	450 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	451 FloatComplexLU fact (m);
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	452
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	453 switch (nargout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	454 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	455 case 0:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	456 case 1:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	457 retval(0) = fact.Y ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	458 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	459
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	460 case 2:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	461 {
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	462 PermMatrix P = fact.P ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	463 FloatComplexMatrix L = P.transpose () * fact.L ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	464 retval(1) = get_lu_u (fact);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	465 retval(0) = L;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	466 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	467 break;
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	468
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	469 case 3:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	470 default:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	471 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	472 if (vecout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	473 retval(2) = fact.P_vec ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	474 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	475 retval(2) = fact.P ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	476 retval(1) = get_lu_u (fact);
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	477 retval(0) = get_lu_l (fact);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	478 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	479 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	480 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	481 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	482 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	483 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	484 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	485 ComplexMatrix m = arg.complex_matrix_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	486
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	487 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	488 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	489 ComplexLU fact (m);
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	490
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	491 switch (nargout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	492 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	493 case 0:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	494 case 1:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	495 retval(0) = fact.Y ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	496 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	497
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	498 case 2:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	499 {
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	500 PermMatrix P = fact.P ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	501 ComplexMatrix L = P.transpose () * fact.L ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	502 retval(1) = get_lu_u (fact);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	503 retval(0) = L;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	504 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	505 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	506
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	507 case 3:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	508 default:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	509 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	510 if (vecout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	511 retval(2) = fact.P_vec ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	512 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	513 retval(2) = fact.P ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	514 retval(1) = get_lu_u (fact);
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	515 retval(0) = get_lu_l (fact);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	516 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	517 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	518 }
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	519 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	520 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	521 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	522 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	523 gripe_wrong_type_arg ("lu", arg);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	524 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	525
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	526 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	527 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	528
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	529 /*
7814 87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	530
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	531 %!assert(lu ([1, 2; 3, 4]), [3, 4; 1/3, 2/3], eps);
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	532
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	533 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	534 %! [l, u] = lu ([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	535 %! assert(l, [1/3, 1; 1, 0], sqrt (eps));
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	536 %! assert(u, [3, 4; 0, 2/3], sqrt (eps));
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	537
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	538 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	539 %! [l, u, p] = lu ([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	540 %! assert(l, [1, 0; 1/3, 1], sqrt (eps));
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	541 %! assert(u, [3, 4; 0, 2/3], sqrt (eps));
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	542 %! assert(p(:,:), [0, 1; 1, 0], sqrt (eps));
7814 87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	543
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	544 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	545 %! [l, u, p] = lu ([1, 2; 3, 4],'vector');
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	546 %! assert(l, [1, 0; 1/3, 1], sqrt (eps));
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	547 %! assert(u, [3, 4; 0, 2/3], sqrt (eps));
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	548 %! assert(p, [2;1], sqrt (eps));
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	549
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	550 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	551 %! [l u p] = lu ([1, 2; 3, 4; 5, 6]);
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	552 %! assert(l, [1, 0; 1/5, 1; 3/5, 1/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	553 %! assert(u, [5, 6; 0, 4/5], sqrt (eps));
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	554 %! assert(p(:,:), [0, 0, 1; 1, 0, 0; 0 1 0], sqrt (eps));
7814 87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	555
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	556 %!assert(lu (single([1, 2; 3, 4])), single([3, 4; 1/3, 2/3]), eps('single'));
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	557
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	558 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	559 %! [l, u] = lu (single([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	560 %! assert(l, single([1/3, 1; 1, 0]), 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	561 %! assert(u, single([3, 4; 0, 2/3]), 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	562
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	563 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	564 %! [l, u, p] = lu (single([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	565 %! assert(l, single([1, 0; 1/3, 1]), 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	566 %! assert(u, single([3, 4; 0, 2/3]), sqrt (eps('single')));
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	567 %! assert(p(:,:), single([0, 1; 1, 0]), sqrt (eps('single')));
7814 87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	568
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	569 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	570 %! [l, u, p] = lu (single([1, 2; 3, 4]),'vector');
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	571 %! assert(l, single([1, 0; 1/3, 1]), 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	572 %! assert(u, single([3, 4; 0, 2/3]), 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	573 %! assert(p, single([2;1]), 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	574
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	575 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	576 %! [l u p] = lu (single([1, 2; 3, 4; 5, 6]));
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	577 %! assert(l, single([1, 0; 1/5, 1; 3/5, 1/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	578 %! assert(u, single([5, 6; 0, 4/5]), sqrt (eps('single')));
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	579 %! assert(p(:,:), single([0, 0, 1; 1, 0, 0; 0 1 0]), sqrt (eps('single')));
7814 87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	580
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	581 %!error <Invalid call to lu.*> lu ();
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	582 %!error lu ([1, 2; 3, 4], 2);
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	583
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	584 */
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	585
9708 6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	586 static
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	587 bool check_lu_dims (const octave_value& l, const octave_value& u,
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	588 const octave_value& p)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	589 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	590 octave_idx_type m = l.rows (), k = u.rows (), n = u.columns ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	591 return ((l.ndims () == 2 && u.ndims () == 2 && k == l.columns ())
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	592 && k == std::min (m, n) &&
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	593 (p.is_undefined () \|\| p.rows () == m));
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	594 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	595
9709 f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	596 DEFUN_DLD (luupdate, args, ,
9708 6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	597 "-- texinfo --\n\
9724 f22bbc5d56e9 Fix various incorrect usages of TeXinfo deffn and deftypefn macros Rik <rdrider0-list@yahoo.com> parents: 9715 diff changeset	598 @deftypefn {Loadable Function} {[@var{L}, @var{U}] =} luupdate (@var{l}, @var{u}, @var{x}, @var{y})\n\
f22bbc5d56e9 Fix various incorrect usages of TeXinfo deffn and deftypefn macros Rik <rdrider0-list@yahoo.com> parents: 9715 diff changeset	599 @deftypefnx {Loadable Function} {[@var{L}, @var{U}, @var{P}] =} luupdate (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
9709 f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	600 Given an LU@tie{}factorization of a real or complex matrix\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	601 @w{@var{A} = @var{L}*@var{U}}, @var{L}@tie{}lower unit trapezoidal and\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	602 @var{U}@tie{}upper trapezoidal, return the LU@tie{}factorization\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	603 of @w{@var{A} + @var{x}*@var{y}.'}, where @var{x} and @var{y} are\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	604 column vectors (rank-1 update) or matrices with equal number of columns\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	605 (rank-k update).\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	606 Optionally, row-pivoted updating can be used by supplying\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	607 a row permutation (pivoting) matrix @var{P};\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	608 in that case, an updated permutation matrix is returned.\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	609 Note that if @var{L}, @var{U}, @var{P} is a pivoted LU@tie{}factorization\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	610 as obtained by @code{lu}:\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	611 \n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	612 @example\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	613 [@var{L}, @var{U}, @var{P}] = lu (@var{A});\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	614 @end example\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	615 \n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	616 then a factorization of @code{@var{a}+@var{x}*@var{y}.'} can be obtained either as\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	617 \n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	618 @example\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	619 [@var{L1}, @var{U1}] = lu (@var{L}, @var{U}, @var{P}*@var{x}, @var{y})\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	620 @end example\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	621 \n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	622 or\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	623 \n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	624 @example\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	625 [@var{L1}, @var{U1}, @var{P1}] = lu (@var{L}, @var{U}, @var{P}, @var{x}, @var{y})\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	626 @end example\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	627 \n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	628 The first form uses the unpivoted algorithm, which is faster, but less stable.\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	629 The second form uses a slower pivoted algorithm, which is more stable.\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	630 \n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	631 Note that the matrix case is done as a sequence of rank-1 updates;\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	632 thus, for k large enough, it will be both faster and more accurate to recompute\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	633 the factorization from scratch.\n\
f8e2e9fdaa8f document luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9708 diff changeset	634 @seealso{lu,qrupdate,cholupdate}\n\
9724 f22bbc5d56e9 Fix various incorrect usages of TeXinfo deffn and deftypefn macros Rik <rdrider0-list@yahoo.com> parents: 9715 diff changeset	635 @end deftypefn")
9708 6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	636 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	637 octave_idx_type nargin = args.length ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	638 octave_value_list retval;
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	639
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	640 bool pivoted = nargin == 5;
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	641
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	642 if (nargin != 4 && nargin != 5)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	643 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	644 print_usage ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	645 return retval;
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	646 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	647
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	648 octave_value argl = args(0);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	649 octave_value argu = args(1);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	650 octave_value argp = pivoted ? args(2) : octave_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	651 octave_value argx = args(2 + pivoted);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	652 octave_value argy = args(3 + pivoted);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	653
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	654 if (argl.is_numeric_type () && argu.is_numeric_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	655 && argx.is_numeric_type () && argy.is_numeric_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	656 && (! pivoted \|\| argp.is_perm_matrix ()))
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	657 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	658 if (check_lu_dims (argl, argu, argp))
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	659 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	660 PermMatrix P = (pivoted
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	661 ? argp.perm_matrix_value ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	662 : PermMatrix::eye (argl.rows ()));
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	663
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	664 if (argl.is_real_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	665 && argu.is_real_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	666 && argx.is_real_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	667 && argy.is_real_type ())
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	668 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	669 // all real case
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	670 if (argl.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	671 \|\| argu.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	672 \|\| argx.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	673 \|\| argy.is_single_type ())
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	674 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	675 FloatMatrix L = argl.float_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	676 FloatMatrix U = argu.float_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	677 FloatMatrix x = argx.float_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	678 FloatMatrix y = argy.float_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	679
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	680 FloatLU fact (L, U, P);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	681 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	682 fact.update_piv (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	683 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	684 fact.update (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	685
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	686 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	687 retval(2) = fact.P ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	688 retval(1) = get_lu_u (fact);
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	689 retval(0) = get_lu_l (fact);
9708 6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	690 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	691 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	692 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	693 Matrix L = argl.matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	694 Matrix U = argu.matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	695 Matrix x = argx.matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	696 Matrix y = argy.matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	697
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	698 LU fact (L, U, P);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	699 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	700 fact.update_piv (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	701 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	702 fact.update (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	703
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	704 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	705 retval(2) = fact.P ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	706 retval(1) = get_lu_u (fact);
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	707 retval(0) = get_lu_l (fact);
9708 6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	708 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	709 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	710 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	711 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	712 // complex case
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	713 if (argl.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	714 \|\| argu.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	715 \|\| argx.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	716 \|\| argy.is_single_type ())
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	717 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	718 FloatComplexMatrix L = argl.float_complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	719 FloatComplexMatrix U = argu.float_complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	720 FloatComplexMatrix x = argx.float_complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	721 FloatComplexMatrix y = argy.float_complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	722
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	723 FloatComplexLU fact (L, U, P);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	724 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	725 fact.update_piv (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	726 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	727 fact.update (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	728
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	729 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	730 retval(2) = fact.P ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	731 retval(1) = get_lu_u (fact);
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	732 retval(0) = get_lu_l (fact);
9708 6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	733 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	734 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	735 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	736 ComplexMatrix L = argl.complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	737 ComplexMatrix U = argu.complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	738 ComplexMatrix x = argx.complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	739 ComplexMatrix y = argy.complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	740
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	741 ComplexLU fact (L, U, P);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	742 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	743 fact.update_piv (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	744 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	745 fact.update (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	746
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	747 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	748 retval(2) = fact.P ();
9715 9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	749 retval(1) = get_lu_u (fact);
9f27172fbd1e auto-set MatrixType from certain functions Jaroslav Hajek <highegg@gmail.com> parents: 9709 diff changeset	750 retval(0) = get_lu_l (fact);
9708 6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	751 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	752 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	753 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	754 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	755 error ("luupdate: dimensions mismatch");
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	756 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	757 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	758 error ("luupdate: expecting numeric arguments");
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	759
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	760 return retval;
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	761 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	762
7814 87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	763 /*
10080 cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	764 %!shared A, u, v, Ac, uc, vc
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	765 %! A = [0.091364 0.613038 0.999083;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	766 %! 0.594638 0.425302 0.603537;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	767 %! 0.383594 0.291238 0.085574;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	768 %! 0.265712 0.268003 0.238409;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	769 %! 0.669966 0.743851 0.445057 ];
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	770 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	771 %! u = [0.85082;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	772 %! 0.76426;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	773 %! 0.42883;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	774 %! 0.53010;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	775 %! 0.80683 ];
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	776 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	777 %! v = [0.98810;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	778 %! 0.24295;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	779 %! 0.43167 ];
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	780 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	781 %! Ac = [0.620405 + 0.956953i 0.480013 + 0.048806i 0.402627 + 0.338171i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	782 %! 0.589077 + 0.658457i 0.013205 + 0.279323i 0.229284 + 0.721929i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	783 %! 0.092758 + 0.345687i 0.928679 + 0.241052i 0.764536 + 0.832406i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	784 %! 0.912098 + 0.721024i 0.049018 + 0.269452i 0.730029 + 0.796517i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	785 %! 0.112849 + 0.603871i 0.486352 + 0.142337i 0.355646 + 0.151496i ];
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	786 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	787 %! uc = [0.20351 + 0.05401i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	788 %! 0.13141 + 0.43708i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	789 %! 0.29808 + 0.08789i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	790 %! 0.69821 + 0.38844i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	791 %! 0.74871 + 0.25821i ];
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	792 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	793 %! vc = [0.85839 + 0.29468i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	794 %! 0.20820 + 0.93090i;
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	795 %! 0.86184 + 0.34689i ];
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	796 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	797
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	798 %!test
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	799 %! [L,U,P] = lu(A);
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	800 %! [L,U] = luupdate(L,U,P*u,v);
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	801 %! assert(norm(vec(tril(L)-L),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	802 %! assert(norm(vec(triu(U)-U),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	803 %! assert(norm(vec(P'LU - A - uv.'),Inf) < norm(A)1e1*eps)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	804 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	805 %!test
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	806 %! [L,U,P] = lu(Ac);
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	807 %! [L,U] = luupdate(L,U,P*uc,vc);
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	808 %! assert(norm(vec(tril(L)-L),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	809 %! assert(norm(vec(triu(U)-U),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	810 %! assert(norm(vec(P'LU - Ac - ucvc.'),Inf) < norm(Ac)1e1*eps)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	811
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	812 %!test
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	813 %! [L,U,P] = lu(single(A));
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	814 %! [L,U] = luupdate(L,U,P*single(u),single(v));
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	815 %! assert(norm(vec(tril(L)-L),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	816 %! assert(norm(vec(triu(U)-U),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	817 %! assert(norm(vec(P'LU - single(A) - single(u)single(v).'),Inf) < norm(single(A))1e1*eps('single'))
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	818 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	819 %!test
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	820 %! [L,U,P] = lu(single(Ac));
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	821 %! [L,U] = luupdate(L,U,P*single(uc),single(vc));
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	822 %! assert(norm(vec(tril(L)-L),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	823 %! assert(norm(vec(triu(U)-U),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	824 %! assert(norm(vec(P'LU - single(Ac) - single(uc)single(vc).'),Inf) < norm(single(Ac))1e1*eps('single'))
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	825
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	826 %!test
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	827 %! [L,U,P] = lu(A);
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	828 %! [L,U,P] = luupdate(L,U,P,u,v);
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	829 %! assert(norm(vec(tril(L)-L),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	830 %! assert(norm(vec(triu(U)-U),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	831 %! assert(norm(vec(P'LU - A - uv.'),Inf) < norm(A)1e1*eps)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	832 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	833 %!test
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	834 %! [L,U,P] = lu(Ac);
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	835 %! [L,U,P] = luupdate(L,U,P,uc,vc);
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	836 %! assert(norm(vec(tril(L)-L),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	837 %! assert(norm(vec(triu(U)-U),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	838 %! assert(norm(vec(P'LU - Ac - ucvc.'),Inf) < norm(Ac)1e1*eps)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	839
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	840 %!test
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	841 %! [L,U,P] = lu(single(A));
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	842 %! [L,U,P] = luupdate(L,U,P,single(u),single(v));
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	843 %! assert(norm(vec(tril(L)-L),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	844 %! assert(norm(vec(triu(U)-U),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	845 %! assert(norm(vec(P'LU - single(A) - single(u)single(v).'),Inf) < norm(single(A))1e1*eps('single'))
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	846 %!
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	847 %!test
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	848 %! [L,U,P] = lu(single(Ac));
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	849 %! [L,U,P] = luupdate(L,U,P,single(uc),single(vc));
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	850 %! assert(norm(vec(tril(L)-L),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	851 %! assert(norm(vec(triu(U)-U),Inf) == 0)
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	852 %! assert(norm(vec(P'LU - single(Ac) - single(uc)single(vc).'),Inf) < norm(single(Ac))1e1*eps('single'))
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	853 */
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	854
cf70ee43077c add tests for LU updating Jaroslav Hajek <highegg@gmail.com> parents: 9724 diff changeset	855 /*
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	856 ;;; Local Variables: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	857 ;;; mode: C++ ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	858 ;;; End: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	859 */

Mercurial > hg > octave-lyh

annotate src/DLD-FUNCTIONS/lu.cc @ 10080:cf70ee43077c