Mercurial > octave

changeset 21576:82092a17fa50

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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.
author Rik <rik@octave.org>
date Fri, 01 Apr 2016 14:06:45 -0700
parents bc9aa534bc29
children 31823239207e
files libinterp/corefcn/__luinc__.cc libinterp/corefcn/luinc.cc libinterp/corefcn/module.mk
diffstat 3 files changed, 286 insertions(+), 335 deletions(-) [+]
line wrap: on
line diff
--- /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
+<http://www.gnu.org/licenses/>.
+
+*/
+
+#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<SparseMatrix> 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<SparseMatrix> 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<SparseMatrix> 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<SparseComplexMatrix> 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<SparseComplexMatrix> 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<SparseComplexMatrix> 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;
+}
+
--- 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
-<http://www.gnu.org/licenses/>.
-
-*/
-
-#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<SparseMatrix> 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<SparseMatrix> 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<SparseMatrix> 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<SparseComplexMatrix> 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<SparseComplexMatrix> 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<SparseComplexMatrix> 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;
-}
-
--- 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 \