--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/DLD-FUNCTIONS/colamd.cc	Fri Feb 25 19:55:28 2005 +0000
@@ -0,0 +1,726 @@
+/*
+
+Copyright (C) 2004 David Bateman
+Copyright (C) 1998-2004 Andy Adler
+
+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 2, 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 this program; see the file COPYING.  If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+
+*/
+
+// This is the octave interface to colamd, which bore the copyright given
+// in the help of the functions.
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <cstdlib>
+
+#include <string>
+#include <vector>
+
+#include "ov.h"
+#include "defun-dld.h"
+#include "pager.h"
+#include "ov-re-mat.h"
+
+#include "ov-re-sparse.h"
+#include "ov-cx-sparse.h"
+
+// External COLAMD functions in C
+extern "C" {
+#include "COLAMD/colamd.h"
+}
+
+// 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)
+{
+  OCTAVE_LOCAL_BUFFER (int, Flag, n);
+  OCTAVE_LOCAL_BUFFER (int, Pinv, (P ? n : 0));
+  if (P)
+    // If P is present then compute Pinv, the inverse of P
+    for (int k = 0 ; k < n ; k++)
+      Pinv [P [k]] = k ;
+
+  for (int 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++)
+	{
+	  // A (i,k) is nonzero (original or permuted A)
+	  int i = (Pinv) ? (Pinv [ridx [p]]) : (ridx [p]) ;
+	  if (i < k)
+	    {
+	      // follow path from i to root of etree, stop at flagged node 
+	      for ( ; Flag [i] != k ; i = Parent [i])
+		{
+		  // find parent of i if not yet determined
+		  if (Parent [i] == n)
+		    Parent [i] = k ;
+		  Flag [i] = k ;	// mark i as visited
+		}
+	    }
+	}
+    }
+}
+
+// The elimination tree post-ordering code below is taken from SuperLU
+static inline
+int make_set (int i, int *pp)
+{
+  pp[i] = i;
+  return i;
+}
+
+static inline
+int link (int s, int t, int *pp)
+{
+  pp[s] = t;
+  return t;
+}
+
+static inline
+int find (int i, int *pp)
+{
+  register int p, gp;
+    
+  p = pp[i];
+  gp = pp[p];
+  while (gp != p) {
+    pp[i] = gp;
+    i = gp;
+    p = pp[i];
+    gp = pp[p];
+  }
+  return (p);
+}
+
+static
+int etdfs (int v, int *first_kid, int *next_kid, int *post, int postnum)
+{
+  for (int w = first_kid[v]; w != -1; w = next_kid[w]) {
+    postnum = etdfs (w, first_kid, next_kid, post, postnum);
+  }
+  post[postnum++] = v;
+
+  return postnum;
+}
+
+static
+void TreePostorder(int n, int *parent, int *post)
+{
+  // Allocate storage for working arrays and results
+  OCTAVE_LOCAL_BUFFER (int, first_kid, n+1);
+  OCTAVE_LOCAL_BUFFER (int, 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--) 
+    {
+      int dad = parent[v];
+      next_kid[v] = first_kid[dad];
+      first_kid[dad] = v;
+    }
+
+  // Depth-first search from dummy root vertex #n
+  etdfs (n, first_kid, next_kid, post, 0);
+}
+
+static 
+void coletree (const int *ridx, const int *colbeg, int *colend, 
+	       int *parent, int nr, int nc)
+{
+  OCTAVE_LOCAL_BUFFER (int, root, nc);
+  OCTAVE_LOCAL_BUFFER (int, pp, nc);
+  OCTAVE_LOCAL_BUFFER (int, 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++) 
+      {
+	int row = ridx[p];
+	if (firstcol[row] > col)
+	  firstcol[row] = col;
+      }
+
+  // Compute etree by Liu's algorithm for symmetric matrices,
+  // 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++) 
+    {
+      int cset = make_set (col, pp);
+      root[cset] = col;
+      parent[col] = nc; 
+      for (int p = colbeg[col]; p < colend[col]; p++) 
+	{
+	  int row = firstcol[ridx[p]];
+	  if (row >= col) 
+	    continue;
+	  int rset = find (row, pp);
+	  int rroot = root[rset];
+	  if (rroot != col) 
+	    {
+	      parent[rroot] = col;
+	      cset = link (cset, rset, pp);
+	      root[cset] = col;
+	    }
+	}
+    }
+}
+
+DEFUN_DLD (colamd, args, nargout,
+    "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{p} =} colamd (@var{s})\n\
+@deftypefnx {Loadable Function} {@var{p} =} colamd (@var{s}, @var{knobs})\n\
+@deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} colamd (@var{s})\n\
+@deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} colamd (@var{s}, @var{knobs})\n\
+\n\
+Column approximate minimum degree permutation. @code{@var{p} = colamd\n\
+(@var{s})} returns the column approximate minimum degree permutation\n\
+vector for the sparse matrix @var{s}. For a non-symmetric matrix @var{s},\n\
+@code{@var{s} (:,@var{p})} tends to have sparser LU factors than @var{s}.\n\
+The Cholesky factorization of @code{@var{s} (:,@var{p})' * @var{s}\n\
+(:,@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\
+\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\
+statistics are in @code{@var{stats} (1:3)}. @code{@var{stats} (1)} and\n\
+@code{@var{stats} (2)} are the number of dense or empty rows and columns\n\
+ignored by COLAMD and @code{@var{stats} (3)} is the number of garbage\n\
+collections performed on the internal data structure used by COLAMD\n\
+(roughly of size @code{2.2 * nnz(@var{s}) + 4 * @var{m} + 7 * @var{n}}\n\
+integers).\n\
+\n\
+Octave built-in functions are intended to generate valid sparse matrices,\n\
+with no duplicate entries, with ascending row indices of the nonzeros\n\
+in each column, with a non-negative number of entries in each column (!)\n\
+and so on.  If a matrix is invalid, then COLAMD may or may not be able\n\
+to continue.  If there are duplicate entries (a row index appears two or\n\
+more times in the same column) or if the row indices in a column are out\n\
+of order, then COLAMD can correct these errors by ignoring the duplicate\n\
+entries and sorting each column of its internal copy of the matrix\n\
+@var{s} (the input matrix @var{s} is not repaired, however). If a matrix\n\
+is invalid in other ways then COLAMD cannot continue, an error message is\n\
+printed, and no output arguments (@var{p} or @var{stats}) are returned.\n\
+COLAMD is thus a simple way to check a sparse matrix to see if it's\n\
+valid.\n\
+\n\
+@code{@var{stats} (4:7)} provide information if COLAMD was able to\n\
+continue. The matrix is OK if @code{@var{stats} (4)} is zero, or 1 if\n\
+invalid. @code{@var{stats} (5)} is the rightmost column index that is\n\
+unsorted or contains duplicate entries, or zero if no such column exists.\n\
+@code{@var{stats} (6)} is the last seen duplicate or out-of-order row\n\
+index in the column index given by @code{@var{stats} (5)}, or zero if no\n\
+such row index exists. @code{@var{stats} (7)} is the number of duplicate\n\
+or out-of-order row indices. @code{@var{stats} (8:20)} is always zero in\n\
+the current version of COLAMD (reserved for future use).\n\
+\n\
+The ordering is followed by a column elimination tree post-ordering.\n\
+\n\
+The authors of the code itself are Stefan I. Larimore and Timothy A.\n\
+Davis (davis@@cise.ufl.edu), University of Florida.  The algorithm was\n\
+developed in collaboration with John Gilbert, Xerox PARC, and Esmond\n\
+Ng, Oak Ridge National Laboratory. (see\n\
+@url{http://www.cise.ufl.edu/research/sparse/colamd})\n\
+@end deftypefn\n\
+@seealso{colperm, symamd}")
+{
+  octave_value_list retval;
+  int nargin = args.length ();
+  int spumoni = 0;
+ 
+  if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2)
+    usage ("colamd: incorrect number of input and/or output arguments");
+  else
+    {
+      // Get knobs
+      OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
+      colamd_set_defaults (knobs);
+
+      // Check for user-passed knobs
+      if (nargin == 2)
+	{
+	  NDArray User_knobs = args(1).array_value ();
+	  int nel_User_knobs = User_knobs.length ();
+	  
+	  if (nel_User_knobs > 0) 
+	    knobs [COLAMD_DENSE_ROW] = User_knobs (COLAMD_DENSE_ROW);
+	  if (nel_User_knobs > 1) 
+	    knobs [COLAMD_DENSE_COL] = User_knobs (COLAMD_DENSE_COL) ;
+	  if (nel_User_knobs > 2) 
+	    spumoni = (int) User_knobs (2);
+	}
+
+      // 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;
+	}
+      
+      int n_row, n_col, nnz;
+      int *ridx, *cidx;
+      SparseComplexMatrix scm;
+      SparseMatrix sm;
+
+      if (args(0).class_name () == "sparse")
+	{
+	  if (args(0).is_complex_type ())
+	    {
+	      scm = args(0). sparse_complex_matrix_value ();
+	      n_row = scm.rows ();
+	      n_col = scm.cols ();
+	      nnz = scm.nnz ();
+	      ridx = scm.xridx ();
+	      cidx = scm.xcidx ();
+	    }
+	  else
+	    {
+	      sm = args(0).sparse_matrix_value ();
+
+	      n_row = sm.rows ();
+	      n_col = sm.cols ();
+	      nnz = sm.nnz ();
+	      ridx = sm.xridx ();
+	      cidx = sm.xcidx ();
+	    }
+	}
+      else
+	{
+	  if (args(0).is_complex_type ())
+	    sm = SparseMatrix (real (args(0).complex_matrix_value ()));
+	  else
+	    sm = SparseMatrix (args(0).matrix_value ());
+
+	  n_row = sm.rows ();
+	  n_col = sm.cols ();
+	  nnz = sm.nnz ();
+	  ridx = sm.xridx ();
+	  cidx = sm.xcidx ();
+	}
+
+      // Allocate workspace for colamd
+      OCTAVE_LOCAL_BUFFER (int, p, n_col+1);
+      for (int 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++)
+	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))
+	{
+	  colamd_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);
+
+      for (int i = 0; i < n_col; i++)
+	{
+	  colbeg[i] = cidx[p[i]];
+	  colend[i] = cidx[p[i]+1];
+	}
+
+      coletree (ridx, colbeg, colend, etree, n_row, n_col);
+
+      // Calculate the tree post-ordering
+      TreePostorder (n_col, etree, colbeg);
+
+      // return the permutation vector
+      NDArray out_perm (dim_vector (1, n_col));
+      for (int 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) ;
+
+      // Return the stats vector
+      if (nargout == 2)
+	{
+	  NDArray out_stats (dim_vector (1, COLAMD_STATS));
+	  for (int i = 0 ; i < COLAMD_STATS ; i++)
+	    out_stats (i) = stats [i] ;
+	  retval(1) = out_stats;
+
+	  // fix stats (5) and (6), for 1-based information on 
+	  // jumbled matrix.  note that this correction doesn't 
+	  // occur if symamd returns FALSE
+	  out_stats (COLAMD_INFO1) ++ ; 
+	  out_stats (COLAMD_INFO2) ++ ; 
+	}
+    }
+
+  return retval;
+}
+
+DEFUN_DLD (symamd, args, nargout,
+    "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{p} =} symamd (@var{s})\n\
+@deftypefnx {Loadable Function} {@var{p} =} symamd (@var{s}, @var{knobs})\n\
+@deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} symamd (@var{s})\n\
+@deftypefnx {Loadable Function} {[@var{p}, @var{stats}] =} symamd (@var{s}, @var{knobs})\n\
+\n\
+For a symmetric positive definite matrix @var{s}, returns the permutation\n\
+vector p such that @code{@var{s} (@var{p}, @var{p})} tends to have a\n\
+sparser Cholesky factor than @var{s}. Sometimes SYMAMD works well for\n\
+symmetric indefinite matrices too. The matrix @var{s} is assumed to be\n\
+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\
+\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\
+statistics are in @code{@var{stats} (1:3)}. @code{@var{stats} (1) =\n\
+@var{stats} (2)} is the number of dense or empty rows and columns\n\
+ignored by SYMAMD and @code{@var{stats} (3)} is the number of garbage\n\
+collections performed on the internal data structure used by SYMAMD\n\
+(roughly of size @code{8.4 * nnz (tril (@var{s}, -1)) + 9 * @var{n}}\n\
+integers).\n\
+\n\
+Octave built-in functions are intended to generate valid sparse matrices,\n\
+with no duplicate entries, with ascending row indices of the nonzeros\n\
+in each column, with a non-negative number of entries in each column (!)\n\
+and so on.  If a matrix is invalid, then SYMAMD may or may not be able\n\
+to continue.  If there are duplicate entries (a row index appears two or\n\
+more times in the same column) or if the row indices in a column are out\n\
+of order, then SYMAMD can correct these errors by ignoring the duplicate\n\
+entries and sorting each column of its internal copy of the matrix S (the\n\
+input matrix S is not repaired, however).  If a matrix is invalid in\n\
+other ways then SYMAMD cannot continue, an error message is printed, and\n\
+no output arguments (@var{p} or @var{stats}) are returned.  SYMAMD is\n\
+thus a simple way to check a sparse matrix to see if it's valid.\n\
+\n\
+@code{@var{stats} (4:7)} provide information if SYMAMD was able to\n\
+continue. The matrix is OK if @code{@var{stats} (4)} is zero, or 1\n\
+if invalid. @code{@var{stats} (5)} is the rightmost column index that\n\
+is unsorted or contains duplicate entries, or zero if no such column\n\
+exists. @code{@var{stats} (6)} is the last seen duplicate or out-of-order\n\
+row index in the column index given by @code{@var{stats} (5)}, or zero\n\
+if no such row index exists. @code{@var{stats} (7)} is the number of\n\
+duplicate or out-of-order row indices. @code{@var{stats} (8:20)} is\n\
+always zero in the current version of SYMAMD (reserved for future use).\n\
+\n\
+The ordering is followed by a column elimination tree post-ordering.\n\
+\n\
+\n\
+The authors of the code itself are Stefan I. Larimore and Timothy A.\n\
+Davis (davis@@cise.ufl.edu), University of Florida.  The algorithm was\n\
+developed in collaboration with John Gilbert, Xerox PARC, and Esmond\n\
+Ng, Oak Ridge National Laboratory. (see\n\
+@url{http://www.cise.ufl.edu/research/sparse/colamd})\n\
+@end deftypefn\n\
+@seealso{colperm, colamd}")
+{
+  octave_value_list retval;
+  int nargin = args.length ();
+  int spumoni = 0;
+ 
+  if (nargout < 0 || nargout > 2 || nargin < 0 || nargin > 2)
+    usage ("symamd: incorrect number of input and/or output arguments");
+  else
+    {
+      // Get knobs
+      OCTAVE_LOCAL_BUFFER (double, knobs, COLAMD_KNOBS);
+      colamd_set_defaults (knobs);
+
+      // Check for user-passed knobs
+      if (nargin == 2)
+	{
+	  NDArray User_knobs = args(1).array_value ();
+	  int nel_User_knobs = User_knobs.length ();
+	  
+	  if (nel_User_knobs > 0) 
+	    knobs [COLAMD_DENSE_ROW] = User_knobs (COLAMD_DENSE_ROW);
+	  if (nel_User_knobs > 1) 
+	    spumoni = (int) User_knobs (1);
+	}
+
+      // print knob settings if spumoni is set
+      if (spumoni > 0)
+	octave_stdout << "symamd: dense row/col fraction: " 
+		      << knobs [COLAMD_DENSE_ROW] << std::endl;
+      
+      int n_row, n_col, nnz;
+      int *ridx, *cidx;
+      SparseMatrix sm;
+      SparseComplexMatrix scm;
+
+      if (args(0).class_name () == "sparse")
+	{
+	  if (args(0).is_complex_type ())
+	    {
+	      scm = args(0).sparse_complex_matrix_value ();
+	      n_row = scm.rows ();
+	      n_col = scm.cols ();
+	      nnz = scm.nnz ();
+	      ridx = scm.xridx ();
+	      cidx = scm.xcidx ();
+	    }
+	  else
+	    {
+	      sm = args(0).sparse_matrix_value ();
+	      n_row = sm.rows ();
+	      n_col = sm.cols ();
+	      nnz = sm.nnz ();
+	      ridx = sm.xridx ();
+	      cidx = sm.xcidx ();
+	    }
+	}
+      else
+	{
+	  if (args(0).is_complex_type ())
+	    sm = SparseMatrix (real (args(0).complex_matrix_value ()));
+	  else
+	    sm = SparseMatrix (args(0).matrix_value ());
+	  
+	  n_row = sm.rows ();
+	  n_col = sm.cols ();
+	  nnz = sm.nnz ();
+	  ridx = sm.xridx ();
+	  cidx = sm.xcidx ();
+	}
+
+      if (n_row != n_col)
+	{
+	  error ("symamd: matrix must be square");
+	  return retval;
+	}
+
+      // 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))
+	{
+	  symamd_report (stats) ;
+	  error ("symamd: internal error!") ;
+	  return retval;
+	}
+
+      // column elimination tree post-ordering
+      OCTAVE_LOCAL_BUFFER (int, etree, n_col + 1);
+      symetree (ridx, cidx, etree, perm, n_col);
+
+      // Calculate the tree post-ordering
+      OCTAVE_LOCAL_BUFFER (int, 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++)
+	out_perm(i) = perm [post [i]] + 1;
+
+      retval (0) = out_perm;
+
+      // print stats if spumoni > 0
+      if (spumoni > 0)
+	symamd_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++)
+	    out_stats (i) = stats [i] ;
+	  retval(1) = out_stats;
+
+	  // fix stats (5) and (6), for 1-based information on 
+	  // jumbled matrix.  note that this correction doesn't 
+	  // occur if symamd returns FALSE
+	  out_stats (COLAMD_INFO1) ++ ; 
+	  out_stats (COLAMD_INFO2) ++ ; 
+	}
+    }
+
+  return retval;
+}
+
+DEFUN_DLD (etree, args, nargout,
+    "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{p} =} etree (@var{s})\n\
+@deftypefnx {Loadable Function} {@var{p} =} etree (@var{s}, @var{typ})\n\
+@deftypefnx {Loadable Function} {[@var{p}, @var{q}] =} etree (@var{s}, @var{typ})\n\
+\n\
+Returns the elimination tree for the matrix @var{s}. By default @var{s}\n\
+is assumed to be symmetric and the symmetric elimination tree is\n\
+returned. The argument @var{typ} controls whether a symmetric or\n\
+column elimination tree is returned. Valid values of @var{typ} are\n\
+'sym' or 'col', for symmetric or column elimination tree respectively\n\
+\n\
+Called with a second argument, @dfn{etree} also returns the postorder\n\
+permutations on the tree.\n\
+@end deftypefn")
+{
+  octave_value_list retval;
+  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;
+      bool is_sym = true;
+      SparseMatrix sm;
+      SparseComplexMatrix scm;
+
+      if (args(0).class_name () == "sparse")
+	{
+	  if (args(0).is_complex_type ())
+	    {
+	      scm = args(0).sparse_complex_matrix_value ();
+	      n_row = scm.rows ();
+	      n_col = scm.cols ();
+	      nnz = scm.nnz ();
+	      ridx = scm.xridx ();
+	      cidx = scm.xcidx ();
+	    }
+	  else
+	    {
+	      sm = args(0).sparse_matrix_value ();
+	      n_row = sm.rows ();
+	      n_col = sm.cols ();
+	      nnz = sm.nnz ();
+	      ridx = sm.xridx ();
+	      cidx = sm.xcidx ();
+	    }
+
+	}
+      else
+	{
+	  error ("etree: must be called with a sparse matrix");
+	  return retval;
+	}
+
+      if (nargin == 2)
+	if (args(1).is_string ())
+	  {
+	    std::string str = args(1).string_value ();
+	    if (str.find("C") == 0 || str.find("c") == 0)
+	      is_sym = false;
+	  }
+	else
+	  {
+	    error ("etree: second argument must be a string");
+	    return retval;
+	  }
+
+      // column elimination tree post-ordering (reuse variables)
+      OCTAVE_LOCAL_BUFFER (int, etree, n_col + 1);
+      
+
+      if (is_sym)
+	{
+	  if (n_row != n_col)
+	    {
+	      error ("etree: matrix is marked as symmetric, but not square");
+	      return retval;
+	    }
+	  symetree (ridx, cidx, etree, NULL, n_col);
+	}
+      else
+	{
+	  OCTAVE_LOCAL_BUFFER (int, colbeg, n_col);
+	  OCTAVE_LOCAL_BUFFER (int, colend, n_col);
+
+	  for (int i = 0; i < n_col; i++)
+	    {
+	      colbeg[i] = cidx[i];
+	      colend[i] = cidx[i+1];
+	    }
+
+	  coletree (ridx, colbeg, colend, etree, n_row, n_col);
+	}
+
+      NDArray tree (dim_vector (1, n_col));
+      for (int 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)
+	  tree (i) = 0;
+	else
+	  tree (i) = etree[i] + 1;
+
+      retval (0) = tree;
+
+      if (nargout == 2)
+	{
+	  // Calculate the tree post-ordering
+	  OCTAVE_LOCAL_BUFFER (int, post, n_col + 1);
+	  TreePostorder (n_col, etree, post);
+
+	  NDArray postorder (dim_vector (1, n_col));
+	  for (int i = 0; i < n_col; i++)
+	    postorder (i) = post[i] + 1;
+
+	  retval (1) = postorder;
+	}
+    }
+
+  return retval;
+}
+
+DEFUN_DLD (symbfact, args, nargout,
+    "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{r}]} = symbfact (@var{s}, @var{typ})\n\
+\n\
+Performs a symbolic factorization analysis on the sparse matrix @var{s}.\n\
+@end deftypefn")
+{
+  error ("symbfact: not implemented yet");
+  return octave_value ();
+}
+
+/*
+;;; Local Variables: ***
+;;; mode: C++ ***
+;;; End: ***
+*/