changeset 5440:b73d469ef0c9

[project @ 2005-09-04 12:31:45 by dbateman] ChangeLog
author dbateman
date Sun, 04 Sep 2005 12:31:45 +0000
parents ca7ebff13a16
children c0c81dc78776
files src/DLD-FUNCTIONS/colamd.cc
diffstat 1 files changed, 135 insertions(+), 130 deletions(-) [+]
line wrap: on
line diff
--- a/src/DLD-FUNCTIONS/colamd.cc
+++ b/src/DLD-FUNCTIONS/colamd.cc
@@ -40,36 +40,43 @@
 #include "ov-re-sparse.h"
 #include "ov-cx-sparse.h"
 
-#if SIZEOF_INT == SIZEOF_OCTAVE_IDX_TYPE
-
 // External COLAMD functions in C
 extern "C" {
 #include "COLAMD/colamd.h"
 }
 
+#ifdef IDX_TYPE_LONG
+#define COLAMD_NAME(name) colamd_l ## name
+#define SYMAMD_NAME(name) symamd_l ## name
+#else
+#define COLAMD_NAME(name) colamd ## name
+#define SYMAMD_NAME(name) symamd ## name
+#endif
+
 // The symmetric column elimination tree code take from the Davis LDL code. 
 // Copyright given elsewhere in this file.
 static 
-void symetree (const int *ridx, const int *cidx, int *Parent, int *P, int n)
+void symetree (const octave_idx_type *ridx, const octave_idx_type *cidx, 
+	       octave_idx_type *Parent, octave_idx_type *P, octave_idx_type n)
 {
-  OCTAVE_LOCAL_BUFFER (int, Flag, n);
-  OCTAVE_LOCAL_BUFFER (int, Pinv, (P ? n : 0));
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Flag, n);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, Pinv, (P ? n : 0));
   if (P)
     // If P is present then compute Pinv, the inverse of P
-    for (int k = 0 ; k < n ; k++)
+    for (octave_idx_type k = 0 ; k < n ; k++)
       Pinv [P [k]] = k ;
 
