Add optional matrix::rank() algorighm selection.

Before, matrix::rank() was hardcoded to use Bareiss elimination,
which for some matrices is way too slow (and I had to work around
this limitation in my own code). I propose adding an argument that
will allow users to select which elimination scheme to use. The
argument is the same as matrix::solve() takes (so, solve_algo::*).

I've tried to mimic how a similar argument was added to
matrix::inverse() by keeping the current set of rank() functions
and defining separate ones with the 'algo' argument, with the
old ones calling the new ones with solve_algo::automatic. This
is instead of using default argument values.

I've also put the code that makes the automatic selection of
the elimination method into a separate function so that both
matrix::solve() and matrix::rank() could share it. There's a bit
of similar code in matrix::determinant(), but since the self
of determinant algorithms is different from the elimination
algorithms, it's not clear if these two could be unified.

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 3ec34c35032d25a78942fa615849b37e1808c353..3742567eca0ba3d331409d7060ea6caafeb0f6d9 100644 (file)
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -2172,15 +2172,15 @@ computing determinants, traces, characteristic polynomials and ranks:
 ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
 ex matrix::trace() const;
 ex matrix::charpoly(const ex & lambda) const;
-unsigned matrix::rank() const;
+unsigned matrix::rank(unsigned algo=solve_algo::automatic) const;
 @end example
 
-The optional @samp{algo} argument of @code{determinant()} allows to
-select between different algorithms for calculating the determinant.
-The asymptotic speed (as parametrized by the matrix size) can greatly
-differ between those algorithms, depending on the nature of the
-matrix' entries.  The possible values are defined in the
-@file{flags.h} header file.  By default, GiNaC uses a heuristic to
+The optional @samp{algo} argument of @code{determinant()} and @code{rank()}
+functions allows to select between different algorithms for calculating the
+determinant and rank respectively. The asymptotic speed (as parametrized
+by the matrix size) can greatly differ between those algorithms, depending
+on the nature of the matrix' entries. The possible values are defined in
+the @file{flags.h} header file. By default, GiNaC uses a heuristic to
 automatically select an algorithm that is likely (but not guaranteed)
 to give the result most quickly.
 

diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp
index 7b0b34889f396cfb66feda3f0011e7fb45378393..6464ed8d52cd90731511ddf04e976d915fb734b9 100644 (file)
--- a/ginac/matrix.cpp
+++ b/ginac/matrix.cpp
@@ -986,7 +986,7 @@ matrix matrix::inverse(unsigned algo) const
 /** Solve a linear system consisting of a m x n matrix and a m x p right hand
  *  side by applying an elimination scheme to the augmented matrix.
  *
- *  @param vars n x p matrix, all elements must be symbols 
+ *  @param vars n x p matrix, all elements must be symbols
  *  @param rhs m x p matrix
  *  @param algo selects the solving algorithm
  *  @return n x p solution matrix
@@ -1002,7 +1002,7 @@ matrix matrix::solve(const matrix & vars,
        const unsigned n = this->cols();
        const unsigned p = rhs.cols();
        
-       // syntax checks    
+       // syntax checks
        if ((rhs.rows() != m) || (vars.rows() != n) || (vars.cols() != p))
                throw (std::logic_error("matrix::solve(): incompatible matrices"));
        for (unsigned ro=0; ro<n; ++ro)
@@ -1018,50 +1018,9 @@ matrix matrix::solve(const matrix & vars,
                for (unsigned c=0; c<p; ++c)
                        aug.m[r*(n+p)+c+n] = rhs.m[r*p+c];
        }
-       
-       // Gather some statistical information about the augmented matrix:
-       bool numeric_flag = true;
-       for (auto & r : aug.m) {
-               if (!r.info(info_flags::numeric)) {
-                       numeric_flag = false;
-                       break;
-               }
-       }
-       
-       // Here is the heuristics in case this routine has to decide:
-       if (algo == solve_algo::automatic) {
-               // Bareiss (fraction-free) elimination is generally a good guess:
-               algo = solve_algo::bareiss;
-               // For m<3, Bareiss elimination is equivalent to division free
-               // elimination but has more logistic overhead
-               if (m<3)
-                       algo = solve_algo::divfree;
-               // This overrides any prior decisions.
-               if (numeric_flag)
-                       algo = solve_algo::gauss;
-       }
-       
+
        // Eliminate the augmented matrix:
-       std::vector<unsigned> colid(aug.cols());
-       for (unsigned c = 0; c < aug.cols(); c++) {
-               colid[c] = c;
-       }
-       switch(algo) {
-               case solve_algo::gauss:
-                       aug.gauss_elimination();
-                       break;
-               case solve_algo::divfree:
-                       aug.division_free_elimination();
-                       break;
-               case solve_algo::bareiss:
-                       aug.fraction_free_elimination();
-                       break;
-               case solve_algo::markowitz:
-                       colid = aug.markowitz_elimination(n);
-                       break;
-               default:
-                       throw std::invalid_argument("matrix::solve(): 'algo' is not one of the solve_algo enum");
-       }
+       auto colid = aug.echelon_form(algo, n);
        
        // assemble the solution matrix:
        matrix sol(n,p);
@@ -1097,20 +1056,23 @@ matrix matrix::solve(const matrix & vars,
        return sol;
 }
 
-
 /** Compute the rank of this matrix. */
 unsigned matrix::rank() const
