X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=f62e2d409244d8a50d44890b8eae74f8ac5e4eb1;hp=2c191ff4f69c04d993b2724b16f13e13d445156e;hb=1c85bdb25fe8fed52caf39e5ec35d6d14e6d837c;hpb=3db194853d3dd2719231b83e4956cf6bcbe53627 diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 2c191ff4..f62e2d40 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -3,7 +3,7 @@ * Implementation of symbolic matrices */ /* - * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -17,7 +17,7 @@ * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software - * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ #include @@ -45,7 +45,7 @@ namespace GiNaC { GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(matrix, basic, print_func(&matrix::do_print). print_func(&matrix::do_print_latex). - print_func(&basic::do_print_tree). + print_func(&matrix::do_print_tree). print_func(&matrix::do_print_python_repr)) ////////// @@ -53,7 +53,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(matrix, basic, ////////// /** Default ctor. Initializes to 1 x 1-dimensional zero-matrix. */ -matrix::matrix() : inherited(TINFO_matrix), row(1), col(1), m(1, _ex0) +matrix::matrix() : inherited(&matrix::tinfo_static), row(1), col(1), m(1, _ex0) { setflag(status_flags::not_shareable); } @@ -69,7 +69,7 @@ matrix::matrix() : inherited(TINFO_matrix), row(1), col(1), m(1, _ex0) * @param r number of rows * @param c number of cols */ matrix::matrix(unsigned r, unsigned c) - : inherited(TINFO_matrix), row(r), col(c), m(r*c, _ex0) + : inherited(&matrix::tinfo_static), row(r), col(c), m(r*c, _ex0) { setflag(status_flags::not_shareable); } @@ -78,7 +78,7 @@ matrix::matrix(unsigned r, unsigned c) /** Ctor from representation, for internal use only. */ matrix::matrix(unsigned r, unsigned c, const exvector & m2) - : inherited(TINFO_matrix), row(r), col(c), m(m2) + : inherited(&matrix::tinfo_static), row(r), col(c), m(m2) { setflag(status_flags::not_shareable); } @@ -88,7 +88,7 @@ matrix::matrix(unsigned r, unsigned c, const exvector & m2) * If the list has more elements than the matrix, the excessive elements are * thrown away. */ matrix::matrix(unsigned r, unsigned c, const lst & l) - : inherited(TINFO_matrix), row(r), col(c), m(r*c, _ex0) + : inherited(&matrix::tinfo_static), row(r), col(c), m(r*c, _ex0) { setflag(status_flags::not_shareable); @@ -234,6 +234,34 @@ ex matrix::subs(const exmap & mp, unsigned options) const return matrix(row, col, m2).subs_one_level(mp, options); } +/** Complex conjugate every matrix entry. */ +ex matrix::conjugate() const +{ + exvector * ev = 0; + for (exvector::const_iterator i=m.begin(); i!=m.end(); ++i) { + ex x = i->conjugate(); + if (ev) { + ev->push_back(x); + continue; + } + if (are_ex_trivially_equal(x, *i)) { + continue; + } + ev = new exvector; + ev->reserve(m.size()); + for (exvector::const_iterator j=m.begin(); j!=i; ++j) { + ev->push_back(*j); + } + ev->push_back(x); + } + if (ev) { + ex result = matrix(row, col, *ev); + delete ev; + return result; + } + return *this; +} + // protected int matrix::compare_same_type(const basic & other) const @@ -552,10 +580,11 @@ matrix matrix::mul(const matrix & other) const for (unsigned r1=0; r1rows(); ++r1) { for (unsigned c=0; ccols(); ++c) { + // Quick test: can we shortcut? if (m[r1*col+c].is_zero()) continue; for (unsigned r2=0; r2info(info_flags::numeric)) + numeric_flag = false; exmap srl; // symbol replacement list ex rtest = r->to_rational(srl); if (!rtest.is_zero()) ++sparse_count; - if (!rtest.info(info_flags::numeric)) - numeric_flag = false; if (!rtest.info(info_flags::crational_polynomial) && rtest.info(info_flags::rational_function)) normal_flag = true; @@ -736,7 +765,7 @@ ex matrix::determinant(unsigned algo) const else return m[0].expand(); } - + // Compute the determinant switch(algo) { case determinant_algo::gauss: { @@ -832,7 +861,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(); @@ -945,8 +974,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(); @@ -1039,6 +1068,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 @@ -1154,7 +1206,7 @@ ex matrix::determinant_minor() const Pkey[j] = Pkey[j-1]+1; } while(fc); // next column, so change the role of A and B: - A = B; + A.swap(B); B.clear(); } @@ -1180,8 +1232,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) { @@ -1212,7 +1264,12 @@ int matrix::gauss_elimination(const bool det) ++r0; } } - + // clear remaining rows + for (unsigned r=r0+1; rm[r*n+c] = _ex0; + } + return sign; } @@ -1234,8 +1291,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) { @@ -1259,7 +1316,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; } @@ -1332,8 +1394,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(); @@ -1531,4 +1599,52 @@ ex symbolic_matrix(unsigned r, unsigned c, const std::string & base_name, const return M; } +ex reduced_matrix(const matrix& m, unsigned r, unsigned c) +{ + if (r+1>m.rows() || c+1>m.cols() || m.cols()<2 || m.rows()<2) + throw std::runtime_error("minor_matrix(): index out of bounds"); + + const unsigned rows = m.rows()-1; + const unsigned cols = m.cols()-1; + matrix &M = *new matrix(rows, cols); + M.setflag(status_flags::dynallocated | status_flags::evaluated); + + unsigned ro = 0; + unsigned ro2 = 0; + while (ro2m.rows() || c+nc>m.cols()) + throw std::runtime_error("sub_matrix(): index out of bounds"); + + matrix &M = *new matrix(nr, nc); + M.setflag(status_flags::dynallocated | status_flags::evaluated); + + for (unsigned ro=0; ro