changeset 6066:27fa671e5599

[project @ 2006-10-20 08:33:19 by dbateman]
author dbateman
date Fri, 20 Oct 2006 08:33:20 +0000
parents 814f20da2cdb
children b83e4f9540a0
files src/ChangeLog src/DLD-FUNCTIONS/spqr.cc
diffstat 2 files changed, 100 insertions(+), 34 deletions(-) [+]
line wrap: on
line diff
--- a/src/ChangeLog	Fri Oct 20 03:03:12 2006 +0000
+++ b/src/ChangeLog	Fri Oct 20 08:33:20 2006 +0000
@@ -1,3 +1,11 @@
+2006-10-20  David Bateman  <dbateman@free.fr>
+
+	* DLD-FUNCTION/spqr.cc (dmperm_internal): New function with core
+	of Fdmperm.
+	(Fdmperm): Call dmperm_internal rather then calculating loally.
+	(Fsprank): New function to calculate the strutural rank also using
+	dmperm_internal.
+
 2006-10-19  John W. Eaton  <jwe@octave.org>
 
 	* ov-struct.cc (octave_struct::as_mxArrary):
--- a/src/DLD-FUNCTIONS/spqr.cc	Fri Oct 20 03:03:12 2006 +0000
+++ b/src/DLD-FUNCTIONS/spqr.cc	Fri Oct 20 08:33:20 2006 +0000
@@ -230,39 +230,11 @@
   return ret;
 }
 
-DEFUN_DLD (dmperm, args, nargout,
-  "-*- texinfo -*-\n\
-@deftypefn {Loadable Function} {@var{p} =} dmperm (@var{s})\n\
-@deftypefnx {Loadable Function} {[@var{p}, @var{q}, @var{r}, @var{s}] =} dmperm (@var{s})\n\
-\n\
-@cindex Dulmage-Mendelsohn decomposition\n\
-Perform a Deulmage-Mendelsohn permutation on the sparse matrix @var{s}.\n\
-With a single output argument @dfn{dmperm} performs the row permutations\n\
-@var{p} such that @code{@var{s} (@var{p},:)} has no zero elements on the\n\
-diagonal.\n\
-\n\
-Called with two or more output arguments, returns the row and column\n\
-permutations, such that @code{@var{s} (@var{p}, @var{q})} is in block\n\
-triangular form. The values of @var{r} and @var{s} define the boundaries\n\
-of the blocks. If @var{s} is square then @code{@var{r} == @var{s}}.\n\
-\n\
-The method used is described in: A. Pothen & C.-J. Fan. Computing the block\n\
-triangular form of a sparse matrix. ACM Trans. Math. Software,\n\
-16(4):303-324, 1990.\n\
-@seealso{colamd, ccolamd}\n\
-@end deftypefn")
+#if HAVE_CXSPARSE
+static octave_value_list
+dmperm_internal (bool rank, const octave_value arg, int nargout)
 {
-  int nargin = args.length();
   octave_value_list retval;
-  
-#if HAVE_CXSPARSE
-  if (nargin != 1)
-    {
-      print_usage ();
-      return retval;
-    }
-
-  octave_value arg = args(0);
   octave_idx_type nr = arg.rows ();
   octave_idx_type nc = arg.columns ();
   SparseMatrix m;
@@ -290,14 +262,23 @@
 
   if (!error_state)
     {
-      if (nargout <= 1)
+      if (nargout <= 1 || rank)
 	{
 #if defined(CS_VER) && (CS_VER >= 2)
 	  octave_idx_type *jmatch = CXSPARSE_NAME (_maxtrans) (&csm, 0);
 #else
 	  octave_idx_type *jmatch = CXSPARSE_NAME (_maxtrans) (&csm);
 #endif
-	  retval(0) = put_int (jmatch + nr, nc);
+	  if (rank)
+	    {
+	      octave_idx_type r = 0;
+	      for (octave_idx_type i = 0; i < nc; i++)
+		if (jmatch[nr+i] >= 0)
+		  r++;
+	      retval(0) = static_cast<double>(r);
+	    }
+	  else
+	    retval(0) = put_int (jmatch + nr, nc);
 	  CXSPARSE_NAME (_free) (jmatch);
 	}
       else
@@ -307,6 +288,7 @@
 #else
 	  CXSPARSE_NAME (d) *dm = CXSPARSE_NAME(_dmperm) (&csm);
 #endif
+
 	  //retval(5) = put_int (dm->rr, 5);
 	  //retval(4) = put_int (dm->cc, 5);
 #if defined(CS_VER) && (CS_VER >= 2)
@@ -323,6 +305,43 @@
 	  CXSPARSE_NAME (_dfree) (dm);
 	}
     }
