sprank function

[David Bateman]

Please consider the attached patch that adds an sprank function to
octave that is compatible with matlab. Note that sprank returns the
structural rank of the matrix, rather than the numerical rank, and only
depends of the structure of the sparse matrix that is used, which is the
equivalent of calculating the rank without taking into account any
floating point errors. Therefore sprank(A) >= rank(A)....

Regards
David

2006-10-20  David Bateman  <address@hidden>

        * 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.

Index: src/DLD-FUNCTIONS/spqr.cc
===================================================================
RCS file: /cvs/octave/src/DLD-FUNCTIONS/spqr.cc,v
retrieving revision 1.7
diff -c -r1.7 spqr.cc
*** src/DLD-FUNCTIONS/spqr.cc   19 May 2006 05:32:19 -0000      1.7
--- src/DLD-FUNCTIONS/spqr.cc   19 Oct 2006 21:19:49 -0000
***************
*** 230,268 ****
    return ret;
  }
  
! DEFUN_DLD (dmperm, args, nargout,
!   "-*- texinfo -*-\n\
! @deftypefn {Loadable Function} address@hidden =} dmperm (@var{s})\n\
! @deftypefnx {Loadable Function} address@hidden, @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 @address@hidden (@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 @address@hidden (@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 @address@hidden == @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 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;
--- 230,240 ----
    return ret;
  }
  
! #if HAVE_CXSPARSE
! static octave_value_list
! dmperm_internal (bool rank, const octave_value arg, int nargout)
  {
    octave_value_list retval;
    octave_idx_type nr = arg.rows ();
    octave_idx_type nc = arg.columns ();
    SparseMatrix m;
***************
*** 290,303 ****
  
    if (!error_state)
      {
!       if (nargout <= 1)
        {
  #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);
          CXSPARSE_NAME (_free) (jmatch);
        }
        else
--- 262,284 ----
  
    if (!error_state)
      {
!       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
!         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,312 ****
--- 288,294 ----
  #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,328 ****
--- 305,347 ----
          CXSPARSE_NAME (_dfree) (dm);
        }
      }
+   return retval;
+ }
+ #endif
+ 
+ DEFUN_DLD (dmperm, args, nargout,
+   "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} address@hidden =} dmperm (@var{s})\n\
+ @deftypefnx {Loadable Function} address@hidden, @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 @address@hidden (@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 @address@hidden (@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 @address@hidden == @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,336 ****
    return retval;
  }
  
! /*
  
  %!test
  %! n=20;
--- 349,355 ----
    return retval;
  }
  
! /* 
  
  %!test
  %! n=20;
***************
*** 347,352 ****
--- 366,410 ----
  
  */
  
+ DEFUN_DLD (sprank, args, nargout,
+   "-*- texinfo -*-\n\
+ @deftypefn {Loadable Function} address@hidden =} 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++ ***