-  for (int k = 0 ; k < n ; k++)
+  for (octave_idx_type k = 0 ; k < n ; k++)
     {
       // L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k)
       Parent [k] = n ;	              // parent of k is not yet known 
       Flag [k] = k ;		      // mark node k as visited 
-      int kk = (P) ? (P [k]) : (k) ;  // kth original, or permuted, column
-      int p2 = cidx [kk+1] ;
-      for (int p = cidx [kk] ; p < p2 ; p++)
+      octave_idx_type kk = (P) ? (P [k]) : (k) ;  // kth original, or permuted, column
+      octave_idx_type p2 = cidx [kk+1] ;
+      for (octave_idx_type p = cidx [kk] ; p < p2 ; p++)
 	{
 	  // A (i,k) is nonzero (original or permuted A)
-	  int i = (Pinv) ? (Pinv [ridx [p]]) : (ridx [p]) ;
+	  octave_idx_type i = (Pinv) ? (Pinv [ridx [p]]) : (ridx [p]) ;
 	  if (i < k)
 	    {
 	      // follow path from i to root of etree, stop at flagged node 
@@ -87,23 +94,23 @@
 
 // The elimination tree post-ordering code below is taken from SuperLU
 static inline
-int make_set (int i, int *pp)
+octave_idx_type make_set (octave_idx_type i, octave_idx_type *pp)
 {
   pp[i] = i;
   return i;
 }
 
 static inline
-int link (int s, int t, int *pp)
+octave_idx_type link (octave_idx_type s, octave_idx_type t, octave_idx_type *pp)
 {
   pp[s] = t;
   return t;
 }
 
 static inline
-int find (int i, int *pp)
+octave_idx_type find (octave_idx_type i, octave_idx_type *pp)
 {
-  register int p, gp;
+  register octave_idx_type p, gp;
     
   p = pp[i];
   gp = pp[p];
@@ -117,9 +124,11 @@
 }
 
 static
-int etdfs (int v, int *first_kid, int *next_kid, int *post, int postnum)
+octave_idx_type etdfs (octave_idx_type v, octave_idx_type *first_kid, 
+		       octave_idx_type *next_kid, octave_idx_type *post, 
+		       octave_idx_type postnum)
 {
-  for (int w = first_kid[v]; w != -1; w = next_kid[w]) {
+  for (octave_idx_type w = first_kid[v]; w != -1; w = next_kid[w]) {
     postnum = etdfs (w, first_kid, next_kid, post, postnum);
   }
   post[postnum++] = v;
@@ -128,17 +137,17 @@
 }
 
 static
-void TreePostorder(int n, int *parent, int *post)
+void TreePostorder(octave_idx_type n, octave_idx_type *parent, octave_idx_type *post)
 {
   // Allocate storage for working arrays and results
-  OCTAVE_LOCAL_BUFFER (int, first_kid, n+1);
-  OCTAVE_LOCAL_BUFFER (int, next_kid, n+1);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, first_kid, n+1);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, next_kid, n+1);
 
   // Set up structure describing children
-  for (int v = 0; v <= n; first_kid[v++] = -1);
-  for (int v = n-1; v >= 0; v--) 
+  for (octave_idx_type v = 0; v <= n; first_kid[v++] = -1);
+  for (octave_idx_type v = n-1; v >= 0; v--) 
     {
-      int dad = parent[v];
+      octave_idx_type dad = parent[v];
       next_kid[v] = first_kid[dad];
       first_kid[dad] = v;
     }
@@ -148,19 +157,20 @@
 }
 
 static 
-void coletree (const int *ridx, const int *colbeg, int *colend, 
-	       int *parent, int nr, int nc)
+void coletree (const octave_idx_type *ridx, const octave_idx_type *colbeg,
+	       octave_idx_type *colend, octave_idx_type *parent, 
+	       octave_idx_type nr, octave_idx_type nc)
 {
-  OCTAVE_LOCAL_BUFFER (int, root, nc);
-  OCTAVE_LOCAL_BUFFER (int, pp, nc);
-  OCTAVE_LOCAL_BUFFER (int, firstcol, nr);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, root, nc);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, pp, nc);
+  OCTAVE_LOCAL_BUFFER (octave_idx_type, firstcol, nr);
 
   // Compute firstcol[row] = first nonzero column in row
-  for (int row = 0; row < nr; firstcol[row++] = nc);
-  for (int col = 0; col < nc; col++) 
-    for (int p = colbeg[col]; p < colend[col]; p++) 
+  for (octave_idx_type row = 0; row < nr; firstcol[row++] = nc);
+  for (octave_idx_type col = 0; col < nc; col++) 
+    for (octave_idx_type p = colbeg[col]; p < colend[col]; p++) 
       {
-	int row = ridx[p];
+	octave_idx_type row = ridx[p];
 	if (firstcol[row] > col)
 	  firstcol[row] = col;
       }
@@ -169,18 +179,18 @@
   // except use (firstcol[r],c) in place of an edge (r,c) of A.
   // Thus each row clique in A'*A is replaced by a star
   // centered at its first vertex, which has the same fill.
-  for (int col = 0; col < nc; col++) 
+  for (octave_idx_type col = 0; col < nc; col++) 
     {
-      int cset = make_set (col, pp);
+      octave_idx_type cset = make_set (col, pp);
       root[cset] = col;
       parent[col] = nc; 
-      for (int p = colbeg[col]; p < colend[col]; p++) 
+      for (octave_idx_type p = colbeg[col]; p < colend[col]; p++) 
 	{
-	  int row = firstcol[ridx[p]];
+	  octave_idx_type row = firstcol[ridx[p]];
 	  if (row >= col) 
 	    continue;
-	  int rset = find (row, pp);
-	  int rroot = root[rset];
+	  octave_idx_type rset = find (row, pp);
+	  octave_idx_type rroot = root[rset];
 	  if (rroot != col) 
 	    {
 	      parent[rroot] = col;
@@ -191,8 +201,6 @@
     }
 }
 
-#endif
-
 DEFUN_DLD (colamd, args, nargout,
     "-*- texinfo -*-\n\
 @deftypefn {Loadable Function} {@var{p} =} colamd (@var{s})\n\
@@ -208,14 +216,15 @@
 (:,@var{p})} also tends to be sparser than that of @code{@var{s}' *\n\
 @var{s}}.\n\
 \n\
-@var{knobs} is an optional two-element input vector. If @var{s} is\n\
-m-by-n, then rows with more than @code{(@var{knobs} (1)) * @var{n}}\n\
-entries are ignored.  Columns with more than @code{(@var{knobs} (2)) *\n\
-@var{m}} entries are removed prior to ordering, and ordered last in the\n\
-output permutation @var{p}. If the knobs parameter is not present, then\n\
-0.5 is used instead, for both @code{@var{knobs} (1)} and\n\
-@code{@var{knobs} (2)}. @code{@var{knobs} (3)} controls the printing of\n\
-statistics and error messages.\n\
+@var{knobs} is an optional one- to three-element input vector.  If @var{s} is\n\
+m-by-n, then rows with more than @code{max(16,@var{knobs}(1)*sqrt(n))} entries\n\
+are ignored. Columns with more than @code{max(16,knobs(2)*sqrt(min(m,n)))}\n\
+entries are removed prior to ordering, and ordered last in the output\n\
+permutation @var{p}. Only completely dense rows or columns are removed\n\
+if @code{@var{knobs} (1)} and @code{@var{knobs} (2)} are < 0, respectively.\n\
+If @code{@var{knobs} (3)} is nonzero, @var{stats} and @var{knobs} are\n\
+printed.  The default is @code{@var{knobs} = [10 10 0]}.  Note that\n\
+@var{knobs} differs from earlier versions of colamd\n\
 \n\
 @var{stats} is an optional 20-element output vector that provides data\n\
 about the ordering and the validity of the input matrix @var{s}. Ordering\n\
@@ -261,9 +270,6 @@
 @seealso{colperm, symamd}")
 {
   octave_value_list retval;
-
-#if SIZEOF_INT == SIZEOF_OCTAVE_IDX_TYPE
-
   int nargin = args.length ();
   int spumoni = 0;
  
@@ -272,8 +278,8 @@
   else
     {
       // Get knobs
-      OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
-      colamd_set_defaults (knobs);
+      OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);      
+      COLAMD_NAME (_set_defaults) (knobs);
 
       // Check for user-passed knobs
       if (nargin == 2)
@@ -282,24 +288,45 @@
 	  int nel_User_knobs = User_knobs.length ();
 	  
 	  if (nel_User_knobs > 0) 
-	    knobs [COLAMD_DENSE_ROW] = User_knobs (COLAMD_DENSE_ROW);
+	    knobs [COLAMD_DENSE_ROW] = User_knobs (0);
 	  if (nel_User_knobs > 1) 
-	    knobs [COLAMD_DENSE_COL] = User_knobs (COLAMD_DENSE_COL) ;
+	    knobs [COLAMD_DENSE_COL] = User_knobs (1) ;
 	  if (nel_User_knobs > 2) 
 	    spumoni = (int) User_knobs (2);
-	}
+
+	  // print knob settings if spumoni is set
+	  if (spumoni)
+	    {
+
+	      octave_stdout << "\ncolamd version " << COLAMD_MAIN_VERSION << "."
+			    <<  COLAMD_SUB_VERSION << ", " << COLAMD_DATE << ":\n";
 
-      // print knob settings if spumoni is set
-      if (spumoni > 0)
-	{
-	  octave_stdout << "colamd: dense row fraction: " 
-			<< knobs [COLAMD_DENSE_ROW] << std::endl;
-	  octave_stdout << "colamd: dense col fraction: " 
-			<< knobs [COLAMD_DENSE_COL] << std::endl;
+	      if (knobs [COLAMD_DENSE_ROW] >= 0)
+		octave_stdout << "knobs(1): " << User_knobs (0) 
+			      << ", rows with > max(16,"
+			      << knobs [COLAMD_DENSE_ROW] << "*sqrt(size(A,2)))"
+			      << " entries removed\n";
+	      else
+		octave_stdout << "knobs(1): " << User_knobs (0)
+			      << ", only completely dense rows removed\n";
+
+	      if (knobs [COLAMD_DENSE_COL] >= 0)
+		octave_stdout << "knobs(2): " << User_knobs (1) 
+			      << ", cols with > max(16,"
+			      << knobs [COLAMD_DENSE_COL] << "*sqrt(size(A)))"
+			      << " entries removed\n";
+	      else
+		octave_stdout << "knobs(2): " << User_knobs (1)
+			      << ", only completely dense columns removed\n";
+
+	      octave_stdout << "knobs(3): " << User_knobs (2) 
+			    << ", statistics and knobs printed\n";
+
+	    }
 	}
       
-      int n_row, n_col, nnz;
-      int *ridx, *cidx;
+      octave_idx_type n_row, n_col, nnz;
+      octave_idx_type *ridx, *cidx;
       SparseComplexMatrix scm;
       SparseMatrix sm;
 
@@ -340,30 +367,30 @@
 	}
 
       // Allocate workspace for colamd
-      OCTAVE_LOCAL_BUFFER (int, p, n_col+1);
-      for (int i = 0; i < n_col+1; i++)
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, p, n_col+1);
+      for (octave_idx_type i = 0; i < n_col+1; i++)
 	p[i] = cidx [i];
 
-      int Alen = colamd_recommended (nnz, n_row, n_col);
-      OCTAVE_LOCAL_BUFFER (int, A, Alen);
-      for (int i = 0; i < nnz; i++)
+      octave_idx_type Alen = COLAMD_NAME (_recommended) (nnz, n_row, n_col);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, A, Alen);
+      for (octave_idx_type i = 0; i < nnz; i++)
 	A[i] = ridx [i];
 
       // Order the columns (destroys A)
