* Implementation of symbolic matrices */
/*
- * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2001 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
#include "symbol.h"
#include "normal.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)
+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());
}
-matrix::~matrix()
-{
- debugmsg("matrix destructor",LOGLEVEL_DESTRUCT);
-}
-
-matrix::matrix(const matrix & other)
-{
- debugmsg("matrix copy constructor",LOGLEVEL_CONSTRUCT);
- copy(other);
-}
-
-const matrix & matrix::operator=(const matrix & other)
-{
- debugmsg("matrix operator=",LOGLEVEL_ASSIGNMENT);
- if (this != &other) {
- destroy(1);
- copy(other);
- }
- return *this;
-}
-
// protected
+/** For use by copy ctor and assignment operator. */
void matrix::copy(const matrix & other)
{
inherited::copy(other);
}
//////////
-// other constructors
+// other ctors
//////////
// public
* @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)
{
- debugmsg("matrix constructor from unsigned,unsigned",LOGLEVEL_CONSTRUCT);
+ debugmsg("matrix ctor from unsigned,unsigned",LOGLEVEL_CONSTRUCT);
m.resize(r*c, _ex0());
}
/** 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)
{
- debugmsg("matrix constructor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT);
+ debugmsg("matrix ctor from unsigned,unsigned,exvector",LOGLEVEL_CONSTRUCT);
}
//////////
/** 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);
// public
-basic * matrix::duplicate() const
-{
- debugmsg("matrix duplicate",LOGLEVEL_DUPLICATE);
- return new matrix(*this);
-}
-
void matrix::print(std::ostream & os, unsigned upper_precedence) const
{
debugmsg("matrix print",LOGLEVEL_PRINT);
* @exception range_error (index out of range) */
const ex & matrix::operator() (unsigned ro, unsigned co) const
{
- if (ro<0 || ro>=row || co<0 || co>=col)
+ if (ro>=row || co>=col)
throw (std::range_error("matrix::operator(): index out of range"));
return m[ro*col+co];
* @exception range_error (index out of range) */
matrix & matrix::set(unsigned ro, unsigned co, ex value)
{
- if (ro<0 || ro>=row || co<0 || co>=col)
+ if (ro>=row || co>=col)
throw (std::range_error("matrix::set(): index out of range"));
-
+
ensure_if_modifiable();
m[ro*col+co] = value;
return *this;
throw (std::logic_error("matrix::inverse(): matrix not square"));
// NOTE: the Gauss-Jordan elimination used here can in principle be
- // replaced this by two clever calls to gauss_elimination() and some to
+ // 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);
}
for (unsigned r2=0; r2<row; ++r2) {
if (r2 != r1) {
- ex a2 = cpy.m[r2*col+r1];
- 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();
+ 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();
+ }
}
}
}
return (m[0]*m[3]-m[2]*m[1]).expand();
if (n==3)
return (m[0]*m[4]*m[8]-m[0]*m[5]*m[7]-
- m[1]*m[3]*m[8]+m[2]*m[3]*m[7]+
- m[1]*m[5]*m[6]-m[2]*m[4]*m[6]).expand();
+ m[1]*m[3]*m[8]+m[2]*m[3]*m[7]+
+ m[1]*m[5]*m[6]-m[2]*m[4]*m[6]).expand();
// This algorithm can best be understood by looking at a naive
// implementation of Laplace-expansion, like this one:
if (Pkey[fc-1]<fc+c)
break;
}
- if (fc<n-c)
+ if (fc<n-c && fc>0)
for (unsigned j=fc; j<n-c; ++j)
Pkey[j] = Pkey[j-1]+1;
} while(fc);
if (indx > 0)
sign = -sign;
for (unsigned r2=r0+1; r2<m; ++r2) {
- ex piv = this->m[r2*n+r1] / this->m[r0*n+r1];
- for (unsigned c=r1+1; c<n; ++c) {
- this->m[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();
+ if (!this->m[r2*n+r1].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; c<n; ++c) {
+ this->m[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 r2=r0+1; r2<m; ++r2) {
for (unsigned c=r1+1; c<n; ++c) {
dividend_n = (tmp_n.m[r0*n+r1]*tmp_n.m[r2*n+c]*
- tmp_d.m[r2*n+r1]*tmp_d.m[r0*n+c]
- -tmp_n.m[r2*n+r1]*tmp_n.m[r0*n+c]*
- tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
+ tmp_d.m[r2*n+r1]*tmp_d.m[r0*n+c]
+ -tmp_n.m[r2*n+r1]*tmp_n.m[r0*n+c]*
+ tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
dividend_d = (tmp_d.m[r2*n+r1]*tmp_d.m[r0*n+c]*
- tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
+ tmp_d.m[r0*n+r1]*tmp_d.m[r2*n+c]).expand();
bool check = divide(dividend_n, divisor_n,
- tmp_n.m[r2*n+c], true);
+ tmp_n.m[r2*n+c], true);
check &= divide(dividend_d, divisor_d,
- tmp_d.m[r2*n+c], true);
+ tmp_d.m[r2*n+c], true);
GINAC_ASSERT(check);
}
// fill up left hand side with zeros
// matrix needs pivoting, so swap rows k and ro
ensure_if_modifiable();
for (unsigned c=0; c<col; ++c)
- m[k*col+c].swap(m[ro*col+c]);
+ this->m[k*col+c].swap(this->m[ro*col+c]);
return k;
}
return m;
}
-//////////
-// global constants
-//////////
-
-const matrix some_matrix;
-const type_info & typeid_matrix=typeid(some_matrix);
-
-#ifndef NO_NAMESPACE_GINAC
} // namespace GiNaC
-#endif // ndef NO_NAMESPACE_GINAC
git://www.ginac.de / ginac.git / blobdiff