X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=bc58245f8f0e7cc8cf33134ce30b80e2cddd8590;hp=b9193a6cd1caa488eba441a902c15373803547a3;hb=36ffb97c48fbbd3acb2b1aa305872e00a4456dc4;hpb=ff3d5ccc8add5bfe76502479768c89d768df8384 diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index b9193a6c..bc58245f 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,52 @@ 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; +} + +ex matrix::real_part() const +{ + exvector v; + v.reserve(m.size()); + for (exvector::const_iterator i=m.begin(); i!=m.end(); ++i) + v.push_back(i->real_part()); + return matrix(row, col, v); +} + +ex matrix::imag_part() const +{ + exvector v; + v.reserve(m.size()); + for (exvector::const_iterator i=m.begin(); i!=m.end(); ++i) + v.push_back(i->imag_part()); + return matrix(row, col, v); +} + // protected int matrix::compare_same_type(const basic & other) const @@ -552,10 +598,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 +783,7 @@ ex matrix::determinant(unsigned algo) const else return m[0].expand(); } - + // Compute the determinant switch(algo) { case determinant_algo::gauss: { @@ -832,7 +879,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(); @@ -938,14 +985,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(); @@ -1038,6 +1086,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 @@ -1153,7 +1224,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(); } @@ -1179,8 +1250,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) { @@ -1211,7 +1282,12 @@ int matrix::gauss_elimination(const bool det) ++r0; } } - + // clear remaining rows + for (unsigned r=r0+1; rm[r*n+c] = _ex0; + } + return sign; } @@ -1233,8 +1309,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) { @@ -1258,7 +1334,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; } @@ -1320,7 +1401,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) { @@ -1331,8 +1412,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(); @@ -1440,6 +1527,16 @@ int matrix::pivot(unsigned ro, unsigned co, bool symbolic) return k; } +/** Function to check that all elements of the matrix are zero. + */ +bool matrix::is_zero_matrix() const +{ + for (exvector::const_iterator i=m.begin(); i!=m.end(); ++i) + if(!(i->is_zero())) + return false; + return true; +} + ex lst_to_matrix(const lst & l) { lst::const_iterator itr, itc; @@ -1530,4 +1627,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