-      OCTAVE_LOCAL_BUFFER (int, stats, COLAMD_STATS);
-      if (!colamd (n_row, n_col, Alen, A, p, knobs, stats))
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS);
+      if (! COLAMD_NAME () (n_row, n_col, Alen, A, p, knobs, stats))
 	{
-	  colamd_report (stats) ;
+	  COLAMD_NAME (_report) (stats) ;
 	  error ("colamd: internal error!");
 	  return retval;
 	}
 
       // column elimination tree post-ordering (reuse variables)
-      OCTAVE_LOCAL_BUFFER (int, colbeg, n_col + 1);
-      OCTAVE_LOCAL_BUFFER (int, colend, n_col + 1);
-      OCTAVE_LOCAL_BUFFER (int, etree, n_col + 1);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col + 1);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col + 1);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
 
-      for (int i = 0; i < n_col; i++)
+      for (octave_idx_type i = 0; i < n_col; i++)
 	{
 	  colbeg[i] = cidx[p[i]];
 	  colend[i] = cidx[p[i]+1];
@@ -376,20 +403,20 @@
 
       // return the permutation vector
       NDArray out_perm (dim_vector (1, n_col));
-      for (int i = 0; i < n_col; i++)
+      for (octave_idx_type i = 0; i < n_col; i++)
 	out_perm(i) = p [colbeg [i]] + 1;
 
       retval (0) = out_perm;
 
       // print stats if spumoni > 0
       if (spumoni > 0)
-	colamd_report (stats) ;
+	COLAMD_NAME (_report) (stats) ;
 
       // Return the stats vector
       if (nargout == 2)
 	{
 	  NDArray out_stats (dim_vector (1, COLAMD_STATS));
-	  for (int i = 0 ; i < COLAMD_STATS ; i++)
+	  for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
 	    out_stats (i) = stats [i] ;
 	  retval(1) = out_stats;
 
@@ -401,12 +428,6 @@
 	}
     }
 
