From: Vitaly Magerya Date: Wed, 13 Jun 2018 21:07:45 +0000 (+0200) Subject: Add optional matrix::rank() algorighm selection. X-Git-Tag: release_1-7-5~8 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=26840ecc367d9f9dd4e3758233be1bdd3563cd8a;ds=sidebyside 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 3ec34c35..3742567e 100644 --- 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 7b0b3488..6464ed8d 100644 --- 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 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 +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 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 cf7dd039..44351d65 100644 --- 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 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