From: Richard Kreckel Date: Tue, 15 Jun 2004 01:48:54 +0000 (+0000) Subject: - Added method matrix::rank(). X-Git-Tag: release_1-3-0~64 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=b5483c25af5d84095ca647390914e3fb63ad96a0 - Added method matrix::rank(). - Fixes to matrix elimination routines, needed for rank computations. --- diff --git a/check/exam_matrices.cpp b/check/exam_matrices.cpp index 67c52be8..449b6d38 100644 --- a/check/exam_matrices.cpp +++ b/check/exam_matrices.cpp @@ -241,6 +241,46 @@ static unsigned matrix_evalm() return result; } +static unsigned matrix_rank() +{ + unsigned result = 0; + symbol x("x"), y("y"); + matrix m(3,3); + + // the zero matrix always has rank 0 + if (m.rank() != 0) { + clog << "The rank of " << m << " was not computed correctly." << endl; + ++result; + } + + // a trivial rank one example + m = 1, 0, 0, + 2, 0, 0, + 3, 0, 0; + if (m.rank() != 1) { + clog << "The rank of " << m << " was not computed correctly." << endl; + ++result; + } + + // an example from Maple's help with rank two + m = x, 1, 0, + 0, 0, 1, + x*y, y, 1; + if (m.rank() != 2) { + clog << "The rank of " << m << " was not computed correctly." << endl; + ++result; + } + + // the 3x3 unit matrix has rank 3 + m = ex_to(unit_matrix(3,3)); + if (m.rank() != 3) { + clog << "The rank of " << m << " was not computed correctly." << endl; + ++result; + } + + return result; +} + static unsigned matrix_misc() { unsigned result = 0; @@ -305,6 +345,7 @@ unsigned exam_matrices() result += matrix_invert3(); cout << '.' << flush; result += matrix_solve2(); cout << '.' << flush; result += matrix_evalm(); cout << "." << flush; + result += matrix_rank(); cout << "." << flush; result += matrix_misc(); cout << '.' << flush; if (!result) { diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index 3f981831..b58d631b 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -2013,15 +2013,17 @@ more information about using matrices with indices, and about indices in general. The @code{matrix} class provides a couple of additional methods for -computing determinants, traces, and characteristic polynomials: +computing determinants, traces, characteristic polynomials and ranks: @cindex @code{determinant()} @cindex @code{trace()} @cindex @code{charpoly()} +@cindex @code{rank()} @example ex matrix::determinant(unsigned algo=determinant_algo::automatic) const; ex matrix::trace() const; ex matrix::charpoly(const ex & lambda) const; +unsigned matrix::rank() const; @end example The @samp{algo} argument of @code{determinant()} allows to select diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 2c6e5bb4..e24f6447 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -764,7 +764,7 @@ ex matrix::determinant(unsigned algo) const else return m[0].expand(); } - + // Compute the determinant switch(algo) { case determinant_algo::gauss: { @@ -860,7 +860,7 @@ ex matrix::trace() const tr += m[r*col+r]; if (tr.info(info_flags::rational_function) && - !tr.info(info_flags::crational_polynomial)) + !tr.info(info_flags::crational_polynomial)) return tr.normal(); else return tr.expand(); @@ -973,8 +973,8 @@ matrix matrix::inverse() const * @exception runtime_error (inconsistent linear system) * @see solve_algo */ matrix matrix::solve(const matrix & vars, - const matrix & rhs, - unsigned algo) const + const matrix & rhs, + unsigned algo) const { const unsigned m = this->rows(); const unsigned n = this->cols(); @@ -1067,6 +1067,29 @@ matrix matrix::solve(const matrix & vars, } +/** Compute the rank of this matrix. */ +unsigned matrix::rank() 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(); + + unsigned r = row*col; // index of last non-zero element + while (r--) { + if (!to_eliminate.m[r].is_zero()) + return 1+r/col; + } + return 0; +} + + // protected /** Recursive determinant for small matrices having at least one symbolic @@ -1208,8 +1231,8 @@ int matrix::gauss_elimination(const bool det) int sign = 1; unsigned r0 = 0; - for (unsigned r1=0; (r1 0) sign = -sign; for (unsigned r2=r0+1; r2m[r2*n+r1].is_zero()) { + if (!this->m[r2*n+c0].is_zero()) { // yes, there is something to do in this row - ex piv = this->m[r2*n+r1] / this->m[r0*n+r1]; - for (unsigned c=r1+1; cm[r2*n+c0] / this->m[r0*n+c0]; + for (unsigned c=c0+1; cm[r2*n+c] -= piv * this->m[r0*n+c]; if (!this->m[r2*n+c].info(info_flags::numeric)) this->m[r2*n+c] = this->m[r2*n+c].normal(); } } // fill up left hand side with zeros - for (unsigned c=0; c<=r1; ++c) + for (unsigned c=r0; c<=c0; ++c) this->m[r2*n+c] = _ex0; } if (det) { @@ -1240,7 +1263,12 @@ int matrix::gauss_elimination(const bool det) ++r0; } } - + // clear remaining rows + for (unsigned r=r0+1; rm[r*n+c] = _ex0; + } + return sign; } @@ -1262,8 +1290,8 @@ int matrix::division_free_elimination(const bool det) int sign = 1; unsigned r0 = 0; - for (unsigned r1=0; (r10) sign = -sign; for (unsigned r2=r0+1; r2m[r2*n+c] = (this->m[r0*n+r1]*this->m[r2*n+c] - this->m[r2*n+r1]*this->m[r0*n+c]).expand(); + for (unsigned c=c0+1; cm[r2*n+c] = (this->m[r0*n+c0]*this->m[r2*n+c] - this->m[r2*n+c0]*this->m[r0*n+c]).expand(); // fill up left hand side with zeros - for (unsigned c=0; c<=r1; ++c) + for (unsigned c=r0; c<=c0; ++c) this->m[r2*n+c] = _ex0; } if (det) { @@ -1287,7 +1315,12 @@ int matrix::division_free_elimination(const bool det) ++r0; } } - + // clear remaining rows + for (unsigned r=r0+1; rm[r*n+c] = _ex0; + } + return sign; } @@ -1360,8 +1393,8 @@ int matrix::fraction_free_elimination(const bool det) } unsigned r0 = 0; - for (unsigned r1=0; (r10) { sign = -sign; // tmp_n's rows r0 and indx were swapped, do the same in tmp_d: - for (unsigned c=r1; cm.begin(), itend = this->m.end(); tmp_n_it = tmp_n.m.begin(); diff --git a/ginac/matrix.h b/ginac/matrix.h index ee306013..29af2a07 100644 --- a/ginac/matrix.h +++ b/ginac/matrix.h @@ -148,6 +148,7 @@ public: matrix inverse() const; matrix solve(const matrix & vars, const matrix & rhs, unsigned algo = solve_algo::automatic) const; + unsigned rank() const; protected: ex determinant_minor() const; int gauss_elimination(const bool det = false); @@ -203,6 +204,9 @@ inline ex charpoly(const matrix & m, const ex & lambda) inline matrix inverse(const matrix & m) { return m.inverse(); } +inline unsigned rank(const matrix & m) +{ return m.rank(); } + // utility functions /** Specialization of is_exactly_a(obj) for matrix objects. */ diff --git a/ginsh/ginsh.1.in b/ginsh/ginsh.1.in index 01ecbb59..1662ce29 100644 --- a/ginsh/ginsh.1.in +++ b/ginsh/ginsh.1.in @@ -350,6 +350,9 @@ detail here. Please refer to the GiNaC documentation. .BI quo( expression ", " expression ", " symbol ) \- quotient of polynomials .br +.BI rank( matrix ) +\- rank of a matrix +.br .BI rem( expression ", " expression ", " symbol ) \- remainder of polynomials .br diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy index c0458267..ed2fcda6 100644 --- a/ginsh/ginsh_parser.yy +++ b/ginsh/ginsh_parser.yy @@ -469,6 +469,12 @@ static ex f_quo(const exprseq &e) return quo(e[0], e[1], e[2]); } +static ex f_rank(const exprseq &e) +{ + CHECK_ARG(0, matrix, rank); + return ex_to(e[0]).rank(); +} + static ex f_rem(const exprseq &e) { return rem(e[0], e[1], e[2]); @@ -581,6 +587,7 @@ static const fcn_init builtin_fcns[] = { {"print_csrc", f_dummy, 0}, // for Tab-completion {"print_latex", f_dummy, 0}, // for Tab-completion {"quo", f_quo, 3}, + {"rank", f_rank, 1}, {"rem", f_rem, 3}, {"series", f_series, 3}, {"sprem", f_sprem, 3},