#include <stdexcept>
#include "matrix.h"
-#include "archive.h"
#include "numeric.h"
#include "lst.h"
-#include "utils.h"
-#include "debugmsg.h"
+#include "idx.h"
+#include "indexed.h"
#include "power.h"
#include "symbol.h"
#include "normal.h"
+#include "print.h"
+#include "archive.h"
+#include "utils.h"
+#include "debugmsg.h"
-#ifndef NO_NAMESPACE_GINAC
namespace GiNaC {
-#endif // ndef NO_NAMESPACE_GINAC
GINAC_IMPLEMENT_REGISTERED_CLASS(matrix, basic)
//////////
-// default constructor, destructor, copy constructor, assignment operator
-// and helpers:
+// default ctor, dtor, copy ctor, assignment operator and helpers:
//////////
-// public
-
/** Default ctor. Initializes to 1 x 1-dimensional zero-matrix. */
matrix::matrix() : inherited(TINFO_matrix), row(1), col(1)
{
- debugmsg("matrix default constructor",LOGLEVEL_CONSTRUCT);
+ debugmsg("matrix default ctor",LOGLEVEL_CONSTRUCT);
m.push_back(_ex0());
}
-// protected
-
void matrix::copy(const matrix & other)
{
inherited::copy(other);
m = other.m; // STL's vector copying invoked here
}
-void matrix::destroy(bool call_parent)
-{
- if (call_parent) inherited::destroy(call_parent);
-}
+DEFAULT_DESTROY(matrix)
//////////
-// other constructors
+// other ctors
//////////
// public
matrix::matrix(unsigned r, unsigned c)
: inherited(TINFO_matrix), row(r), col(c)
{
- debugmsg("matrix constructor from unsigned,unsigned",LOGLEVEL_CONSTRUCT);
+ debugmsg("matrix ctor from unsigned,unsigned",LOGLEVEL_CONSTRUCT);
m.resize(r*c, _ex0());
}
matrix::matrix(unsigned r, unsigned c, const exvector & m2)
: inherited(TINFO_matrix), row(r), col(c), m(m2)
{
- debugmsg("matrix constructor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT);
+ debugmsg("matrix ctor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT);
+}
+
+/** 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)
+{
+ debugmsg("matrix ctor from unsigned,unsigned,lst",LOGLEVEL_CONSTRUCT);
+ m.resize(r*c, _ex0());
+
+ for (unsigned i=0; i<l.nops(); i++) {
+ unsigned x = i % c;
+ unsigned y = i / c;
+ if (y >= r)
+ break; // matrix smaller than list: throw away excessive elements
+ m[y*c+x] = l.op(i);
+ }
}
//////////
// archiving
//////////
-/** Construct object from archive_node. */
matrix::matrix(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
{
- debugmsg("matrix constructor from archive_node", LOGLEVEL_CONSTRUCT);
+ debugmsg("matrix ctor from archive_node", LOGLEVEL_CONSTRUCT);
if (!(n.find_unsigned("row", row)) || !(n.find_unsigned("col", col)))
throw (std::runtime_error("unknown matrix dimensions in archive"));
m.reserve(row * col);
}
}
-/** Unarchive the object. */
-ex matrix::unarchive(const archive_node &n, const lst &sym_lst)
-{
- return (new matrix(n, sym_lst))->setflag(status_flags::dynallocated);
-}
-
-/** Archive the object. */
void matrix::archive(archive_node &n) const
{
inherited::archive(n);
}
}
+DEFAULT_UNARCHIVE(matrix)
+
//////////
// functions overriding virtual functions from bases classes
//////////
// public
-void matrix::print(std::ostream & os, unsigned upper_precedence) const
+void matrix::print(const print_context & c, unsigned level) const
{
- debugmsg("matrix print",LOGLEVEL_PRINT);
- os << "[[ ";
- for (unsigned r=0; r<row-1; ++r) {
- os << "[[";
- for (unsigned c=0; c<col-1; ++c)
- os << m[r*col+c] << ",";
- os << m[col*(r+1)-1] << "]], ";
- }
- os << "[[";
- for (unsigned c=0; c<col-1; ++c)
- os << m[(row-1)*col+c] << ",";
- os << m[row*col-1] << "]] ]]";
-}
+ debugmsg("matrix print", LOGLEVEL_PRINT);
+
+ if (is_of_type(c, print_tree)) {
+
+ inherited::print(c, level);
+
+ } else {
+
+ c.s << "[";
+ for (unsigned y=0; y<row-1; ++y) {
+ c.s << "[";
+ for (unsigned x=0; x<col-1; ++x) {
+ m[y*col+x].print(c);
+ c.s << ",";
+ }
+ m[col*(y+1)-1].print(c);
+ c.s << "],";
+ }
+ c.s << "[";
+ for (unsigned x=0; x<col-1; ++x) {
+ m[(row-1)*col+x].print(c);
+ c.s << ",";
+ }
+ m[row*col-1].print(c);
+ c.s << "]]";
-void matrix::printraw(std::ostream & os) const
-{
- debugmsg("matrix printraw",LOGLEVEL_PRINT);
- os << "matrix(" << row << "," << col <<",";
- for (unsigned r=0; r<row-1; ++r) {
- os << "(";
- for (unsigned c=0; c<col-1; ++c)
- os << m[r*col+c] << ",";
- os << m[col*(r-1)-1] << "),";
}
- os << "(";
- for (unsigned c=0; c<col-1; ++c)
- os << m[(row-1)*col+c] << ",";
- os << m[row*col-1] << "))";
}
/** nops is defined to be rows x columns. */
return matrix(row, col, tmp);
}
-/** Search ocurrences. A matrix 'has' an expression if it is the expression
- * itself or one of the elements 'has' it. */
-bool matrix::has(const ex & other) const
-{
- GINAC_ASSERT(other.bp!=0);
-
- // tautology: it is the expression itself
- if (is_equal(*other.bp)) return true;
-
- // search all the elements
- for (exvector::const_iterator r=m.begin(); r!=m.end(); ++r)
- if ((*r).has(other)) return true;
-
- return false;
-}
-
-/** evaluate matrix entry by entry. */
+/** Evaluate matrix entry by entry. */
ex matrix::eval(int level) const
{
debugmsg("matrix eval",LOGLEVEL_MEMBER_FUNCTION);
status_flags::evaluated );
}
-/** evaluate matrix numerically entry by entry. */
+/** Evaluate matrix numerically entry by entry. */
ex matrix::evalf(int level) const
{
debugmsg("matrix evalf",LOGLEVEL_MEMBER_FUNCTION);
return matrix(row, col, m2);
}
+ex matrix::subs(const lst & ls, const lst & lr, bool no_pattern) const
+{
+ exvector m2(row * col);
+ for (unsigned r=0; r<row; ++r)
+ for (unsigned c=0; c<col; ++c)
+ m2[r*col+c] = m[r*col+c].subs(ls, lr, no_pattern);
+
+ return ex(matrix(row, col, m2)).bp->basic::subs(ls, lr, no_pattern);
+}
+
// protected
int matrix::compare_same_type(const basic & other) const
return 0;
}
+/** Automatic symbolic evaluation of an indexed matrix. */
+ex matrix::eval_indexed(const basic & i) const
+{
+ GINAC_ASSERT(is_of_type(i, indexed));
+ GINAC_ASSERT(is_ex_of_type(i.op(0), matrix));
+
+ bool all_indices_unsigned = static_cast<const indexed &>(i).all_index_values_are(info_flags::nonnegint);
+
+ // Check indices
+ if (i.nops() == 2) {
+
+ // One index, must be one-dimensional vector
+ if (row != 1 && col != 1)
+ throw (std::runtime_error("matrix::eval_indexed(): vector must have exactly 1 index"));
+
+ const idx & i1 = ex_to_idx(i.op(1));
+
+ if (col == 1) {
+
+ // Column vector
+ if (!i1.get_dim().is_equal(row))
+ throw (std::runtime_error("matrix::eval_indexed(): dimension of index must match number of vector elements"));
+
+ // Index numeric -> return vector element
+ if (all_indices_unsigned) {
+ unsigned n1 = ex_to_numeric(i1.get_value()).to_int();
+ if (n1 >= row)
+ throw (std::runtime_error("matrix::eval_indexed(): value of index exceeds number of vector elements"));
+ return (*this)(n1, 0);
+ }
+
+ } else {
+
+ // Row vector
+ if (!i1.get_dim().is_equal(col))
+ throw (std::runtime_error("matrix::eval_indexed(): dimension of index must match number of vector elements"));
+
+ // Index numeric -> return vector element
+ if (all_indices_unsigned) {
+ unsigned n1 = ex_to_numeric(i1.get_value()).to_int();
+ if (n1 >= col)
+ throw (std::runtime_error("matrix::eval_indexed(): value of index exceeds number of vector elements"));
+ return (*this)(0, n1);
+ }
+ }
+
+ } else if (i.nops() == 3) {
+
+ // Two indices
+ const idx & i1 = ex_to_idx(i.op(1));
+ const idx & i2 = ex_to_idx(i.op(2));
+
+ if (!i1.get_dim().is_equal(row))
+ throw (std::runtime_error("matrix::eval_indexed(): dimension of first index must match number of rows"));
+ if (!i2.get_dim().is_equal(col))
+ throw (std::runtime_error("matrix::eval_indexed(): dimension of second index must match number of columns"));
+
+ // Pair of dummy indices -> compute trace
+ if (is_dummy_pair(i1, i2))
+ return trace();
+
+ // Both indices numeric -> return matrix element
+ if (all_indices_unsigned) {
+ unsigned n1 = ex_to_numeric(i1.get_value()).to_int(), n2 = ex_to_numeric(i2.get_value()).to_int();
+ if (n1 >= row)
+ throw (std::runtime_error("matrix::eval_indexed(): value of first index exceeds number of rows"));
+ if (n2 >= col)
+ throw (std::runtime_error("matrix::eval_indexed(): value of second index exceeds number of columns"));
+ return (*this)(n1, n2);
+ }
+
+ } else
+ throw (std::runtime_error("matrix::eval_indexed(): matrix must have exactly 2 indices"));
+
+ return i.hold();
+}
+
+/** Sum of two indexed matrices. */
+ex matrix::add_indexed(const ex & self, const ex & other) const
+{
+ GINAC_ASSERT(is_ex_of_type(self, indexed));
+ GINAC_ASSERT(is_ex_of_type(self.op(0), matrix));
+ GINAC_ASSERT(is_ex_of_type(other, indexed));
+ GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
+
+ // Only add two matrices
+ if (is_ex_of_type(other.op(0), matrix)) {
+ GINAC_ASSERT(other.nops() == 2 || other.nops() == 3);
+
+ const matrix &self_matrix = ex_to_matrix(self.op(0));
+ const matrix &other_matrix = ex_to_matrix(other.op(0));
+
+ if (self.nops() == 2 && other.nops() == 2) { // vector + vector
+
+ if (self_matrix.row == other_matrix.row)
+ return indexed(self_matrix.add(other_matrix), self.op(1));
+ else if (self_matrix.row == other_matrix.col)
+ return indexed(self_matrix.add(other_matrix.transpose()), self.op(1));
+
+ } else if (self.nops() == 3 && other.nops() == 3) { // matrix + matrix
+
+ if (self.op(1).is_equal(other.op(1)) && self.op(2).is_equal(other.op(2)))
+ return indexed(self_matrix.add(other_matrix), self.op(1), self.op(2));
+ else if (self.op(1).is_equal(other.op(2)) && self.op(2).is_equal(other.op(1)))
+ return indexed(self_matrix.add(other_matrix.transpose()), self.op(1), self.op(2));
+
+ }
+ }
+
+ // Don't know what to do, return unevaluated sum
+ return self + other;
+}
+
+/** Product of an indexed matrix with a number. */
+ex matrix::scalar_mul_indexed(const ex & self, const numeric & other) const
+{
+ GINAC_ASSERT(is_ex_of_type(self, indexed));
+ GINAC_ASSERT(is_ex_of_type(self.op(0), matrix));
+ GINAC_ASSERT(self.nops() == 2 || self.nops() == 3);
+
+ const matrix &self_matrix = ex_to_matrix(self.op(0));
+
+ if (self.nops() == 2)
+ return indexed(self_matrix.mul(other), self.op(1));
+ else // self.nops() == 3
+ return indexed(self_matrix.mul(other), self.op(1), self.op(2));
+}
+
+/** Contraction of an indexed matrix with something else. */
+bool matrix::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+ GINAC_ASSERT(is_ex_of_type(*self, indexed));
+ GINAC_ASSERT(is_ex_of_type(*other, indexed));
+ GINAC_ASSERT(self->nops() == 2 || self->nops() == 3);
+ GINAC_ASSERT(is_ex_of_type(self->op(0), matrix));
+
+ // Only contract with other matrices
+ if (!is_ex_of_type(other->op(0), matrix))
+ return false;
+
+ GINAC_ASSERT(other->nops() == 2 || other->nops() == 3);
+
+ const matrix &self_matrix = ex_to_matrix(self->op(0));
+ const matrix &other_matrix = ex_to_matrix(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) {
+ // Column vector * column vector, transpose first vector
+ *self = self_matrix.transpose().mul(other_matrix)(0, 0);
+ } else {
+ // Column vector * row vector, swap factors
+ *self = other_matrix.mul(self_matrix)(0, 0);
+ }
+ } else {
+ if (other_matrix.col == 1) {
+ // Row vector * column vector, perfect
+ *self = self_matrix.mul(other_matrix)(0, 0);
+ } else {
+ // Row vector * row vector, transpose second vector
+ *self = self_matrix.mul(other_matrix.transpose())(0, 0);
+ }
+ }
+ *other = _ex1();
+ return true;
+
+ } else { // vector * matrix
+
+ // B_i * A_ij = (B*A)_j (B is row vector)
+ if (is_dummy_pair(self->op(1), other->op(1))) {
+ if (self_matrix.row == 1)
+ *self = indexed(self_matrix.mul(other_matrix), other->op(2));
+ else
+ *self = indexed(self_matrix.transpose().mul(other_matrix), other->op(2));
+ *other = _ex1();
+ return true;
+ }
+
+ // B_j * A_ij = (A*B)_i (B is column vector)
+ if (is_dummy_pair(self->op(1), other->op(2))) {
+ if (self_matrix.col == 1)
+ *self = indexed(other_matrix.mul(self_matrix), other->op(1));
+ else
+ *self = indexed(other_matrix.mul(self_matrix.transpose()), other->op(1));
+ *other = _ex1();
+ return true;
+ }
+ }
+
+ } else if (other->nops() == 3) { // matrix * matrix
+
+ // A_ij * B_jk = (A*B)_ik
+ if (is_dummy_pair(self->op(2), other->op(1))) {
+ *self = indexed(self_matrix.mul(other_matrix), self->op(1), other->op(2));
+ *other = _ex1();
+ return true;
+ }
+
+ // A_ij * B_kj = (A*Btrans)_ik
+ if (is_dummy_pair(self->op(2), other->op(2))) {
+ *self = indexed(self_matrix.mul(other_matrix.transpose()), self->op(1), other->op(1));
+ *other = _ex1();
+ return true;
+ }
+
+ // A_ji * B_jk = (Atrans*B)_ik
+ if (is_dummy_pair(self->op(1), other->op(1))) {
+ *self = indexed(self_matrix.transpose().mul(other_matrix), self->op(2), other->op(2));
+ *other = _ex1();
+ return true;
+ }
+
+ // A_ji * B_kj = (B*A)_ki
+ if (is_dummy_pair(self->op(1), other->op(2))) {
+ *self = indexed(other_matrix.mul(self_matrix), other->op(1), self->op(2));
+ *other = _ex1();
+ return true;
+ }
+ }
+
+ return false;
+}
+
+
//////////
// non-virtual functions in this class
//////////
matrix matrix::add(const matrix & other) const
{
if (col != other.col || row != other.row)
- throw (std::logic_error("matrix::add(): incompatible matrices"));
+ throw std::logic_error("matrix::add(): incompatible matrices");
exvector sum(this->m);
exvector::iterator i;
matrix matrix::sub(const matrix & other) const
{
if (col != other.col || row != other.row)
- throw (std::logic_error("matrix::sub(): incompatible matrices"));
+ throw std::logic_error("matrix::sub(): incompatible matrices");
exvector dif(this->m);
exvector::iterator i;
matrix matrix::mul(const matrix & other) const
{
if (this->cols() != other.rows())
- throw (std::logic_error("matrix::mul(): incompatible matrices"));
+ throw std::logic_error("matrix::mul(): incompatible matrices");
exvector prod(this->rows()*other.cols());
}
-/** operator() to access elements.
+/** Product of matrix and scalar. */
+matrix matrix::mul(const numeric & other) const
+{
+ exvector prod(row * col);
+
+ for (unsigned r=0; r<row; ++r)
+ for (unsigned c=0; c<col; ++c)
+ prod[r*col+c] = m[r*col+c] * other;
+
+ return matrix(row, col, prod);
+}
+
+
+/** Product of matrix and scalar expression. */
+matrix matrix::mul_scalar(const ex & other) const
+{
+ if (other.return_type() != return_types::commutative)
+ throw std::runtime_error("matrix::mul_scalar(): non-commutative scalar");
+
+ exvector prod(row * col);
+
+ for (unsigned r=0; r<row; ++r)
+ for (unsigned c=0; c<col; ++c)
+ prod[r*col+c] = m[r*col+c] * other;
+
+ return matrix(row, col, prod);
+}
+
+
+/** Power of a matrix. Currently handles integer exponents only. */
+matrix matrix::pow(const ex & expn) const
+{
+ if (col!=row)
+ throw (std::logic_error("matrix::pow(): matrix not square"));
+
+ if (is_ex_exactly_of_type(expn, numeric)) {
+ // Integer cases are computed by successive multiplication, using the
+ // obvious shortcut of storing temporaries, like A^4 == (A*A)*(A*A).
+ if (expn.info(info_flags::integer)) {
+ numeric k;
+ matrix prod(row,col);
+ if (expn.info(info_flags::negative)) {
+ k = -ex_to_numeric(expn);
+ prod = this->inverse();
+ } else {
+ k = ex_to_numeric(expn);
+ prod = *this;
+ }
+ matrix result(row,col);
+ for (unsigned r=0; r<row; ++r)
+ result(r,r) = _ex1();
+ numeric b(1);
+ // this loop computes the representation of k in base 2 and
+ // multiplies the factors whenever needed:
+ while (b.compare(k)<=0) {
+ b *= numeric(2);
+ numeric r(mod(k,b));
+ if (!r.is_zero()) {
+ k -= r;
+ result = result.mul(prod);
+ }
+ if (b.compare(k)<=0)
+ prod = prod.mul(prod);
+ }
+ return result;
+ }
+ }
+ throw (std::runtime_error("matrix::pow(): don't know how to handle exponent"));
+}
+
+
+/** operator() to access elements for reading.
*
* @param ro row of element
* @param co column of element
}
-/** Set individual elements manually.
+/** operator() to access elements for writing.
*
+ * @param ro row of element
+ * @param co column of element
* @exception range_error (index out of range) */
-matrix & matrix::set(unsigned ro, unsigned co, ex value)
+ex & matrix::operator() (unsigned ro, unsigned co)
{
if (ro>=row || co>=col)
- throw (std::range_error("matrix::set(): index out of range"));
-
+ throw (std::range_error("matrix::operator(): index out of range"));
+
ensure_if_modifiable();
- m[ro*col+co] = value;
- return *this;
+ clearflag(status_flags::hash_calculated);
+ return m[ro*col+co];
}
return matrix(this->cols(),this->rows(),trans);
}
-
/** Determinant of square matrix. This routine doesn't actually calculate the
* determinant, it only implements some heuristics about which algorithm to
* run. If all the elements of the matrix are elements of an integral domain
std::vector<unsigned> pre_sort;
for (std::vector<uintpair>::iterator i=c_zeros.begin(); i!=c_zeros.end(); ++i)
pre_sort.push_back(i->second);
- int sign = permutation_sign(pre_sort);
+ std::vector<unsigned> 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<unsigned>::iterator i=pre_sort.begin();
if (row != col)
throw (std::logic_error("matrix::inverse(): matrix not square"));
- // NOTE: the Gauss-Jordan elimination used here can in principle be
- // replaced by two clever calls to gauss_elimination() and some to
- // transpose(). Wouldn't be more efficient (maybe less?), just more
- // orthogonal.
- matrix tmp(row,col);
- // set tmp to the unit matrix
- for (unsigned i=0; i<col; ++i)
- tmp.m[i*col+i] = _ex1();
+ // This routine actually doesn't do anything fancy at all. We compute the
+ // inverse of the matrix A by solving the system A * A^{-1} == Id.
- // create a copy of this matrix
- matrix cpy(*this);
- for (unsigned r1=0; r1<row; ++r1) {
- int indx = cpy.pivot(r1, r1);
- if (indx == -1) {
+ // First populate the identity matrix supposed to become the right hand side.
+ matrix identity(row,col);
+ for (unsigned i=0; i<row; ++i)
+ identity(i,i) = _ex1();
+
+ // Populate a dummy matrix of variables, just because of compatibility with
+ // matrix::solve() which wants this (for compatibility with under-determined
+ // systems of equations).
+ matrix vars(row,col);
+ for (unsigned r=0; r<row; ++r)
+ for (unsigned c=0; c<col; ++c)
+ vars(r,c) = symbol();
+
+ matrix sol(row,col);
+ try {
+ sol = this->solve(vars,identity);
+ } catch (const std::runtime_error & e) {
+ if (e.what()==std::string("matrix::solve(): inconsistent linear system"))
throw (std::runtime_error("matrix::inverse(): singular matrix"));
- }
- if (indx != 0) { // swap rows r and indx of matrix tmp
- for (unsigned i=0; i<col; ++i)
- tmp.m[r1*col+i].swap(tmp.m[indx*col+i]);
- }
- ex a1 = cpy.m[r1*col+r1];
- for (unsigned c=0; c<col; ++c) {
- cpy.m[r1*col+c] /= a1;
- tmp.m[r1*col+c] /= a1;
- }
- for (unsigned r2=0; r2<row; ++r2) {
- if (r2 != r1) {
- if (!cpy.m[r2*col+r1].is_zero()) {
- ex a2 = cpy.m[r2*col+r1];
- // yes, there is something to do in this column
- for (unsigned c=0; c<col; ++c) {
- cpy.m[r2*col+c] -= a2 * cpy.m[r1*col+c];
- if (!cpy.m[r2*col+c].info(info_flags::numeric))
- cpy.m[r2*col+c] = cpy.m[r2*col+c].normal();
- tmp.m[r2*col+c] -= a2 * tmp.m[r1*col+c];
- if (!tmp.m[r2*col+c].info(info_flags::numeric))
- tmp.m[r2*col+c] = tmp.m[r2*col+c].normal();
- }
- }
- }
- }
+ else
+ throw;
}
-
- return tmp;
+ return sol;
}
switch(algo) {
case solve_algo::gauss:
aug.gauss_elimination();
+ break;
case solve_algo::divfree:
aug.division_free_elimination();
+ break;
case solve_algo::bareiss:
default:
aug.fraction_free_elimination();
// assign solutions for vars between fnz+1 and
// last_assigned_sol-1: free parameters
for (unsigned c=fnz; c<last_assigned_sol-1; ++c)
- sol.set(c,co,vars.m[c*p+co]);
+ sol(c,co) = vars.m[c*p+co];
ex e = aug.m[r*(n+p)+n+co];
for (unsigned c=fnz; c<n; ++c)
e -= aug.m[r*(n+p)+c]*sol.m[c*p+co];
- sol.set(fnz-1,co,
- (e/(aug.m[r*(n+p)+(fnz-1)])).normal());
+ sol(fnz-1,co) = (e/(aug.m[r*(n+p)+(fnz-1)])).normal();
last_assigned_sol = fnz;
}
}
// assign solutions for vars between 1 and
// last_assigned_sol-1: free parameters
for (unsigned ro=0; ro<last_assigned_sol-1; ++ro)
- sol.set(ro,co,vars(ro,co));
+ sol(ro,co) = vars(ro,co);
}
return sol;
// for (unsigned r=0; r<minorM.rows(); ++r) {
// for (unsigned c=0; c<minorM.cols(); ++c) {
// if (r<r1)
- // minorM.set(r,c,m[r*col+c+1]);
+ // minorM(r,c) = m[r*col+c+1];
// else
- // minorM.set(r,c,m[(r+1)*col+c+1]);
+ // minorM(r,c) = m[(r+1)*col+c+1];
// }
// }
// // recurse down and care for sign:
return k;
}
-/** Convert list of lists to matrix. */
-ex lst_to_matrix(const ex &l)
+ex lst_to_matrix(const lst & l)
{
- if (!is_ex_of_type(l, lst))
- throw(std::invalid_argument("argument to lst_to_matrix() must be a lst"));
-
// Find number of rows and columns
unsigned rows = l.nops(), cols = 0, i, j;
for (i=0; i<rows; i++)
if (l.op(i).nops() > cols)
cols = l.op(i).nops();
-
+
// Allocate and fill matrix
matrix &m = *new matrix(rows, cols);
+ m.setflag(status_flags::dynallocated);
for (i=0; i<rows; i++)
for (j=0; j<cols; j++)
if (l.op(i).nops() > j)
- m.set(i, j, l.op(i).op(j));
+ m(i, j) = l.op(i).op(j);
else
- m.set(i, j, ex(0));
+ m(i, j) = _ex0();
+ return m;
+}
+
+ex diag_matrix(const lst & l)
+{
+ unsigned dim = l.nops();
+
+ matrix &m = *new matrix(dim, dim);
+ m.setflag(status_flags::dynallocated);
+ for (unsigned i=0; i<dim; i++)
+ m(i, i) = l.op(i);
+
return m;
}
-#ifndef NO_NAMESPACE_GINAC
} // namespace GiNaC
-#endif // ndef NO_NAMESPACE_GINAC
git://www.ginac.de / ginac.git / blobdiff