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

annotate src/DLD-FUNCTIONS/lu.cc @ 9708:6f3ffe11d926

implement luupdate

author	Jaroslav Hajek <highegg@gmail.com>
date	Thu, 08 Oct 2009 16:05:53 +0200
parents	923c7cb7f13f
children	f8e2e9fdaa8f

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
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	45 maybe_set_triangular (const MT& m, MatrixType::matrix_type t = MatrixType::Upper)
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	46 {
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	47 typedef typename MT::element_type T;
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	48 octave_value retval;
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	49 octave_idx_type r = m.rows (), c = m.columns ();
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	50 if (r == c)
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	51 {
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	52 const T zero = T();
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	53 octave_idx_type i = 0;
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	54 for (;i != r && m(i,i) != zero; i++) ;
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	55 if (i == r)
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	56 retval = octave_value (m, MatrixType (t));
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	57 else
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	58 retval = m;
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	59 }
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	60 else
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	61 retval = m;
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	62
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	63 return retval;
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	64 }
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	65
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	66 DEFUN_DLD (lu, args, nargout,
3548 ab7fa5a8f23f [project @ 2000-02-03 01:17:15 by jwe] jwe parents: 3372 diff changeset	67 "-- texinfo --\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	68 @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	69 @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	70 @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	71 @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	72 @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	73 @deftypefnx {Loadable Function} {[@dots{}] =} lu (@dots{}, 'vector')\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	74 @cindex LU decomposition\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	75 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	76 @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	77 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	78 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	79 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	80 @example\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	81 [l, u, p] = lu (a)\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	82 @end example\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	83 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	84 @noindent\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	85 returns\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	86 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	87 @example\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	88 @group\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	89 l =\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	90 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	91 1.00000 0.00000\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	92 0.33333 1.00000\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	93 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	94 u =\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	95 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	96 3.00000 4.00000\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	97 0.00000 0.66667\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	98 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	99 p =\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	100 \n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	101 0 1\n\
f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	102 1 0\n\
9064 7c02ec148a3c Check grammar on all .cc files Rik <rdrider0-list@yahoo.com> parents: 8969 diff changeset	103 @end group\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	104 @end example\n\
4329 d53c33d93440 [project @ 2003-02-18 20:00:48 by jwe] jwe parents: 3587 diff changeset	105 \n\
5457 c6dc1ccd83a9 [project @ 2005-09-19 19:23:35 by jwe] jwe parents: 5307 diff changeset	106 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	107 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	108 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	109 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	110 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	111 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	112 @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	113 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	114 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	115 @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	116 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	117 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	118 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	119 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	120 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	121 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	122 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	123 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	124 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	125 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	126 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	127 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	128 @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	129 (@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	130 (:, @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	131 \n\
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	132 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	133 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	134 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	135 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	136 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	137 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	138 information.\n\
3372 f16c2ce14886 [project @ 1999-11-23 19:07:09 by jwe] jwe parents: 3014 diff changeset	139 @end deftypefn")
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	140 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	141 octave_value_list retval;
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	142 int nargin = args.length ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	143 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	144 bool scale = (nargout == 5);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	145
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	146 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	147 \|\| (!issparse && (nargin > 2 \|\| nargout > 3)))
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	148 {
5823 080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe] jwe parents: 5457 diff changeset	149 print_usage ();
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	150 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	151 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	152
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	153 bool vecout = false;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	154 Matrix thres;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	155
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	156 int n = 1;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	157 while (n < nargin && ! error_state)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	158 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	159 if (args (n).is_string ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	160 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	161 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	162
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	163 if (! error_state )
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	164 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	165 if (tmp.compare ("vector") == 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	166 vecout = true;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	167 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	168 error ("lu: unrecognized string argument");
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 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	171 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	172 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	173 Matrix tmp = args(n++).matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	174
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	175 if (! error_state )
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	176 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	177 if (!issparse)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	178 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	179 else if (tmp.nelem () == 1)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	180 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	181 thres.resize(1,2);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	182 thres(0) = tmp(0);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	183 thres(1) = tmp(0);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	184 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	185 else if (tmp.nelem () == 2)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	186 thres = tmp;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	187 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	188 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	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 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	192
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	193 octave_value arg = args(0);
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	194
5275 23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4329 diff changeset	195 octave_idx_type nr = arg.rows ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe] jwe parents: 4329 diff changeset	196 octave_idx_type nc = arg.columns ();
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	197
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	198 int arg_is_empty = empty_arg ("lu", nr, nc);
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	199
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	200 if (issparse)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	201 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	202 if (arg_is_empty < 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	203 return retval;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	204 else if (arg_is_empty > 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	205 return octave_value_list (5, SparseMatrix ());
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	206
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	207 ColumnVector Qinit;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	208 if (nargout < 4)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	209 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	210 Qinit.resize (nc);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	211 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	212 Qinit (i) = i;
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
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	215 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	216 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	217 SparseMatrix m = arg.sparse_matrix_value ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	218
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	219 switch (nargout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	220 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	221 case 0:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	222 case 1:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	223 case 2:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	224 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	225 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	226
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	227 if (nargout < 2)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	228 retval (0) = fact.Y ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	229 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	230 {
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	231 PermMatrix P = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	232 SparseMatrix L = P.transpose () * fact.L ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	233 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	234 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	235
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	236 retval(0) = octave_value (L,
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	237 MatrixType (MatrixType::Permuted_Lower,
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	238 nr, fact.row_perm ()));
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 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	241 break;
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	242
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	243 case 3:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	244 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	245 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	246
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	247 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	248 retval (2) = fact.Pr_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	249 else
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	250 retval(2) = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	251
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	252 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	253 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	254 retval(0) = octave_value (fact.L (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	255 MatrixType (MatrixType::Lower));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	256 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	257 break;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	258
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	259 case 4:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	260 default:
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	261 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	262 SparseLU fact (m, thres, scale);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	263
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	264 if (scale)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	265 retval(4) = fact.R ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	266
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	267 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	268 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	269 retval(3) = fact.Pc_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	270 retval(2) = fact.Pr_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	271 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	272 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	273 {
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	274 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	275 retval(2) = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	276 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	277 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	278 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	279 retval(0) = octave_value (fact.L (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	280 MatrixType (MatrixType::Lower));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	281 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	282 break;
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 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	285 else if (arg.is_complex_type ())
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	286 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	287 SparseComplexMatrix m = arg.sparse_complex_matrix_value ();
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	288
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	289 switch (nargout)
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	290 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	291 case 0:
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	292 case 1:
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	293 case 2:
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	294 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	295 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	296
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	297 if (nargout < 2)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	298 retval (0) = fact.Y ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	299 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	300 {
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	301 PermMatrix P = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	302 SparseComplexMatrix L = P.transpose () * fact.L ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	303 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	304 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	305
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	306 retval(0) = octave_value (L,
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	307 MatrixType (MatrixType::Permuted_Lower,
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	308 nr, fact.row_perm ()));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	309 }
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	310 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	311 break;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	312
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	313 case 3:
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	314 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	315 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	316
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	317 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	318 retval (2) = fact.Pr_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	319 else
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	320 retval(2) = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	321
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	322 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	323 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	324 retval(0) = octave_value (fact.L (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	325 MatrixType (MatrixType::Lower));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	326 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	327 break;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	328
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	329 case 4:
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	330 default:
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	331 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	332 SparseComplexLU fact (m, thres, scale);
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	333
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	334 if (scale)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	335 retval(4) = fact.R ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	336
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	337 if (vecout)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	338 {
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	339 retval(3) = fact.Pc_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	340 retval(2) = fact.Pr_vec ();
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	341 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	342 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	343 {
8969 3ecbc236e2e0 Have sparse LU return permutation matrices rather than sparse matrices. Jason Riedy <jason@acm.org> parents: 8920 diff changeset	344 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	345 retval(2) = fact.Pr_mat ();
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	346 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	347 retval(1) = octave_value (fact.U (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	348 MatrixType (MatrixType::Upper));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	349 retval(0) = octave_value (fact.L (),
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	350 MatrixType (MatrixType::Lower));
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	351 }
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	352 break;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	353 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	354 }
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	355 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	356 gripe_wrong_type_arg ("lu", arg);
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	357 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	358 else
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	359 {
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	360 if (arg_is_empty < 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	361 return retval;
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	362 else if (arg_is_empty > 0)
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	363 return octave_value_list (3, Matrix ());
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	364
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	365 if (arg.is_real_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	366 {
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	367 if (arg.is_single_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	368 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	369 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	370
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	371 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	372 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	373 FloatLU fact (m);
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	374
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	375 switch (nargout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	376 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	377 case 0:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	378 case 1:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	379 retval(0) = fact.Y ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	380 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	381
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	382 case 2:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	383 {
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	384 PermMatrix P = fact.P ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	385 FloatMatrix L = P.transpose () * fact.L ();
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	386 retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	387 retval(0) = L;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	388 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	389 break;
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	390
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	391 case 3:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	392 default:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	393 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	394 if (vecout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	395 retval(2) = fact.P_vec ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	396 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	397 retval(2) = fact.P ();
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	398 retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper);
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	399 retval(0) = maybe_set_triangular (fact.L (), MatrixType::Lower);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	400 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	401 break;
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 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	405 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	406 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	407 Matrix m = arg.matrix_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	408
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	409 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	410 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	411 LU fact (m);
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	412
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	413 switch (nargout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	414 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	415 case 0:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	416 case 1:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	417 retval(0) = fact.Y ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	418 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	419
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	420 case 2:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	421 {
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	422 PermMatrix P = fact.P ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	423 Matrix L = P.transpose () * fact.L ();
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	424 retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	425 retval(0) = L;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	426 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	427 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	428
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	429 case 3:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	430 default:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	431 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	432 if (vecout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	433 retval(2) = fact.P_vec ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	434 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	435 retval(2) = fact.P ();
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	436 retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper);
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	437 retval(0) = maybe_set_triangular (fact.L (), MatrixType::Lower);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	438 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	439 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	440 }
7515 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 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	444 else if (arg.is_complex_type ())
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	445 {
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	446 if (arg.is_single_type ())
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	447 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	448 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	449
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	450 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	451 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	452 FloatComplexLU fact (m);
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	453
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	454 switch (nargout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	455 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	456 case 0:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	457 case 1:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	458 retval(0) = fact.Y ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	459 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	460
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	461 case 2:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	462 {
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	463 PermMatrix P = fact.P ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	464 FloatComplexMatrix L = P.transpose () * fact.L ();
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	465 retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	466 retval(0) = L;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	467 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	468 break;
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	469
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	470 case 3:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	471 default:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	472 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	473 if (vecout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	474 retval(2) = fact.P_vec ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	475 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	476 retval(2) = fact.P ();
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	477 retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper);
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	478 retval(0) = maybe_set_triangular (fact.L (), MatrixType::Lower);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	479 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	480 break;
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 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	484 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	485 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	486 ComplexMatrix m = arg.complex_matrix_value ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	487
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	488 if (! error_state)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	489 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	490 ComplexLU fact (m);
7515 f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	491
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	492 switch (nargout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	493 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	494 case 0:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	495 case 1:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	496 retval(0) = fact.Y ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	497 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	498
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	499 case 2:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	500 {
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	501 PermMatrix P = fact.P ();
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	502 ComplexMatrix L = P.transpose () * fact.L ();
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	503 retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	504 retval(0) = L;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	505 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	506 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	507
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	508 case 3:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	509 default:
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	510 {
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	511 if (vecout)
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	512 retval(2) = fact.P_vec ();
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	513 else
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	514 retval(2) = fact.P ();
8869 c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	515 retval(1) = maybe_set_triangular (fact.U (), MatrixType::Upper);
c3b743b1b1c6 preset triangular type if possible for lu and qr outputs Jaroslav Hajek <highegg@gmail.com> parents: 8367 diff changeset	516 retval(0) = maybe_set_triangular (fact.L (), MatrixType::Lower);
7789 82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	517 }
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	518 break;
82be108cc558 First attempt at single precision tyeps David Bateman <dbateman@free.fr> parents: 7515 diff changeset	519 }
7515 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 }
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	523 else
f3c00dc0912b Eliminate the rest of the dispatched sparse functions David Bateman <dbateman@free.fr> parents: 7017 diff changeset	524 gripe_wrong_type_arg ("lu", arg);
2928 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
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	527 return retval;
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
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	530 /*
7814 87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	531
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	532 %!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	533
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	534 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	535 %! [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	536 %! 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	537 %! 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	538
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	539 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	540 %! [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	541 %! 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	542 %! assert(u, [3, 4; 0, 2/3], sqrt (eps));
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	543 %! 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	544
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	545 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	546 %! [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	547 %! 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	548 %! 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	549 %! 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	550
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	551 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	552 %! [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	553 %! 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	554 %! assert(u, [5, 6; 0, 4/5], sqrt (eps));
8367 445d27d79f4e support permutation matrix objects Jaroslav Hajek <highegg@gmail.com> parents: 7814 diff changeset	555 %! 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	556
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	557 %!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	558
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	559 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	560 %! [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	561 %! 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	562 %! 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	563
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	564 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	565 %! [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	566 %! 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	567 %! 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	568 %! 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	569
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	570 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	571 %! [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	572 %! 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	573 %! 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	574 %! 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	575
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	576 %!test
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	577 %! [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	578 %! 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	579 %! 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	580 %! 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	581
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	582 %!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	583 %!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	584
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	585 */
87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	586
9708 6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	587 static
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	588 bool check_lu_dims (const octave_value& l, const octave_value& u,
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	589 const octave_value& p)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	590 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	591 octave_idx_type m = l.rows (), k = u.rows (), n = u.columns ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	592 return ((l.ndims () == 2 && u.ndims () == 2 && k == l.columns ())
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	593 && k == std::min (m, n) &&
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	594 (p.is_undefined () \|\| p.rows () == m));
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	595 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	596
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	597 DEFUN_DLD (luupdate, args, nargout,
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	598 "-- texinfo --\n\
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	599 @deftypefn {Loadable Function} {[@var{l}, @var{u}] =} luupdate (@var{l}, @var{u}, @var{x}, @var{y})\n\
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	600 @deftypefn {Loadable Function} {[@var{l}, @var{u}, @var{p}] =}\
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	601 luupdate (@var{l}, @var{u}, @var{p}, @var{x}, @var{y})\n\
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	602 ")
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	603 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	604 octave_idx_type nargin = args.length ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	605 octave_value_list retval;
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	606
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	607 bool pivoted = nargin == 5;
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	608
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	609 if (nargin != 4 && nargin != 5)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	610 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	611 print_usage ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	612 return retval;
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	613 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	614
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	615 octave_value argl = args(0);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	616 octave_value argu = args(1);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	617 octave_value argp = pivoted ? args(2) : octave_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	618 octave_value argx = args(2 + pivoted);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	619 octave_value argy = args(3 + pivoted);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	620
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	621 if (argl.is_numeric_type () && argu.is_numeric_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	622 && argx.is_numeric_type () && argy.is_numeric_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	623 && (! pivoted \|\| argp.is_perm_matrix ()))
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	624 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	625 if (check_lu_dims (argl, argu, argp))
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	626 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	627 PermMatrix P = (pivoted
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	628 ? argp.perm_matrix_value ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	629 : PermMatrix::eye (argl.rows ()));
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	630
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	631 if (argl.is_real_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	632 && argu.is_real_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	633 && argx.is_real_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	634 && argy.is_real_type ())
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	635 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	636 // all real case
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	637 if (argl.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	638 \|\| argu.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	639 \|\| argx.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	640 \|\| argy.is_single_type ())
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 FloatMatrix L = argl.float_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	643 FloatMatrix U = argu.float_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	644 FloatMatrix x = argx.float_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	645 FloatMatrix y = argy.float_matrix_value ();
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 FloatLU fact (L, U, P);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	648 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	649 fact.update_piv (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	650 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	651 fact.update (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	652
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	653 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	654 retval(2) = fact.P ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	655 retval(1) = fact.U ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	656 retval(0) = fact.L ();
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 else
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 Matrix L = argl.matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	661 Matrix U = argu.matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	662 Matrix x = argx.matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	663 Matrix y = argy.matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	664
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	665 LU fact (L, U, P);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	666 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	667 fact.update_piv (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	668 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	669 fact.update (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	670
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	671 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	672 retval(2) = fact.P ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	673 retval(1) = fact.U ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	674 retval(0) = fact.L ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	675 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	676 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	677 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	678 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	679 // complex case
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	680 if (argl.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	681 \|\| argu.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	682 \|\| argx.is_single_type ()
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	683 \|\| argy.is_single_type ())
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	684 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	685 FloatComplexMatrix L = argl.float_complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	686 FloatComplexMatrix U = argu.float_complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	687 FloatComplexMatrix x = argx.float_complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	688 FloatComplexMatrix y = argy.float_complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	689
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	690 FloatComplexLU fact (L, U, P);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	691 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	692 fact.update_piv (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	693 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	694 fact.update (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	695
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	696 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	697 retval(2) = fact.P ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	698 retval(1) = fact.U ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	699 retval(0) = fact.L ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	700 }
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 {
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	703 ComplexMatrix L = argl.complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	704 ComplexMatrix U = argu.complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	705 ComplexMatrix x = argx.complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	706 ComplexMatrix y = argy.complex_matrix_value ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	707
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	708 ComplexLU fact (L, U, P);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	709 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	710 fact.update_piv (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	711 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	712 fact.update (x, y);
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	713
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	714 if (pivoted)
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	715 retval(2) = fact.P ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	716 retval(1) = fact.U ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	717 retval(0) = fact.L ();
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	718 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	719 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	720 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	721 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	722 error ("luupdate: dimensions mismatch");
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	723 }
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	724 else
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	725 error ("luupdate: expecting numeric arguments");
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	726
6f3ffe11d926 implement luupdate Jaroslav Hajek <highegg@gmail.com> parents: 9209 diff changeset	727 return retval;
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
7814 87865ed7405f Second set of single precision test code and fix of resulting bugs David Bateman <dbateman@free.fr> parents: 7789 diff changeset	730 /*
2928 295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	731 ;;; Local Variables: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	732 ;;; mode: C++ ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	733 ;;; End: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe] jwe parents: diff changeset	734 */

Mercurial > octave-nkf

annotate src/DLD-FUNCTIONS/lu.cc @ 9708:6f3ffe11d926