# HG changeset patch # User Rik # Date 1459785622 25200 # Node ID bad3ed83330d0967722a16f8f8d6341d624244ea # Parent ee1a009dd60f87c7881be8118bf2caf8eeb2e78d 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. diff -r ee1a009dd60f -r bad3ed83330d libinterp/dldfcn/symbfact.cc --- 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 (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 (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(R->p); - octave_idx_type *Ri = static_cast(R->i); + = CHOLMOD_NAME(allocate_sparse) (n, 1, n, false, true, + 0, CHOLMOD_PATTERN, cm); + octave_idx_type *Rp = static_cast (R->p); + octave_idx_type *Ri = static_cast (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 (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"); + +*/