From bc01339d2eb85f4efa433bc10b97c4343d477afb Mon Sep 17 00:00:00 2001 From: Vitaly Magerya Date: Thu, 7 Jun 2018 17:14:12 +0200 Subject: [PATCH] Add 'algo' parameter to matrix::rank() Also share the automatic elimination algorithm selection between matrix::solve() and matrix::rank(). --- ginac/matrix.cpp | 108 +++++++++++++++++++++++++++++-------------------------- ginac/matrix.h | 4 +++ 2 files changed, 62 insertions(+), 50 deletions(-) diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 7b0b3488..3a01da2d 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 -- 2.13.6