-#else
-
-  error ("colamd: not available in this version of Octave");
-
-#endif
-
   return retval;
 }
 
@@ -424,12 +445,13 @@
 symmetric; only the strictly lower triangular part is referenced. @var{s}\n\
 must be square.\n\
 \n\
-@var{knobs} is an optional input argument. If @var{s} is n-by-n, then\n\
-rows and columns with more than @code{@var{knobs} (1) * @var{n}} entries\n\
-are removed prior to ordering, and ordered last in the output permutation\n\
-@var{p}. If the @var{knobs} parameter is not present, then the default of\n\
-0.5 is used instead. @code{@var{knobs} (2)} controls the printing of\n\
-statistics and error messages.\n\
+@var{knobs} is an optional one- to two-element input vector.  If @var{s} is\n\
+n-by-n, then rows and columns with more than\n\
+@code{max(16,@var{knobs}(1)*sqrt(n))} entries are removed prior to ordering,\n\
+and ordered last in the output permutation @var{p}. No rows/columns are\n\
+removed if @code{@var{knobs}(1) < 0}.  If @code{@var{knobs} (2)} is nonzero,\n\
+@code{stats} and @var{knobs} are printed.  The default is @code{@var{knobs} \n\
+= [10 0]}.  Note that @var{knobs} differs from earlier versions of symamd.\n\
 \n\
 @var{stats} is an optional 20-element output vector that provides data\n\
 about the ordering and the validity of the input matrix @var{s}. Ordering\n\
