Mercurial > octave

changeset 21585:bad3ed83330d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
symbfact.cc: Overhaul dldfcn. * symbfact.cc: Redo docstring. Eliminate nargout input validation. Improve input validation for TYP and MODE variables. Add label cleanup: and use goto statements to guarantee that memory allocated in cholmod library is cleaned up. Use std::max rather than hand-coded max routine. Add BIST tests.
author Rik <rik@octave.org>
date Mon, 04 Apr 2016 09:00:22 -0700
parents ee1a009dd60f
children 3b9524234d89
files libinterp/dldfcn/symbfact.cc
diffstat 1 files changed, 136 insertions(+), 72 deletions(-) [+]
line wrap: on
line diff
--- a/libinterp/dldfcn/symbfact.cc	Sun Apr 03 18:09:56 2016 -0700
+++ b/libinterp/dldfcn/symbfact.cc	Mon Apr 04 09:00:22 2016 -0700
@@ -41,49 +41,53 @@
 
 DEFUN_DLD (symbfact, args, nargout,
            "-*- texinfo -*-\n\
-@deftypefn  {} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{r}] =} symbfact (@var{S})\n\
+@deftypefn  {} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{R}] =} symbfact (@var{S})\n\
 @deftypefnx {} {[@dots{}] =} symbfact (@var{S}, @var{typ})\n\
 @deftypefnx {} {[@dots{}] =} symbfact (@var{S}, @var{typ}, @var{mode})\n\
 \n\
-Perform a symbolic factorization analysis on the sparse matrix @var{S}.\n\
+Perform a symbolic factorization analysis of the sparse matrix @var{S}.\n\
 \n\
 The input variables are\n\
 \n\
 @table @var\n\
 @item S\n\
-@var{S} is a complex or real sparse matrix.\n\
+@var{S} is a real or complex sparse matrix.\n\
 \n\
 @item typ\n\
 Is the type of the factorization and can be one of\n\
 \n\
-@table @samp\n\
-@item sym\n\
-Factorize @var{S}.  This is the default.\n\
+@table @asis\n\
+@item @qcode{\"sym\"} (default)\n\
+Factorize @var{S}.  Assumes @var{S} is symmetric and uses the upper\n\
+triangular portion of the matrix.\n\
 \n\
-@item col\n\
-Factorize @code{@var{S}' * @var{S}}.\n\
+@item @qcode{\"col\"}\n\
+Factorize @tcode{@var{S}' * @var{S}}.\n\
 \n\
-@item row\n\
+@item @qcode{\"row\"}\n\
 Factorize @tcode{@var{S} * @var{S}'}.\n\
 \n\
-@item lo\n\
-Factorize @tcode{@var{S}'}\n\
+@item @qcode{\"lo\"}\n\
+Factorize @tcode{@var{S}'}.  Assumes @var{S} is symmetric and uses the lower\n\
+triangular portion of the matrix.\n\
 @end table\n\
 \n\
 @item mode\n\
-The default is to return the Cholesky@tie{}factorization for @var{r}, and if\n\
-@var{mode} is @qcode{'L'}, the conjugate transpose of the\n\
-Cholesky@tie{}factorization is returned.  The conjugate transpose version is\n\
-faster and uses less memory, but returns the same values for @var{count},\n\
-@var{h}, @var{parent} and @var{post} outputs.\n\
+When @var{mode} is unspecified return the Cholesky@tie{}factorization for\n\
+@var{R}.  If @var{mode} is @qcode{\"lower\"} or @qcode{\"L\"} then return\n\
+the conjugate transpose @tcode{@var{R}'} which is a lower triangular factor.\n\
+The conjugate transpose version is faster and uses less memory, but still\n\
+returns the same values for all other outputs: @var{count}, @var{h},\n\
+@var{parent}, and @var{post}.\n\
 @end table\n\
 \n\
-The output variables are\n\
+The output variables are:\n\
 \n\
 @table @var\n\
 @item count\n\
 The row counts of the Cholesky@tie{}factorization as determined by\n\
-@var{typ}.\n\
+@var{typ}.  The computational difficulty of performing the true\n\
+factorization using @code{chol} is @code{sum (@var{count} .^ 2)}.\n\
 \n\
 @item h\n\
 The height of the elimination tree.\n\
@@ -92,40 +96,21 @@
 The elimination tree itself.\n\
 \n\
 @item post\n\
-A sparse boolean matrix whose structure is that of the Cholesky\n\
-factorization as determined by @var{typ}.\n\
+A sparse boolean matrix whose structure is that of the\n\
+Cholesky@tie{}factorization as determined by @var{typ}.\n\
 @end table\n\
+@seealso{chol, etree, treelayout}\n\
 @end deftypefn")
 {
 #ifdef HAVE_CHOLMOD
 
   int nargin = args.length ();
 
-  if (nargin < 1 || nargin > 3 || nargout > 5)
+  if (nargin < 1 || nargin > 3)
     print_usage ();
 
   octave_value_list retval;
-
-  cholmod_common Common;
-  cholmod_common *cm = &Common;
-  CHOLMOD_NAME(start) (cm);
-
-  double spu = octave_sparse_params::get_key ("spumoni");
-  if (spu == 0.)
-    {
-      cm->print = -1;
-      SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
-    }
-  else
-    {
-      cm->print = static_cast<int> (spu) + 2;
-      SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
-    }
-
-  cm->error_handler = &SparseCholError;
-  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
-  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
-
+ 
   double dummy;
   cholmod_sparse Astore;
   cholmod_sparse *A = &Astore;
@@ -175,23 +160,34 @@
 
   if (nargin > 1)
     {
+      std::string str = args(1).xstring_value ("TYP must be a string");
+      // FIXME: The input validation could be improved to use strncmp
       char ch;
-      std::string str = args(1).string_value ();
-      ch = tolower (str.c_str ()[0]);
-      if (ch == 'r')
+      ch = tolower (str[0]);
+      if (ch == 'r')          // 'row'
         A->stype = 0;
-      else if (ch == 'c')
+      else if (ch == 'c')     // 'col'
         {
           n = A->ncol;
           coletree = true;
           A->stype = 0;
         }
-      else if (ch == 's')
+      else if (ch == 's')     // 'sym' (default)
         A->stype = 1;
-      else if (ch == 's')
+      else if (ch == 'l')     // 'lo'
         A->stype = -1;
       else
-        error ("symbfact: unrecognized TYP in symbolic factorization");
+        error ("symbfact: unrecognized TYP \"%s\"", str.c_str ());
+    }
+
+  if (nargin == 3)
+    {
+      std::string str = args(2).xstring_value ("MODE must be a string");
+      // FIXME: The input validation could be improved to use strncmp
+      char ch;
+      ch = toupper (str[0]);
+      if (ch != 'L')
+        error ("symbfact: unrecognized MODE \"%s\"", str.c_str ());
     }
 
   if (A->stype && A->nrow != A->ncol)
@@ -203,33 +199,65 @@
   OCTAVE_LOCAL_BUFFER (octave_idx_type, First, n);
   OCTAVE_LOCAL_BUFFER (octave_idx_type, Level, n);
 
+  cholmod_common Common;
+  cholmod_common *cm = &Common;
+  CHOLMOD_NAME(start) (cm);
+
+  double spu = octave_sparse_params::get_key ("spumoni");
+  if (spu == 0.)
+    {
+      cm->print = -1;
+      SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, 0);
+    }
+  else
+    {
+      cm->print = static_cast<int> (spu) + 2;
+      SUITESPARSE_ASSIGN_FPTR (printf_func, cm->print_function, &SparseCholPrint);
+    }
+
+  cm->error_handler = &SparseCholError;
+  SUITESPARSE_ASSIGN_FPTR2 (divcomplex_func, cm->complex_divide, divcomplex);
+  SUITESPARSE_ASSIGN_FPTR2 (hypot_func, cm->hypotenuse, hypot);
+
   cholmod_sparse *F = CHOLMOD_NAME(transpose) (A, 0, cm);
   cholmod_sparse *Aup, *Alo;
 
   if (A->stype == 1 || coletree)
     {
-      Aup = A ;
-      Alo = F ;
+      Aup = A;
+      Alo = F;
     }
   else
     {
-      Aup = F ;
-      Alo = A ;
+      Aup = F;
+      Alo = A;
     }
 
   CHOLMOD_NAME(etree) (Aup, Parent, cm);
 
+  ColumnVector tmp (n);    // Declaration must precede any goto cleanup.
+  std::string err_msg;
+
   if (cm->status < CHOLMOD_OK)
-    error ("symbfact: matrix corrupted");
+    {
+      err_msg = "symbfact: matrix corrupted";
+      goto cleanup;
+    }
 
   if (CHOLMOD_NAME(postorder) (Parent, n, 0, Post, cm) != n)
-    error ("symbfact: postorder failed");
+    {
+      err_msg = "symbfact: postorder failed";
+      goto cleanup;
+    }
 
   CHOLMOD_NAME(rowcolcounts) (Alo, 0, 0, Parent, Post, 0,
                               ColCount, First, Level, cm);
 
   if (cm->status < CHOLMOD_OK)
-    error ("symbfact: matrix corrupted");
+    {
+      err_msg = "symbfact: matrix corrupted";
+      goto cleanup;
+    }
 
   if (nargout > 4)
     {
@@ -257,11 +285,10 @@
         }
 
       // count the total number of entries in L
-      octave_idx_type lnz = 0 ;
+      octave_idx_type lnz = 0;
       for (octave_idx_type j = 0 ; j < n ; j++)
         lnz += ColCount[j];
 
-
       // allocate the output matrix L (pattern-only)
       SparseBoolMatrix L (n, n, lnz);
 
@@ -274,7 +301,6 @@
         }
       L.xcidx(n) = lnz;
 
