X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fmatrix.cpp;h=15ebd9aae23af79ceb10aec9f4b4bca73d1b1894;hp=9f9f67a82c83b62e76431de79ea4bf4054a05047;hb=aa6281216091efd92dc5fcc3f96c7189114e80f1;hpb=199b64938ab86af572d0816c15d7838730567b2d diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 9f9f67a8..15ebd9aa 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -141,7 +141,7 @@ void matrix::archive(archive_node &n) const DEFAULT_UNARCHIVE(matrix) ////////// -// functions overriding virtual functions from bases classes +// functions overriding virtual functions from base classes ////////// // public @@ -237,7 +237,7 @@ ex matrix::subs(const lst & ls, const lst & lr, bool no_pattern) const int matrix::compare_same_type(const basic & other) const { GINAC_ASSERT(is_exactly_of_type(other, matrix)); - const matrix & o = static_cast(const_cast(other)); + const matrix & o = static_cast(other); // compare number of rows if (row != o.rows()) @@ -259,6 +259,16 @@ int matrix::compare_same_type(const basic & other) const return 0; } +bool matrix::match_same_type(const basic & other) const +{ + GINAC_ASSERT(is_exactly_of_type(other, matrix)); + const matrix & o = static_cast(other); + + // The number of rows and columns must be the same. This is necessary to + // prevent a 2x3 matrix from matching a 3x2 one. + return row == o.rows() && col == o.cols(); +} + /** Automatic symbolic evaluation of an indexed matrix. */ ex matrix::eval_indexed(const basic & i) const { @@ -405,10 +415,8 @@ bool matrix::contract_with(exvector::iterator self, exvector::iterator other, ex const matrix &other_matrix = ex_to(other->op(0)); if (self->nops() == 2) { - unsigned self_dim = (self_matrix.col == 1) ? self_matrix.row : self_matrix.col; if (other->nops() == 2) { // vector * vector (scalar product) - unsigned other_dim = (other_matrix.col == 1) ? other_matrix.row : other_matrix.col; if (self_matrix.col == 1) { if (other_matrix.col == 1) { @@ -503,10 +511,10 @@ matrix matrix::add(const matrix & other) const throw std::logic_error("matrix::add(): incompatible matrices"); exvector sum(this->m); - exvector::iterator i; - exvector::const_iterator ci; - for (i=sum.begin(), ci=other.m.begin(); i!=sum.end(); ++i, ++ci) - (*i) += (*ci); + exvector::iterator i = sum.begin(), end = sum.end(); + exvector::const_iterator ci = other.m.begin(); + while (i != end) + *i++ += *ci++; return matrix(row,col,sum); } @@ -521,10 +529,10 @@ matrix matrix::sub(const matrix & other) const throw std::logic_error("matrix::sub(): incompatible matrices"); exvector dif(this->m); - exvector::iterator i; - exvector::const_iterator ci; - for (i=dif.begin(), ci=other.m.begin(); i!=dif.end(); ++i, ++ci) - (*i) -= (*ci); + exvector::iterator i = dif.begin(), end = dif.end(); + exvector::const_iterator ci = other.m.begin(); + while (i != end) + *i++ -= *ci++; return matrix(row,col,dif); } @@ -689,9 +697,10 @@ ex matrix::determinant(unsigned algo) const bool numeric_flag = true; bool normal_flag = false; unsigned sparse_count = 0; // counts non-zero elements - for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r) { + exvector::const_iterator r = m.begin(), rend = m.end(); + while (r != rend) { lst srl; // symbol replacement list - ex rtest = (*r).to_rational(srl); + ex rtest = r->to_rational(srl); if (!rtest.is_zero()) ++sparse_count; if (!rtest.info(info_flags::numeric)) @@ -699,6 +708,7 @@ ex matrix::determinant(unsigned algo) const if (!rtest.info(info_flags::crational_polynomial) && rtest.info(info_flags::rational_function)) normal_flag = true; + ++r; } // Here is the heuristics in case this routine has to decide: @@ -778,13 +788,13 @@ ex matrix::determinant(unsigned algo) const } sort(c_zeros.begin(),c_zeros.end()); std::vector pre_sort; - for (std::vector::iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i) + for (std::vector::const_iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i) 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::iterator i=pre_sort.begin(); + for (std::vector::const_iterator i=pre_sort.begin(); i!=pre_sort.end(); ++i,++c) { for (unsigned r=0; rinfo(info_flags::numeric)) numeric_flag = false; - } + ++r; } // The pure numeric case is traditionally rather common. Hence, it is @@ -950,9 +961,11 @@ matrix matrix::solve(const matrix & vars, // Gather some statistical information about the augmented matrix: bool numeric_flag = true; - for (exvector::const_iterator r=aug.m.begin(); r!=aug.m.end(); ++r) { - if (!(*r).info(info_flags::numeric)) + exvector::const_iterator r = aug.m.begin(), rend = aug.m.end(); + while (r != rend) { + if (!r->info(info_flags::numeric)) numeric_flag = false; + ++r; } // Here is the heuristics in case this routine has to decide: @@ -1299,13 +1312,13 @@ int matrix::fraction_free_elimination(const bool det) matrix tmp_n(*this); matrix tmp_d(m,n); // for denominators, if needed lst srl; // symbol replacement list - exvector::iterator it = this->m.begin(); - exvector::iterator tmp_n_it = tmp_n.m.begin(); - exvector::iterator tmp_d_it = tmp_d.m.begin(); - for (; it!= this->m.end(); ++it, ++tmp_n_it, ++tmp_d_it) { - (*tmp_n_it) = (*it).normal().to_rational(srl); - (*tmp_d_it) = (*tmp_n_it).denom(); - (*tmp_n_it) = (*tmp_n_it).numer(); + 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; + *tmp_n_it++ = nd.op(0); + *tmp_d_it++ = nd.op(1); } unsigned r0 = 0; @@ -1357,11 +1370,11 @@ int matrix::fraction_free_elimination(const bool det) } } // repopulate *this matrix: - it = this->m.begin(); + exvector::iterator it = this->m.begin(), itend = this->m.end(); tmp_n_it = tmp_n.m.begin(); tmp_d_it = tmp_d.m.begin(); - for (; it!= this->m.end(); ++it, ++tmp_n_it, ++tmp_d_it) - (*it) = ((*tmp_n_it)/(*tmp_d_it)).subs(srl); + while (it != itend) + *it++ = ((*tmp_n_it++)/(*tmp_d_it++)).subs(srl); return sign; }