@@ -475,9 +497,6 @@
 @seealso{colperm, colamd}")
 {
   octave_value_list retval;
-
-#if SIZEOF_INT == SIZEOF_OCTAVE_IDX_TYPE
-
   int nargin = args.length ();
   int spumoni = 0;
  
@@ -487,7 +506,7 @@
     {
       // Get knobs
       OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
-      colamd_set_defaults (knobs);
+      COLAMD_NAME (_set_defaults) (knobs);
 
       // Check for user-passed knobs
       if (nargin == 2)
@@ -506,8 +525,8 @@
 	octave_stdout << "symamd: dense row/col fraction: " 
 		      << knobs [COLAMD_DENSE_ROW] << std::endl;
       
-      int n_row, n_col, nnz;
-      int *ridx, *cidx;
+      octave_idx_type n_row, n_col, nnz;
+      octave_idx_type *ridx, *cidx;
       SparseMatrix sm;
       SparseComplexMatrix scm;
 
@@ -553,39 +572,39 @@
 	}
 
       // Allocate workspace for symamd
-      OCTAVE_LOCAL_BUFFER (int, perm, n_col+1);
-      OCTAVE_LOCAL_BUFFER (int, stats, COLAMD_STATS);
-      if (!symamd (n_col, ridx, cidx, perm, knobs, stats, &calloc, &free))
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, perm, n_col+1);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, stats, COLAMD_STATS);
+      if (!SYMAMD_NAME () (n_col, ridx, cidx, perm, knobs, stats, &calloc, &free))
 	{
-	  symamd_report (stats) ;
+	  SYMAMD_NAME (_report) (stats) ;
 	  error ("symamd: internal error!") ;
 	  return retval;
 	}
 
       // column elimination tree post-ordering
-      OCTAVE_LOCAL_BUFFER (int, etree, n_col + 1);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
       symetree (ridx, cidx, etree, perm, n_col);
 
       // Calculate the tree post-ordering