-
       // create a copy of the column pointers
       octave_idx_type *W = First;
       for (octave_idx_type j = 0 ; j < n ; j++)
@@ -282,26 +308,25 @@
 
       // get workspace for computing one row of L
       cholmod_sparse *R
-        = CHOLMOD_NAME (allocate_sparse) (n, 1, n, false, true,
-                                          0, CHOLMOD_PATTERN, cm);
-      octave_idx_type *Rp = static_cast<octave_idx_type *>(R->p);
-      octave_idx_type *Ri = static_cast<octave_idx_type *>(R->i);
+        = CHOLMOD_NAME(allocate_sparse) (n, 1, n, false, true,
+                                         0, CHOLMOD_PATTERN, cm);
+      octave_idx_type *Rp = static_cast<octave_idx_type *> (R->p);
+      octave_idx_type *Ri = static_cast<octave_idx_type *> (R->i);
 
       // compute L one row at a time
       for (octave_idx_type k = 0 ; k < n ; k++)
         {
           // get the kth row of L and store in the columns of L
-          CHOLMOD_NAME (row_subtree) (A1, A2, k, Parent, R, cm) ;
+          CHOLMOD_NAME(row_subtree) (A1, A2, k, Parent, R, cm);
           for (octave_idx_type p = 0 ; p < Rp[1] ; p++)
-            L.xridx (W[Ri[p]]++) = k ;
+            L.xridx (W[Ri[p]]++) = k;
 
           // add the diagonal entry
-          L.xridx (W[k]++) = k ;
+          L.xridx (W[k]++) = k;
         }
 
       // free workspace
