X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=6464ed8d52cd90731511ddf04e976d915fb734b9;hp=7b0b34889f396cfb66feda3f0011e7fb45378393;hb=26840ecc367d9f9dd4e3758233be1bdd3563cd8a;hpb=48619ed77871a6bcae23df460f426fc34698cd1e;ds=sidebyside 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