-      OCTAVE_LOCAL_BUFFER (int, post, n_col + 1);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
       TreePostorder (n_col, etree, post);
 
       // return the permutation vector
       NDArray out_perm (dim_vector (1, n_col));
-      for (int i = 0; i < n_col; i++)
+      for (octave_idx_type i = 0; i < n_col; i++)
 	out_perm(i) = perm [post [i]] + 1;
 
       retval (0) = out_perm;
 
       // print stats if spumoni > 0
       if (spumoni > 0)
-	symamd_report (stats) ;
+	SYMAMD_NAME (_report) (stats) ;
 
       // Return the stats vector
       if (nargout == 2)
 	{
 	  NDArray out_stats (dim_vector (1, COLAMD_STATS));
-	  for (int i = 0 ; i < COLAMD_STATS ; i++)
+	  for (octave_idx_type i = 0 ; i < COLAMD_STATS ; i++)
 	    out_stats (i) = stats [i] ;
 	  retval(1) = out_stats;
 
@@ -597,12 +616,6 @@
 	}
     }
 
-#else
-
-  error ("symamd: not available in this version of Octave");
-
-#endif
-
   return retval;
 }
 
@@ -624,16 +637,14 @@
 {
   octave_value_list retval;
 
-#if SIZEOF_INT == SIZEOF_OCTAVE_IDX_TYPE
-
   int nargin = args.length ();
 
   if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2)
     usage ("etree: incorrect number of input and/or output arguments");
   else
     {
-      int n_row, n_col, nnz;
-      int *ridx, *cidx;
+      octave_idx_type n_row, n_col, nnz;
+      octave_idx_type *ridx, *cidx;
       bool is_sym = true;
       SparseMatrix sm;
       SparseComplexMatrix scm;
@@ -680,7 +691,7 @@
 	  }
 
       // column elimination tree post-ordering (reuse variables)
-      OCTAVE_LOCAL_BUFFER (int, etree, n_col + 1);
+      OCTAVE_LOCAL_BUFFER (octave_idx_type, etree, n_col + 1);
       
 
       if (is_sym)
@@ -694,10 +705,10 @@
 	}
       else
 	{
-	  OCTAVE_LOCAL_BUFFER (int, colbeg, n_col);
-	  OCTAVE_LOCAL_BUFFER (int, colend, n_col);
+	  OCTAVE_LOCAL_BUFFER (octave_idx_type, colbeg, n_col);
+	  OCTAVE_LOCAL_BUFFER (octave_idx_type, colend, n_col);
 
-	  for (int i = 0; i < n_col; i++)
+	  for (octave_idx_type i = 0; i < n_col; i++)
 	    {
 	      colbeg[i] = cidx[i];
 	      colend[i] = cidx[i+1];
@@ -707,7 +718,7 @@
 	}
 
       NDArray tree (dim_vector (1, n_col));
-      for (int i = 0; i < n_col; i++)
+      for (octave_idx_type i = 0; i < n_col; i++)
 	// We flag a root with n_col while Matlab does it with zero
 	// Convert for matlab compatiable output
 	if (etree[i] == n_col)
@@ -720,23 +731,17 @@
       if (nargout == 2)
 	{
 	  // Calculate the tree post-ordering
-	  OCTAVE_LOCAL_BUFFER (int, post, n_col + 1);
+	  OCTAVE_LOCAL_BUFFER (octave_idx_type, post, n_col + 1);
 	  TreePostorder (n_col, etree, post);
 
 	  NDArray postorder (dim_vector (1, n_col));
-	  for (int i = 0; i < n_col; i++)
+	  for (octave_idx_type i = 0; i < n_col; i++)
 	    postorder (i) = post[i] + 1;
 
 	  retval (1) = postorder;
 	}
     }
 
-#else
-
-  error ("etree: not available in this version of Octave");
-
-#endif
-
   return retval;
 }