-      CHOLMOD_NAME (free_sparse) (&R, cm) ;
-
+      CHOLMOD_NAME(free_sparse) (&R, cm);
 
       // transpose L to get R, or leave as is
       if (nargin < 3)
@@ -314,7 +339,6 @@
       retval(4) = L;
     }
 
-  ColumnVector tmp (n);
   if (nargout > 3)
     {
       for (octave_idx_type i = 0; i < n; i++)
@@ -332,10 +356,10 @@
   if (nargout > 1)
     {
       // compute the elimination tree height
-      octave_idx_type height = 0 ;
+      octave_idx_type height = 0;
       for (int i = 0 ; i < n ; i++)
-        height = (height > Level[i] ? height : Level[i]);
-      height++ ;
+        height = std::max (height, Level[i]);
+      height++;
       retval(1) = static_cast<double> (height);
     }
 
@@ -343,9 +367,49 @@
     tmp(i) = ColCount[i];
   retval(0) = tmp;
 
+cleanup:
+  CHOLMOD_NAME(free_sparse) (&F, cm);
+  CHOLMOD_NAME(finish) (cm);
+
+  if (! err_msg.empty ())
+    error (err_msg.c_str ());
+
   return retval;
 
 #else
   err_disabled_feature ("symbfact", "CHOLMOD");
 #endif
 }
+
+/*
+%!testif HAVE_CHOLMOD
+%! A = sparse (magic (3));
+%! [count, h, parent, post, r] = symbfact (A);
+%! assert (count, [3; 2; 1]);
+%! assert (h, 3);
+%! assert (parent, [2; 3; 0]);
+%! assert (r, sparse (triu (true (3))));
+
+%!testif HAVE_CHOLMOD
+%! ## Test MODE "lower"
+%! A = sparse (magic (3));
+%! [~, ~, ~, ~, l] = symbfact (A, "sym", "lower");
+%! assert (l, sparse (tril (true (3))));
+
+%!testif HAVE_CHOLMOD
+%! ## Bug #42587, singular matrix 
+%! A = sparse ([1 0 8;0 1 8;8 8 1]);
+%! [count, h, parent, post, r] = symbfact (A);
+
+## Test input validation
+%!testif HAVE_CHOLMOD
+%! fail ("symbfact ()");
+%! fail ("symbfact (1,2,3,4)");
+%! fail ("symbfact ({1})", "wrong type argument 'cell'");
+%! fail ("symbfact (sparse (1), {1})", "TYP must be a string");
+%! fail ("symbfact (sparse (1), 'foobar')", 'unrecognized TYP "foobar"'); 
+%! fail ("symbfact (sparse (1), 'sym', {'L'})", "MODE must be a string");
+%! fail ('symbfact (sparse (1), "sym", "foobar")', 'unrecognized MODE "foobar"'); 
+%! fail ("symbfact (sparse ([1, 2; 3, 4; 5, 6]))", "S must be a square matrix");
+
+*/