# HG changeset patch # User Rik # Date 1459544805 25200 # Node ID 82092a17fa50cc3ed46e9801616ee0a46dd177b7 # Parent bc9aa534bc29908ef768138c1a9e3ec9accc4965 maint: Rename luinc.cc to __luinc__.cc to reflect file contents. * __luinc__.cc: Renamed from luinc.cc * libinterp/corefcn/module.mk: Update build system to use new name. diff -r bc9aa534bc29 -r 82092a17fa50 libinterp/corefcn/__luinc__.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/__luinc__.cc Fri Apr 01 14:06:45 2016 -0700 @@ -0,0 +1,285 @@ +/* + +Copyright (C) 2005-2015 David Bateman + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +. + +*/ + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + +#include "defun.h" +#include "error.h" +#include "errwarn.h" +#include "ovl.h" +#include "utils.h" +#include "oct-map.h" + +#include "MatrixType.h" +#include "sparse-lu.h" +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" + +// FIXME: Deprecated in 4.0 and should be removed in 4.4. +DEFUN (__luinc__, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} __luinc__ (@var{A}, '0')\n\ +@deftypefnx {} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} __luinc__ (@var{A}, @var{droptol})\n\ +@deftypefnx {} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} __luinc__ (@var{A}, @var{opts})\n\ +Internal implementation of @code{luinc}.\n\ +\n\ +See documentation for @code{luinc}.\n\ +@seealso{luinc}\n\ +@end deftypefn") +{ + int nargin = args.length (); + + if (nargin < 2 || nargin > 3) + print_usage (); + + if (! args(0).is_sparse_type ()) + error ("luinc: matrix A must be sparse"); + + bool zero_level = false; + bool milu = false; + bool udiag = false; + Matrix thresh; + double droptol = -1.0; + bool vecout = false; + + if (args(1).is_string ()) + { + if (args(1).string_value () == "0") + zero_level = true; + else + error ("luinc: unrecognized string argument"); + } + else if (args(1).is_map ()) + { + octave_scalar_map map = args(1).xscalar_map_value ("luinc: OPTS must be a scalar structure"); + + octave_value tmp; + + tmp = map.getfield ("droptol"); + if (tmp.is_defined ()) + droptol = tmp.double_value (); + + tmp = map.getfield ("milu"); + if (tmp.is_defined ()) + { + double val = tmp.double_value (); + + milu = (val == 0.0 ? false : true); + } + + tmp = map.getfield ("udiag"); + if (tmp.is_defined ()) + { + double val = tmp.double_value (); + + udiag = (val == 0.0 ? false : true); + } + + tmp = map.getfield ("thresh"); + if (tmp.is_defined ()) + { + thresh = tmp.matrix_value (); + + if (thresh.numel () == 1) + { + thresh.resize (1, 2); + thresh(1) = thresh(0); + } + else if (thresh.numel () != 2) + error ("luinc: THRESH must be a 1 or 2-element vector"); + } + } + else + droptol = args(1).double_value (); + + if (nargin == 3) + { + std::string tmp = args(2).string_value (); + + if (tmp == "vector") + vecout = true; + else + error ("luinc: unrecognized string argument"); + } + + // FIXME: Add code for zero-level factorization + if (zero_level) + error ("luinc: zero-level factorization not implemented"); + + octave_value_list retval; + + if (args(0).is_real_type ()) + { + SparseMatrix sm = args(0).sparse_matrix_value (); + octave_idx_type sm_nr = sm.rows (); + octave_idx_type sm_nc = sm.cols (); + ColumnVector Qinit (sm_nc); + + for (octave_idx_type i = 0; i < sm_nc; i++) + Qinit(i) = i; + + switch (nargout) + { + case 0: + case 1: + case 2: + { + sparse_lu fact (sm, Qinit, thresh, false, true, droptol, + milu, udiag); + + SparseMatrix P = fact.Pr (); + SparseMatrix L = P.transpose () * fact.L (); + + retval(1) + = octave_value (fact.U (), MatrixType (MatrixType::Upper)); + + retval(0) + = octave_value (L, MatrixType (MatrixType::Permuted_Lower, + sm_nr, fact.row_perm ())); + } + break; + + case 3: + { + sparse_lu fact (sm, Qinit, thresh, false, true, droptol, + milu, udiag); + + if (vecout) + retval(2) = fact.Pr_vec (); + else + retval(2) = fact.Pr_mat (); + + retval(1) + = octave_value (fact.U (), MatrixType (MatrixType::Upper)); + + retval(0) + = octave_value (fact.L (), MatrixType (MatrixType::Lower)); + } + break; + + case 4: + default: + { + sparse_lu fact (sm, Qinit, thresh, false, false, droptol, + milu, udiag); + + if (vecout) + { + retval(3) = fact.Pc_vec (); + retval(2) = fact.Pr_vec (); + } + else + { + retval(3) = fact.Pc_mat (); + retval(2) = fact.Pr_mat (); + } + + retval(1) + = octave_value (fact.U (), MatrixType (MatrixType::Upper)); + + retval(0) + = octave_value (fact.L (), MatrixType (MatrixType::Lower)); + } + break; + } + } + else + { + SparseComplexMatrix sm = args(0).sparse_complex_matrix_value (); + octave_idx_type sm_nr = sm.rows (); + octave_idx_type sm_nc = sm.cols (); + ColumnVector Qinit (sm_nc); + + for (octave_idx_type i = 0; i < sm_nc; i++) + Qinit(i) = i; + + switch (nargout) + { + case 0: + case 1: + case 2: + { + sparse_lu fact (sm, Qinit, thresh, false, true, + droptol, milu, udiag); + + SparseMatrix P = fact.Pr (); + SparseComplexMatrix L = P.transpose () * fact.L (); + + retval(1) + = octave_value (fact.U (), MatrixType (MatrixType::Upper)); + + retval(0) + = octave_value (L, MatrixType (MatrixType::Permuted_Lower, + sm_nr, fact.row_perm ())); + } + break; + + case 3: + { + sparse_lu fact (sm, Qinit, thresh, false, true, + droptol, milu, udiag); + + if (vecout) + retval(2) = fact.Pr_vec (); + else + retval(2) = fact.Pr_mat (); + + retval(1) + = octave_value (fact.U (), MatrixType (MatrixType::Upper)); + + retval(0) + = octave_value (fact.L (), MatrixType (MatrixType::Lower)); + } + break; + + case 4: + default: + { + sparse_lu fact (sm, Qinit, thresh, false, false, + droptol, milu, udiag); + + if (vecout) + { + retval(3) = fact.Pc_vec (); + retval(2) = fact.Pr_vec (); + } + else + { + retval(3) = fact.Pc_mat (); + retval(2) = fact.Pr_mat (); + } + + retval(1) + = octave_value (fact.U (), MatrixType (MatrixType::Upper)); + + retval(0) + = octave_value (fact.L (), MatrixType (MatrixType::Lower)); + } + break; + } + } + + return retval; +} + diff -r bc9aa534bc29 -r 82092a17fa50 libinterp/corefcn/luinc.cc --- a/libinterp/corefcn/luinc.cc Fri Apr 01 13:56:09 2016 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,334 +0,0 @@ -/* - -Copyright (C) 2005-2015 David Bateman - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -. - -*/ - -#ifdef HAVE_CONFIG_H -# include "config.h" -#endif - -#include "defun.h" -#include "error.h" -#include "errwarn.h" -#include "ovl.h" -#include "utils.h" -#include "oct-map.h" - -#include "MatrixType.h" -#include "sparse-lu.h" -#include "ov-re-sparse.h" -#include "ov-cx-sparse.h" - -// FIXME: Deprecated in 4.0 and should be removed in 4.4. -DEFUN (__luinc__, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} luinc (@var{A}, '0')\n\ -@deftypefnx {} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} luinc (@var{A}, @var{droptol})\n\ -@deftypefnx {} {[@var{L}, @var{U}, @var{P}, @var{Q}] =} luinc (@var{A}, @var{opts})\n\ -@cindex LU decomposition\n\ -Produce the incomplete LU@tie{}factorization of the sparse matrix @var{A}.\n\ -\n\ -Two types of incomplete factorization are possible, and the type\n\ -is determined by the second argument to @code{luinc}.\n\ -\n\ -Called with a second argument of @qcode{'0'}, the zero-level incomplete\n\ -LU@tie{}factorization is produced. This creates a factorization of @var{A}\n\ -where the position of the nonzero arguments correspond to the same\n\ -positions as in the matrix @var{A}.\n\ -\n\ -Alternatively, the fill-in of the incomplete LU@tie{}factorization can\n\ -be controlled through the variable @var{droptol} or the structure\n\ -@var{opts}. The @sc{umfpack} multifrontal factorization code by Tim A.\n\ -Davis is used for the incomplete LU@tie{}factorization, (availability\n\ -@url{http://www.cise.ufl.edu/research/sparse/umfpack/})\n\ -\n\ -@var{droptol} determines the values below which the values in the\n\ -LU@tie{} factorization are dropped and replaced by zero. It must be a\n\ -positive scalar, and any values in the factorization whose absolute value\n\ -are less than this value are dropped, expect if leaving them increase the\n\ -sparsity of the matrix. Setting @var{droptol} to zero results in a complete\n\ -LU@tie{}factorization which is the default.\n\ -\n\ -@var{opts} is a structure containing one or more of the fields\n\ -\n\ -@table @code\n\ -@item droptol\n\ -The drop tolerance as above. If @var{opts} only contains @code{droptol}\n\ -then this is equivalent to using the variable @var{droptol}.\n\ -\n\ -@item milu\n\ -A logical variable flagging whether to use the modified incomplete\n\ -LU@tie{} factorization. In the case that @code{milu} is true, the dropped\n\ -values are subtracted from the diagonal of the matrix @var{U} of the\n\ -factorization. The default is @code{false}.\n\ -\n\ -@item udiag\n\ -A logical variable that flags whether zero elements on the diagonal of\n\ -@var{U} should be replaced with @var{droptol} to attempt to avoid singular\n\ -factors. The default is @code{false}.\n\ -\n\ -@item thresh\n\ -Defines the pivot threshold in the interval [0,1]. Values outside that\n\ -range are ignored.\n\ -@end table\n\ -\n\ -All other fields in @var{opts} are ignored. The outputs from @code{luinc}\n\ -are the same as for @code{lu}.\n\ -\n\ -Given the string argument @qcode{\"vector\"}, @code{luinc} returns the\n\ -values of @var{p} @var{q} as vector values.\n\ -@seealso{sparse, lu, ilu, ichol}\n\ -@end deftypefn") -{ - int nargin = args.length (); - - if (nargin < 2 || nargin > 3) - print_usage (); - - if (! args(0).is_sparse_type ()) - error ("luinc: matrix A must be sparse"); - - bool zero_level = false; - bool milu = false; - bool udiag = false; - Matrix thresh; - double droptol = -1.0; - bool vecout = false; - - if (args(1).is_string ()) - { - if (args(1).string_value () == "0") - zero_level = true; - else - error ("luinc: unrecognized string argument"); - } - else if (args(1).is_map ()) - { - octave_scalar_map map = args(1).xscalar_map_value ("luinc: OPTS must be a scalar structure"); - - octave_value tmp; - - tmp = map.getfield ("droptol"); - if (tmp.is_defined ()) - droptol = tmp.double_value (); - - tmp = map.getfield ("milu"); - if (tmp.is_defined ()) - { - double val = tmp.double_value (); - - milu = (val == 0.0 ? false : true); - } - - tmp = map.getfield ("udiag"); - if (tmp.is_defined ()) - { - double val = tmp.double_value (); - - udiag = (val == 0.0 ? false : true); - } - - tmp = map.getfield ("thresh"); - if (tmp.is_defined ()) - { - thresh = tmp.matrix_value (); - - if (thresh.numel () == 1) - { - thresh.resize (1, 2); - thresh(1) = thresh(0); - } - else if (thresh.numel () != 2) - error ("luinc: THRESH must be a 1 or 2-element vector"); - } - } - else - droptol = args(1).double_value (); - - if (nargin == 3) - { - std::string tmp = args(2).string_value (); - - if (tmp == "vector") - vecout = true; - else - error ("luinc: unrecognized string argument"); - } - - // FIXME: Add code for zero-level factorization - if (zero_level) - error ("luinc: zero-level factorization not implemented"); - - octave_value_list retval; - - if (args(0).is_real_type ()) - { - SparseMatrix sm = args(0).sparse_matrix_value (); - octave_idx_type sm_nr = sm.rows (); - octave_idx_type sm_nc = sm.cols (); - ColumnVector Qinit (sm_nc); - - for (octave_idx_type i = 0; i < sm_nc; i++) - Qinit(i) = i; - - switch (nargout) - { - case 0: - case 1: - case 2: - { - sparse_lu fact (sm, Qinit, thresh, false, true, droptol, - milu, udiag); - - SparseMatrix P = fact.Pr (); - SparseMatrix L = P.transpose () * fact.L (); - - retval(1) - = octave_value (fact.U (), MatrixType (MatrixType::Upper)); - - retval(0) - = octave_value (L, MatrixType (MatrixType::Permuted_Lower, - sm_nr, fact.row_perm ())); - } - break; - - case 3: - { - sparse_lu fact (sm, Qinit, thresh, false, true, droptol, - milu, udiag); - - if (vecout) - retval(2) = fact.Pr_vec (); - else - retval(2) = fact.Pr_mat (); - - retval(1) - = octave_value (fact.U (), MatrixType (MatrixType::Upper)); - - retval(0) - = octave_value (fact.L (), MatrixType (MatrixType::Lower)); - } - break; - - case 4: - default: - { - sparse_lu fact (sm, Qinit, thresh, false, false, droptol, - milu, udiag); - - if (vecout) - { - retval(3) = fact.Pc_vec (); - retval(2) = fact.Pr_vec (); - } - else - { - retval(3) = fact.Pc_mat (); - retval(2) = fact.Pr_mat (); - } - - retval(1) - = octave_value (fact.U (), MatrixType (MatrixType::Upper)); - - retval(0) - = octave_value (fact.L (), MatrixType (MatrixType::Lower)); - } - break; - } - } - else - { - SparseComplexMatrix sm = args(0).sparse_complex_matrix_value (); - octave_idx_type sm_nr = sm.rows (); - octave_idx_type sm_nc = sm.cols (); - ColumnVector Qinit (sm_nc); - - for (octave_idx_type i = 0; i < sm_nc; i++) - Qinit(i) = i; - - switch (nargout) - { - case 0: - case 1: - case 2: - { - sparse_lu fact (sm, Qinit, thresh, false, true, - droptol, milu, udiag); - - SparseMatrix P = fact.Pr (); - SparseComplexMatrix L = P.transpose () * fact.L (); - - retval(1) - = octave_value (fact.U (), MatrixType (MatrixType::Upper)); - - retval(0) - = octave_value (L, MatrixType (MatrixType::Permuted_Lower, - sm_nr, fact.row_perm ())); - } - break; - - case 3: - { - sparse_lu fact (sm, Qinit, thresh, false, true, - droptol, milu, udiag); - - if (vecout) - retval(2) = fact.Pr_vec (); - else - retval(2) = fact.Pr_mat (); - - retval(1) - = octave_value (fact.U (), MatrixType (MatrixType::Upper)); - - retval(0) - = octave_value (fact.L (), MatrixType (MatrixType::Lower)); - } - break; - - case 4: - default: - { - sparse_lu fact (sm, Qinit, thresh, false, false, - droptol, milu, udiag); - - if (vecout) - { - retval(3) = fact.Pc_vec (); - retval(2) = fact.Pr_vec (); - } - else - { - retval(3) = fact.Pc_mat (); - retval(2) = fact.Pr_mat (); - } - - retval(1) - = octave_value (fact.U (), MatrixType (MatrixType::Upper)); - - retval(0) - = octave_value (fact.L (), MatrixType (MatrixType::Lower)); - } - break; - } - } - - return retval; -} - diff -r bc9aa534bc29 -r 82092a17fa50 libinterp/corefcn/module.mk --- a/libinterp/corefcn/module.mk Fri Apr 01 13:56:09 2016 -0700 +++ b/libinterp/corefcn/module.mk Fri Apr 01 14:06:45 2016 -0700 @@ -121,6 +121,7 @@ libinterp/corefcn/__ichol__.cc \ libinterp/corefcn/__ilu__.cc \ libinterp/corefcn/__lin_interpn__.cc \ + libinterp/corefcn/__luinc__.cc \ libinterp/corefcn/__pchip_deriv__.cc \ libinterp/corefcn/__qp__.cc \ libinterp/corefcn/balance.cc \ @@ -191,7 +192,6 @@ libinterp/corefcn/ls-utils.cc \ libinterp/corefcn/lsode.cc \ libinterp/corefcn/lu.cc \ - libinterp/corefcn/luinc.cc \ libinterp/corefcn/mappers.cc \ libinterp/corefcn/matrix_type.cc \ libinterp/corefcn/max.cc \