X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=bf7f8347c9a1f7b31645595d408ce602cf77081d;hp=b9193a6cd1caa488eba441a902c15373803547a3;hb=d0ff428fb5b7bb565a0aea69e05e5705d840c16d;hpb=ff3d5ccc8add5bfe76502479768c89d768df8384 diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index b9193a6c..bf7f8347 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-2018 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,16 +17,9 @@ * * 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 -#include -#include -#include -#include -#include - #include "matrix.h" #include "numeric.h" #include "lst.h" @@ -40,12 +33,19 @@ #include "archive.h" #include "utils.h" +#include +#include +#include +#include +#include +#include + 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() : row(1), col(1), m(1, _ex0) { setflag(status_flags::not_shareable); } @@ -68,17 +68,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) -{ - 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) +matrix::matrix(unsigned r, unsigned c) : row(r), col(c), m(r*c, _ex0) { setflag(status_flags::not_shareable); } @@ -88,54 +78,89 @@ 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) + : row(r), col(c), m(r*c, _ex0) { setflag(status_flags::not_shareable); size_t i = 0; - for (lst::const_iterator it = l.begin(); it != l.end(); ++it, ++i) { + for (auto & it : l) { size_t x = i % c; size_t y = i / c; if (y >= r) break; // matrix smaller than list: throw away excessive elements - m[y*c+x] = *it; + m[y*c+x] = it; + ++i; + } +} + +/** Construct a matrix from an 2 dimensional initializer list. + * Throws an exception if some row has a different length than all the others. + */ +matrix::matrix(std::initializer_list> l) + : row(l.size()), col(l.begin()->size()) +{ + setflag(status_flags::not_shareable); + + m.reserve(row*col); + for (const auto & r : l) { + unsigned c = 0; + for (const auto & e : r) { + m.push_back(e); + ++c; + } + if (c != col) + throw std::invalid_argument("matrix::matrix{{}}: wrong dimension"); } } +// protected + +/** Ctor from representation, for internal use only. */ +matrix::matrix(unsigned r, unsigned c, const exvector & m2) + : row(r), col(c), m(m2) +{ + setflag(status_flags::not_shareable); +} +matrix::matrix(unsigned r, unsigned c, exvector && m2) + : row(r), col(c), m(std::move(m2)) +{ + setflag(status_flags::not_shareable); +} + ////////// // archiving ////////// -matrix::matrix(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) +void matrix::read_archive(const archive_node &n, lst &sym_lst) { - setflag(status_flags::not_shareable); + inherited::read_archive(n, sym_lst); if (!(n.find_unsigned("row", row)) || !(n.find_unsigned("col", col))) throw (std::runtime_error("unknown matrix dimensions in archive")); m.reserve(row * col); - for (unsigned int i=0; true; i++) { + // XXX: default ctor inserts a zero element, we need to erase it here. + m.pop_back(); + auto first = n.find_first("m"); + auto last = n.find_last("m"); + ++last; + for (auto i=first; i != last; ++i) { ex e; - if (n.find_ex("m", e, sym_lst, i)) - m.push_back(e); - else - break; + n.find_ex_by_loc(i, e, sym_lst); + m.push_back(e); } } +GINAC_BIND_UNARCHIVER(matrix); void matrix::archive(archive_node &n) const { inherited::archive(n); n.add_unsigned("row", row); n.add_unsigned("col", col); - exvector::const_iterator i = m.begin(), iend = m.end(); - while (i != iend) { - n.add_ex("m", *i); - ++i; + for (auto & i : m) { + n.add_ex("m", i); } } -DEFAULT_UNARCHIVE(matrix) - ////////// // functions overriding virtual functions from base classes ////////// @@ -202,28 +227,6 @@ ex & matrix::let_op(size_t i) return m[i]; } -/** Evaluate matrix entry by entry. */ -ex matrix::eval(int level) const -{ - // check if we have to do anything at all - if ((level==1)&&(flags & status_flags::evaluated)) - return *this; - - // emergency break - if (level == -max_recursion_level) - throw (std::runtime_error("matrix::eval(): recursion limit exceeded")); - - // eval() entry by entry - exvector m2(row*col); - --level; - for (unsigned r=0; rsetflag(status_flags::dynallocated | - status_flags::evaluated); -} - ex matrix::subs(const exmap & mp, unsigned options) const { exvector m2(row * col); @@ -231,7 +234,51 @@ ex matrix::subs(const exmap & mp, unsigned options) const for (unsigned c=0; c ev(nullptr); + for (auto 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.reset(new exvector); + ev->reserve(m.size()); + for (auto j=m.begin(); j!=i; ++j) { + ev->push_back(*j); + } + ev->push_back(x); + } + if (ev) { + return matrix(row, col, std::move(*ev)); + } + return *this; +} + +ex matrix::real_part() const +{ + exvector v; + v.reserve(m.size()); + for (auto & i : m) + v.push_back(i.real_part()); + return matrix(row, col, std::move(v)); +} + +ex matrix::imag_part() const +{ + exvector v; + v.reserve(m.size()); + for (auto & i : m) + v.push_back(i.imag_part()); + return matrix(row, col, std::move(v)); } // protected @@ -513,12 +560,11 @@ matrix matrix::add(const matrix & other) const throw std::logic_error("matrix::add(): incompatible matrices"); exvector sum(this->m); - exvector::iterator i = sum.begin(), end = sum.end(); - exvector::const_iterator ci = other.m.begin(); - while (i != end) - *i++ += *ci++; + auto ci = other.m.begin(); + for (auto & i : sum) + i += *ci++; - return matrix(row,col,sum); + return matrix(row, col, std::move(sum)); } @@ -531,12 +577,11 @@ matrix matrix::sub(const matrix & other) const throw std::logic_error("matrix::sub(): incompatible matrices"); exvector dif(this->m); - exvector::iterator i = dif.begin(), end = dif.end(); - exvector::const_iterator ci = other.m.begin(); - while (i != end) - *i++ -= *ci++; + auto ci = other.m.begin(); + for (auto & i : dif) + i -= *ci++; - return matrix(row,col,dif); + return matrix(row, col, std::move(dif)); } @@ -552,13 +597,14 @@ 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; r2rows(); ++c) trans[r*this->rows()+c] = m[c*this->cols()+r]; - return matrix(this->cols(),this->rows(),trans); + return matrix(this->cols(), this->rows(), std::move(trans)); } /** Determinant of square matrix. This routine doesn't actually calculate the @@ -700,18 +746,16 @@ ex matrix::determinant(unsigned algo) const bool numeric_flag = true; bool normal_flag = false; 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 - ex rtest = r->to_rational(srl); + for (auto r : m) { + if (!r.info(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)) + rtest.info(info_flags::rational_function)) normal_flag = true; - ++r; } // Here is the heuristics in case this routine has to decide: @@ -736,7 +780,7 @@ ex matrix::determinant(unsigned algo) const else return m[0].expand(); } - + // Compute the determinant switch(algo) { case determinant_algo::gauss: { @@ -794,23 +838,22 @@ ex matrix::determinant(unsigned algo) const } std::sort(c_zeros.begin(),c_zeros.end()); std::vector pre_sort; - for (std::vector::const_iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i) - pre_sort.push_back(i->second); + for (auto & i : c_zeros) + pre_sort.push_back(i.second); std::vector pre_sort_test(pre_sort); // permutation_sign() modifies the vector so we make a copy here int sign = permutation_sign(pre_sort_test.begin(), pre_sort_test.end()); exvector result(row*col); // represents sorted matrix unsigned c = 0; - for (std::vector::const_iterator i=pre_sort.begin(); - i!=pre_sort.end(); - ++i,++c) { + for (auto & it : pre_sort) { for (unsigned r=0; rinfo(info_flags::numeric)) + for (auto & r : m) { + if (!r.info(info_flags::numeric)) { numeric_flag = false; - ++r; + break; + } } // The pure numeric case is traditionally rather common. Hence, it is @@ -896,10 +939,11 @@ ex matrix::charpoly(const ex & lambda) const /** Inverse of this matrix. * + * @param algo selects the algorithm (one of solve_algo) * @return the inverted matrix * @exception logic_error (matrix not square) * @exception runtime_error (singular matrix) */ -matrix matrix::inverse() const +matrix matrix::inverse(unsigned algo) const { if (row != col) throw (std::logic_error("matrix::inverse(): matrix not square")); @@ -922,7 +966,7 @@ matrix matrix::inverse() const matrix sol(row,col); try { - sol = this->solve(vars,identity); + sol = this->solve(vars, identity, algo); } catch (const std::runtime_error & e) { if (e.what()==std::string("matrix::solve(): inconsistent linear system")) throw (std::runtime_error("matrix::inverse(): singular matrix")); @@ -938,14 +982,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(); @@ -970,11 +1015,11 @@ matrix matrix::solve(const matrix & vars, // Gather some statistical information about the augmented matrix: bool numeric_flag = true; - exvector::const_iterator r = aug.m.begin(), rend = aug.m.end(); - while (r!=rend && numeric_flag==true) { - if (!r->info(info_flags::numeric)) + for (auto & r : aug.m) { + if (!r.info(info_flags::numeric)) { numeric_flag = false; - ++r; + break; + } } // Here is the heuristics in case this routine has to decide: @@ -1009,11 +1054,11 @@ matrix matrix::solve(const matrix & vars, unsigned last_assigned_sol = n+1; for (int r=m-1; r>=0; --r) { unsigned fnz = 1; // first non-zero in row - while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].is_zero())) + while ((fnz<=n) && (aug.m[r*(n+p)+(fnz-1)].normal().is_zero())) ++fnz; if (fnz>n) { // row consists only of zeros, corresponding rhs must be 0, too - if (!aug.m[r*(n+p)+n+co].is_zero()) { + if (!aug.m[r*(n+p)+n+co].normal().is_zero()) { throw (std::runtime_error("matrix::solve(): inconsistent linear system")); } } else { @@ -1038,6 +1083,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 @@ -1046,7 +1114,7 @@ matrix matrix::solve(const matrix & vars, * more than once. According to W.M.Gentleman and S.C.Johnson this algorithm * is better than elimination schemes for matrices of sparse multivariate * polynomials and also for matrices of dense univariate polynomials if the - * matrix' dimesion is larger than 7. + * matrix' dimension is larger than 7. * * @return the determinant as a new expression (in expanded form) * @see matrix::determinant() */ @@ -1152,9 +1220,8 @@ ex matrix::determinant_minor() const for (unsigned j=fc; j 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 +1278,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 +1305,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 +1330,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,39 +1397,46 @@ 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 - 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) { - ex nd = cit->normal().to_rational(srl).numer_denom(); - ++cit; + exmap srl; // symbol replacement list + auto tmp_n_it = tmp_n.m.begin(), tmp_d_it = tmp_d.m.begin(); + for (auto & it : this->m) { + ex nd = it.normal().to_rational(srl).numer_denom(); *tmp_n_it++ = nd.op(0); *tmp_d_it++ = nd.op(1); } unsigned r0 = 0; - for (unsigned r1=0; (r1=0) { - if (indx>0) { + } else { + if (indx>r0) { + // Matrix needs pivoting, swap rows r0 and indx of tmp_n and tmp_d. 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(); tmp_d_it = tmp_d.m.begin(); - while (it != itend) - *it++ = ((*tmp_n_it++)/(*tmp_d_it++)).subs(srl, subs_options::no_pattern); + for (auto & it : this->m) + it = ((*tmp_n_it++)/(*tmp_d_it++)).subs(srl, subs_options::no_pattern); return sign; } @@ -1399,7 +1488,7 @@ int matrix::fraction_free_elimination(const bool det) * @param co is the column to be inspected * @param symbolic signal if we want the first non-zero element to be pivoted * (true) or the one with the largest absolute value (false). - * @return 0 if no interchange occured, -1 if all are zero (usually signaling + * @return 0 if no interchange occurred, -1 if all are zero (usually signaling * a degeneracy) and positive integer k means that rows ro and k were swapped. */ int matrix::pivot(unsigned ro, unsigned co, bool symbolic) @@ -1440,28 +1529,38 @@ int matrix::pivot(unsigned ro, unsigned co, bool symbolic) return k; } -ex lst_to_matrix(const lst & l) +/** Function to check that all elements of the matrix are zero. + */ +bool matrix::is_zero_matrix() const { - lst::const_iterator itr, itc; + for (auto & i : m) + if (!i.is_zero()) + return false; + return true; +} +ex lst_to_matrix(const lst & l) +{ // Find number of rows and columns size_t rows = l.nops(), cols = 0; - for (itr = l.begin(); itr != l.end(); ++itr) { - if (!is_a(*itr)) + for (auto & itr : l) { + if (!is_a(itr)) throw (std::invalid_argument("lst_to_matrix: argument must be a list of lists")); - if (itr->nops() > cols) - cols = itr->nops(); + if (itr.nops() > cols) + cols = itr.nops(); } // Allocate and fill matrix - matrix &M = *new matrix(rows, cols); - M.setflag(status_flags::dynallocated); - - unsigned i; - for (itr = l.begin(), i = 0; itr != l.end(); ++itr, ++i) { - unsigned j; - for (itc = ex_to(*itr).begin(), j = 0; itc != ex_to(*itr).end(); ++itc, ++j) - M(i, j) = *itc; + matrix & M = dynallocate(rows, cols); + + unsigned i = 0; + for (auto & itr : l) { + unsigned j = 0; + for (auto & itc : ex_to(itr)) { + M(i, j) = itc; + ++j; + } + ++i; } return M; @@ -1469,24 +1568,40 @@ ex lst_to_matrix(const lst & l) ex diag_matrix(const lst & l) { - lst::const_iterator it; size_t dim = l.nops(); // Allocate and fill matrix - matrix &M = *new matrix(dim, dim); - M.setflag(status_flags::dynallocated); + matrix & M = dynallocate(dim, dim); - unsigned i; - for (it = l.begin(), i = 0; it != l.end(); ++it, ++i) - M(i, i) = *it; + unsigned i = 0; + for (auto & it : l) { + M(i, i) = it; + ++i; + } + + return M; +} + +ex diag_matrix(std::initializer_list l) +{ + size_t dim = l.size(); + + // Allocate and fill matrix + matrix & M = dynallocate(dim, dim); + + unsigned i = 0; + for (auto & it : l) { + M(i, i) = it; + ++i; + } return M; } ex unit_matrix(unsigned r, unsigned c) { - matrix &Id = *new matrix(r, c); - Id.setflag(status_flags::dynallocated); + matrix & Id = dynallocate(r, c); + Id.setflag(status_flags::evaluated); for (unsigned i=0; i(r, c); + M.setflag(status_flags::evaluated); bool long_format = (r > 10 || c > 10); bool single_row = (r == 1 || c == 1); @@ -1530,4 +1645,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 = dynallocate(rows, cols); + M.setflag(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 = dynallocate(nr, nc); + M.setflag(status_flags::evaluated); + + for (unsigned ro=0; ro