X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=3cc20039e9f66e103ee4c1d19bd1fa36806e0aaf;hb=22abfbe8c78e339188096a5bf749a7c2d4f0a368;hp=4828bb9392db1e87666e382a77e5c4a5368f0cfb;hpb=b301f03a61bc9f72b27940ca7fe1f8d0b343a4e2;p=ginac.git diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 4828bb93..3cc20039 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-2005 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 @@ -37,22 +37,25 @@ #include "symbol.h" #include "operators.h" #include "normal.h" -#include "print.h" #include "archive.h" #include "utils.h" namespace GiNaC { -GINAC_IMPLEMENT_REGISTERED_CLASS(matrix, basic) +GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(matrix, basic, + print_func(&matrix::do_print). + print_func(&matrix::do_print_latex). + print_func(&matrix::do_print_tree). + print_func(&matrix::do_print_python_repr)) ////////// // default constructor ////////// /** Default ctor. Initializes to 1 x 1-dimensional zero-matrix. */ -matrix::matrix() : inherited(TINFO_matrix), row(1), col(1) +matrix::matrix() : inherited(TINFO_matrix), row(1), col(1), m(1, _ex0) { - m.push_back(_ex0); + setflag(status_flags::not_shareable); } ////////// @@ -66,25 +69,28 @@ matrix::matrix() : inherited(TINFO_matrix), row(1), col(1) * @param r number of rows * @param c number of cols */ matrix::matrix(unsigned r, unsigned c) - : inherited(TINFO_matrix), row(r), col(c) + : inherited(TINFO_matrix), row(r), col(c), m(r*c, _ex0) { - m.resize(r*c, _ex0); + setflag(status_flags::not_shareable); } // protected /** 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(TINFO_matrix), row(r), col(c), m(m2) +{ + setflag(status_flags::not_shareable); +} /** Construct matrix from (flat) list of elements. If the list has fewer * elements than the matrix, the remaining matrix elements are set to zero. * 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) + : inherited(TINFO_matrix), row(r), col(c), m(r*c, _ex0) { - m.resize(r*c, _ex0); + setflag(status_flags::not_shareable); size_t i = 0; for (lst::const_iterator it = l.begin(); it != l.end(); ++it, ++i) { @@ -102,6 +108,8 @@ matrix::matrix(unsigned r, unsigned c, const lst & l) matrix::matrix(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) { + setflag(status_flags::not_shareable); + if (!(n.find_unsigned("row", row)) || !(n.find_unsigned("col", col))) throw (std::runtime_error("unknown matrix dimensions in archive")); m.reserve(row * col); @@ -134,54 +142,41 @@ DEFAULT_UNARCHIVE(matrix) // public -void matrix::print(const print_context & c, unsigned level) const +void matrix::print_elements(const print_context & c, const char *row_start, const char *row_end, const char *row_sep, const char *col_sep) const { - if (is_a(c)) { - - inherited::print(c, level); - - } else { - - if (is_a(c)) - c.s << class_name() << '('; - - if (is_a(c)) - c.s << "\\left(\\begin{array}{" << std::string(col,'c') << "}"; - else - c.s << "["; - - for (unsigned ro=0; ro(c)) - c.s << "["; - for (unsigned co=0; co(c)) - c.s << "&"; - else - c.s << ","; - } else { - if (!is_a(c)) - c.s << "]"; - } - } - if (ro(c)) - c.s << "\\\\"; - else - c.s << ","; - } + for (unsigned ro=0; ro(c)) - c.s << "\\end{array}\\right)"; - else - c.s << "]"; +void matrix::do_print(const print_context & c, unsigned level) const +{ + c.s << "["; + print_elements(c, "[", "]", ",", ","); + c.s << "]"; +} - if (is_a(c)) - c.s << ')'; +void matrix::do_print_latex(const print_latex & c, unsigned level) const +{ + c.s << "\\left(\\begin{array}{" << std::string(col,'c') << "}"; + print_elements(c, "", "", "\\\\", "&"); + c.s << "\\end{array}\\right)"; +} - } +void matrix::do_print_python_repr(const print_python_repr & c, unsigned level) const +{ + c.s << class_name() << '('; + print_elements(c, "[", "]", ",", ","); + c.s << ')'; } /** nops is defined to be rows x columns. */ @@ -239,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 @@ -624,12 +647,12 @@ matrix matrix::pow(const ex & expn) const // that this is not entirely optimal but close to optimal and // "better" algorithms are much harder to implement. (See Knuth, // TAoCP2, section "Evaluation of Powers" for a good discussion.) - while (b!=_num1) { + while (b!=*_num1_p) { if (b.is_odd()) { C = C.mul(A); --b; } - b /= _num2; // still integer. + b /= *_num2_p; // still integer. A = A.mul(A); } return A.mul(C); @@ -707,7 +730,7 @@ ex matrix::determinant(unsigned algo) const unsigned sparse_count = 0; // counts non-zero elements exvector::const_iterator r = m.begin(), rend = m.end(); while (r != rend) { - lst srl; // symbol replacement list + exmap srl; // symbol replacement list ex rtest = r->to_rational(srl); if (!rtest.is_zero()) ++sparse_count; @@ -741,7 +764,7 @@ ex matrix::determinant(unsigned algo) const else return m[0].expand(); } - + // Compute the determinant switch(algo) { case determinant_algo::gauss: { @@ -837,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(); @@ -855,7 +878,7 @@ ex matrix::trace() const * @return characteristic polynomial as new expression * @exception logic_error (matrix not square) * @see matrix::determinant() */ -ex matrix::charpoly(const symbol & lambda) const +ex matrix::charpoly(const ex & lambda) const { if (row != col) throw (std::logic_error("matrix::charpoly(): matrix not square")); @@ -875,13 +898,13 @@ ex matrix::charpoly(const symbol & lambda) const matrix B(*this); ex c = B.trace(); - ex poly = power(lambda,row)-c*power(lambda,row-1); + ex poly = power(lambda, row) - c*power(lambda, row-1); for (unsigned i=1; imul(B); c = B.trace() / ex(i+1); - poly -= c*power(lambda,row-i-1); + poly -= c*power(lambda, row-i-1); } if (row%2) return -poly; @@ -943,14 +966,15 @@ matrix matrix::inverse() const * * @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 * @exception logic_error (incompatible matrices) * @exception invalid_argument (1st argument must be matrix of symbols) * @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(); @@ -1043,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 @@ -1158,7 +1205,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(); } @@ -1184,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) { @@ -1216,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; } @@ -1238,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) { @@ -1263,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; } @@ -1325,7 +1382,7 @@ int matrix::fraction_free_elimination(const bool det) // makes things more complicated than they need to be. matrix tmp_n(*this); matrix tmp_d(m,n); // for denominators, if needed - lst srl; // symbol replacement list + exmap srl; // symbol replacement list exvector::const_iterator cit = this->m.begin(), citend = this->m.end(); exvector::iterator tmp_n_it = tmp_n.m.begin(), tmp_d_it = tmp_d.m.begin(); while (cit != citend) { @@ -1336,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();