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;
 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
 
 @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.
 
 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.
  *
 /** 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
  *  @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();
        
        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)
        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];
        }
                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:
        // 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);
        
        // assemble the solution matrix:
        matrix sol(n,p);


@@ -1097,20 +1056,23 @@ matrix matrix::solve(const matrix & vars,
        return sol;
 }
 
        return sol;
 }
 

-

 /** Compute the rank of this matrix. */
 unsigned matrix::rank() const
 /** 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.
 {
        // Method:
        // Transform this matrix into upper echelon form and then count the
        // number of non-zero rows.

-

        GINAC_ASSERT(row*col==m.capacity());
 
        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;
        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--) {
 
        unsigned r = row*col;  // index of last non-zero element
        while (r--) {


@@ -1242,6 +1204,52 @@ ex matrix::determinant_minor() const
        return det;
 }
 
        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
 
 /** 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;
        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;
        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);
        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)
 { return m.rank(); }

+inline unsigned rank(const matrix & m, unsigned solve_algo)
+{ return m.rank(solve_algo); }

 
 // utility functions
 
 
 // utility functions