+{
+       return rank(solve_algo::automatic);
+}
+
+/** Compute the rank of this matrix using the given algorithm,
+ *  which should be a member of enum solve_algo. */
+unsigned matrix::rank(unsigned solve_algo) const
 {
        // Method:
        // Transform this matrix into upper echelon form and then count the
        // number of non-zero rows.
-
        GINAC_ASSERT(row*col==m.capacity());
 
-       // Actually, any elimination scheme will do since we are only
-       // interested in the echelon matrix' zeros.
        matrix to_eliminate = *this;
-       to_eliminate.fraction_free_elimination();
+       to_eliminate.echelon_form(solve_algo, col);
 
        unsigned r = row*col;  // index of last non-zero element
        while (r--) {
@@ -1242,6 +1204,52 @@ ex matrix::determinant_minor() const
        return det;
 }
 
+std::vector<unsigned>
+matrix::echelon_form(unsigned algo, int n)
+{
+       // Here is the heuristics in case this routine has to decide:
+       if (algo == solve_algo::automatic) {
+               // Gather some statistical information about the augmented matrix:
+               bool numeric_flag = true;
+               for (auto & r : m) {
+                       if (!r.info(info_flags::numeric)) {
+                               numeric_flag = false;
+                               break;
+                       }
+               }
+               // Bareiss (fraction-free) elimination is generally a good guess:
+               algo = solve_algo::bareiss;
+               // For row<3, Bareiss elimination is equivalent to division free
+               // elimination but has more logistic overhead
+               if (row<3)
+                       algo = solve_algo::divfree;
+               // This overrides any prior decisions.
+               if (numeric_flag)
+                       algo = solve_algo::gauss;
+       }
+       // Eliminate the augmented matrix:
+       std::vector<unsigned> colid(col);
+       for (unsigned c = 0; c < col; c++) {
+               colid[c] = c;
+       }
+       switch(algo) {
+               case solve_algo::gauss:
+                       gauss_elimination();
+                       break;
+               case solve_algo::divfree:
+                       division_free_elimination();
+                       break;
+               case solve_algo::bareiss:
+                       fraction_free_elimination();
+                       break;
+               case solve_algo::markowitz:
+                       colid = markowitz_elimination(n);
+                       break;
+               default:
+                       throw std::invalid_argument("matrix::echelon_form(): 'algo' is not one of the solve_algo enum");
+       }
+       return colid;
+}
 
 /** Perform the steps of an ordinary Gaussian elimination to bring the m x n
  *  matrix into an upper echelon form.  The algorithm is ok for matrices

diff --git a/ginac/matrix.h b/ginac/matrix.h
index cf7dd03954a2e73331d26a113499380584253769..44351d65885f6f11b755fa08b4476302c05fc087 100644 (file)
--- a/ginac/matrix.h
+++ b/ginac/matrix.h
@@ -153,9 +153,11 @@ public:
        matrix solve(const matrix & vars, const matrix & rhs,
                     unsigned algo = solve_algo::automatic) const;
        unsigned rank() const;
+       unsigned rank(unsigned solve_algo) const;
        bool is_zero_matrix() const;
 protected:
        ex determinant_minor() const;
+       std::vector<unsigned> echelon_form(unsigned algo, int n);
        int gauss_elimination(const bool det = false);
        int division_free_elimination(const bool det = false);
        int fraction_free_elimination(const bool det = false);
@@ -219,6 +221,8 @@ inline matrix inverse(const matrix & m, unsigned algo)
 
 inline unsigned rank(const matrix & m)
 { return m.rank(); }
+inline unsigned rank(const matrix & m, unsigned solve_algo)
+{ return m.rank(solve_algo); }
 
 // utility functions