+  return retval;
+}
+#endif
+
+DEFUN_DLD (dmperm, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{p} =} dmperm (@var{s})\n\
+@deftypefnx {Loadable Function} {[@var{p}, @var{q}, @var{r}, @var{s}] =} dmperm (@var{s})\n\
+\n\
+@cindex Dulmage-Mendelsohn decomposition\n\
+Perform a Dulmage-Mendelsohn permutation on the sparse matrix @var{s}.\n\
+With a single output argument @dfn{dmperm} performs the row permutations\n\
+@var{p} such that @code{@var{s} (@var{p},:)} has no zero elements on the\n\
+diagonal.\n\
+\n\
+Called with two or more output arguments, returns the row and column\n\
+permutations, such that @code{@var{s} (@var{p}, @var{q})} is in block\n\
+triangular form. The values of @var{r} and @var{s} define the boundaries\n\
+of the blocks. If @var{s} is square then @code{@var{r} == @var{s}}.\n\
+\n\
+The method used is described in: A. Pothen & C.-J. Fan. Computing the block\n\
+triangular form of a sparse matrix. ACM Trans. Math. Software,\n\
+16(4):303-324, 1990.\n\
+@seealso{colamd, ccolamd}\n\
+@end deftypefn")
+{
+  int nargin = args.length();
+  octave_value_list retval;
+  
+  if (nargin != 1)
+    {
+      print_usage ();
+      return retval;
+    }
+
+#if HAVE_CXSPARSE
+  retval = dmperm_internal (false, args(0), nargout);
 #else
   error ("dmperm: not available in this version of Octave");
 #endif
@@ -330,7 +349,7 @@
   return retval;
 }
 
-/*
+/* 
 
 %!test
 %! n=20;
@@ -347,6 +366,45 @@
 
 */
 
+DEFUN_DLD (sprank, args, nargout,
+  "-*- texinfo -*-\n\
+@deftypefn {Loadable Function} {@var{p} =} sprank (@var{s})\n\
+\n\
+@cindex Structural Rank\n\
+Calculates the structural rank of a sparse matrix @var{s}. Note that\n\
+only the structure of the matrix is used in this calculation based on\n\
+a Dulmage-Mendelsohn to block triangular form. As such the numerical\n\
+rank of the matrix @var{s} is bounded by @code{sprank (@var{s}) >=\n\
+rank (@var{s})}. Ignoring floating point errors @code{sprank (@var{s}) ==\n\
+rank (@var{s})}.\n\
+@seealso{dmperm}\n\
+@end deftypefn")
+{
+  int nargin = args.length();
+  octave_value_list retval;
+  
+  if (nargin != 1)
+    {
+      print_usage ();
+      return retval;
+    }
+
+#if HAVE_CXSPARSE
+  retval = dmperm_internal (true, args(0), nargout);
+#else
+  error ("sprank: not available in this version of Octave");
+#endif
+
+  return retval;
+}
+
+/* 
+
+%!error(sprank(1,2));
+%!assert(sprank(speye(20)), 20)
+%!assert(sprank([1,0,2,0;2,0,4,0]),2)
+
+*/
 /*
 ;;; Local Variables: ***
 ;;; mode: C++ ***