Homec4science
Diffusion lammps (lammm-master) Browse

lammps/lib/atc/Matrix.h
7855ef6ddac4lammm-master

Last Change
Raw File

Matrix.h
View Options

#ifndef MATRIX_H
#define MATRIX_H
#include "MatrixDef.h"
namespace ATC_matrix {
static const int myPrecision = 15;
/**
* @class Matrix
* @brief Base class for linear algebra subsystem
*/
template<typename T>
class Matrix
{
protected:
Matrix(const Matrix &c);
public:
Matrix() {}
virtual ~Matrix() {}
//* stream output functions
void print(ostream &o) const { o << to_string(); }
void print(ostream &o, const string &name) const;
friend ostream& operator<<(ostream &o, const Matrix<T> &m){m.print(o); return o;}
void print() const;
virtual void print(const string &name) const;
virtual string to_string() const;
// element by element operations
DenseMatrix<T> operator/(const Matrix<T>& B) const;
DenseMatrix<T> pow(int n) const;
DenseMatrix<T> pow(double n) const;
// functions that return a copy
DenseMatrix<T> transpose() const;
void row_partition(const set<int> & rowsIn, set<int> & rows, set<int> & colsC,
DenseMatrix<T> & A1, DenseMatrix<T> & A2, bool complement=true) const;
set<int> row_partition(const set<int> & rows,
DenseMatrix<T> & A1, DenseMatrix<T> & A2) const;
void map(const set<int>& rows, const set<int>& cols, DenseMatrix<T> & A) const;
void insert(const set<int>& rows, const set<int>& cols, const DenseMatrix<T> & A);
void assemble(const set<int>& rows, const set<int>& cols, const DenseMatrix<T> & A);
// matrix to scalar functions
T sum() const;
T stdev() const;
T max() const;
T min() const;
T maxabs() const;
T minabs() const;
T norm() const;
T norm_sq() const;
T mean() const;
T dot(const Matrix<T> &r) const;
T trace() const;
// row and column operations
T row_sum (INDEX i=0) const { return row(*this,i).sum(); }
T row_mean (INDEX i=0) const { return row(*this,i).mean(); }
T row_norm (INDEX i=0) const { return row(*this,i).norm(); }
T row_min (INDEX i=0) const { return row(*this,i).min(); }
T row_max (INDEX i=0) const { return row(*this,i).max(); }
T row_stdev(INDEX i=0) const { return row(*this,i).stdev(); }
T col_sum (INDEX i=0) const { return column(*this,i).sum(); }
T col_mean (INDEX i=0) const { return column(*this,i).mean(); }
T col_norm (INDEX i=0) const { return column(*this,i).norm(); }
T col_min (INDEX i=0) const { return column(*this,i).min(); }
T col_max (INDEX i=0) const { return column(*this,i).max(); }
T col_stdev(INDEX i=0) const { return column(*this,i).stdev(); }
// pure virtual functions (required to implement these) ---------------------
//* reference index operator
virtual T& operator()(INDEX i, INDEX j)=0;
//* value index operator
virtual T operator()(INDEX i, INDEX j)const=0;
//* value flat index operator
virtual T& operator [](INDEX i)=0;
//* reference flat index operator
virtual T operator [](INDEX i) const=0;
//* returns the # of rows
virtual INDEX nRows() const=0;
//* returns the # of columns
virtual INDEX nCols() const=0;
//* returns a pointer to the data (dangerous)
virtual T * ptr() const=0;
//* resizes the matrix, copy what fits default to OFF
virtual void resize(INDEX nRows, INDEX nCols=1, bool copy=false)=0;
//* resizes the matrix, zero it out default to ON
virtual void reset(INDEX nRows, INDEX nCols=1, bool zero=true)=0;
//* resizes and copies data
virtual void copy(const T * ptr, INDEX nRows, INDEX nCols=1)=0;
//* create restart file
virtual void write_restart(FILE *f) const=0;
//* writes a matlab command to recreate this in a variable named s
virtual void matlab(ostream &o, const string &s="M") const;
//* writes a mathematica command to recreate this in a variable named s
virtual void mathematica(ostream &o, const string &s="M") const;
// output to matlab, with variable name s
void matlab(const string &s="M") const;
// output to mathematica, with variable name s
void mathematica(const string &s="M") const;
Matrix<T>& operator+=(const Matrix &r);
Matrix<T>& operator-=(const Matrix &r);
Matrix<T>& operator*=(const Matrix<T>& R);
Matrix<T>& operator/=(const Matrix<T>& R);
Matrix<T>& operator+=(const T v);
Matrix<T>& operator-=(const T v);
Matrix<T>& operator*=(const T v);
Matrix<T>& operator/=(T v);
Matrix<T>& divide_zero_safe(const Matrix<T>& B);
Matrix<T>& operator=(const T &v);
Matrix<T>& operator=(const Matrix<T> &c);
virtual void set_all_elements_to(const T &v);
//* adds a matrix scaled by factor s to this one.
void add_scaled(const Matrix<T> &A, const T& s);
//* sets all elements to zero
Matrix& zero();
//* sets matrix to the identity
Matrix& identity(int nrows=0);
//* returns the total number of elements
virtual INDEX size() const;
//* returns true if (i,j) is within the range of the matrix
bool in_range(INDEX i, INDEX j) const;
//* returns true if the matrix size is rs x cs
bool is_size(INDEX rs, INDEX cs) const;
//* returns true if the matrix is square and not empty
bool is_square() const;
//* returns true if Matrix, m, is the same size as this
bool same_size(const Matrix &m) const;
//* returns true if Matrix a and Matrix b are the same size
static bool same_size(const Matrix<T> &a, const Matrix<T> &b);
//* returns true if Matrix a rows are equal to Matrix b cols
static bool cols_equals_rows(const Matrix<T> &a, const Matrix<T> &b);
//* checks if memory is contiguous, only can be false for clone vector
virtual bool memory_contiguous() const { return true; }
//* checks if all values are within the prescribed range
virtual bool check_range(T min, T max) const;
protected:
virtual void _set_equal(const Matrix<T> &r);
};
//* Matrix operations
//@{
//* Sets C as b*C + a*A[tranpose?]*B[transpose?]
template<typename T>
void MultAB(const Matrix<T> &A, const Matrix<T> &B, DenseMatrix<T> &C,
bool At=0, bool Bt=0, T a=1, T b=0);
//* performs a matrix-vector multiply
template<typename T>
void MultMv(const Matrix<T> &A, const Vector<T> &v, DenseVector<T> &c,
const bool At, T a, T b);
// returns the inverse of a double precision matrix
DenseMatrix<double> inv(const Matrix<double>& A);
// returns the eigensystem of a pair of double precision matrices
DenseMatrix<double> eigensystem(const Matrix<double>& A, const Matrix<double>& B, DenseMatrix<double> & eVals, bool normalize = true);
// returns the polar decomposition of a double precision matrix
DenseMatrix<double> polar_decomposition(const Matrix<double>& A, DenseMatrix<double> & rotation, DenseMatrix<double> & stretch, bool leftRotation = true);
//* returns the trace of a matrix
template<typename T>
T trace(const Matrix<T>& A) { return A.trace(); }
//* computes the determinant of a square matrix
double det(const Matrix<double>& A);
//* Returns the maximum eigenvalue of a matrix.
double max_eigenvalue(const Matrix<double>& A);
//@}
//-----------------------------------------------------------------------------
// computes the sum of the difference squared of each element.
//-----------------------------------------------------------------------------
template<typename T>
double sum_difference_squared(const Matrix<T>& A, const Matrix<T> &B)
{
SSCK(A, B, "sum_difference_squared");
double v=0.0;
for (INDEX i=0; i<A.size(); i++) {
double d = A[i]-B[i];
v += d*d;
}
return v;
}
//-----------------------------------------------------------------------------
//* Operator for Matrix-matrix product
//-----------------------------------------------------------------------------
template<typename T>
DenseMatrix<T> operator*(const Matrix<T> &A, const Matrix<T> &B)
{
DenseMatrix<T> C(0,0,false);
MultAB(A,B,C);
return C;
}
//-----------------------------------------------------------------------------
//* Multiply a Matrix by a scalar
//-----------------------------------------------------------------------------
template<typename T>
DenseMatrix<T> operator*(const Matrix<T> &M, const T s)
{
DenseMatrix<T> R(M);
return R*=s;
}
//-----------------------------------------------------------------------------
//* Multiply a Matrix by a scalar
template<typename T>
DenseMatrix<T> operator*(const T s, const Matrix<T> &M)
{
DenseMatrix<T> R(M);
return R*=s;
}
//-----------------------------------------------------------------------------
//* inverse scaling operator - must always create memory
template<typename T>
DenseMatrix<T> operator/(const Matrix<T> &M, const T s)
{
DenseMatrix<T> R(M);
return R*=(1.0/s); // for integer types this may be worthless
}
//-----------------------------------------------------------------------------
//* Operator for Matrix-matrix sum
template<typename T>
DenseMatrix<T> operator+(const Matrix<T> &A, const Matrix<T> &B)
{
DenseMatrix<T> C(A);
return C+=B;
}
//-----------------------------------------------------------------------------
//* Operator for Matrix-matrix subtraction
template<typename T>
DenseMatrix<T> operator-(const Matrix<T> &A, const Matrix<T> &B)
{
DenseMatrix<T> C(A);
return C-=B;
}
/******************************************************************************
* Template definitions for class Matrix
******************************************************************************/
//-----------------------------------------------------------------------------
//* performs a matrix-matrix multiply with general type implementation
template<typename T>
void MultAB(const Matrix<T> &A, const Matrix<T> &B, DenseMatrix<T> &C,
const bool At, const bool Bt, T a, T b)
{
const INDEX sA[2] = {A.nRows(), A.nCols()}; // m is sA[At] k is sA[!At]
const INDEX sB[2] = {B.nRows(), B.nCols()}; // k is sB[Bt] n is sB[!Bt]
const INDEX M=sA[At], K=sB[Bt], N=sB[!Bt]; // M is the number of rows in A or Atrans (sA[At]),
// K is the number of rows in B or Btrans (sB[Bt], sA[!At]),
// N is the number of columns in B or Btrans (sB[!Bt]).
GCK(A, B, sA[!At]!=K, "MultAB<T> shared index not equal size");
if (!C.is_size(M,N))
{
C.resize(M,N); // set size of C
C.zero();
}
else C *= b; // Zero C
for (INDEX p=0; p<M; p++) {
INDEX p_times_At = p*At;
INDEX p_times_notAt = p*!At;
for (INDEX q=0; q<N; q++) {
INDEX q_times_Bt = q*Bt;
INDEX q_times_notBt = q*!Bt;
for (INDEX r=0; r<K; r++) {
INDEX ai = p_times_notAt+r*At;
INDEX aj = p_times_At+r*!At;
INDEX bi = r*!Bt+q_times_Bt;
INDEX bj = r*Bt+q_times_notBt;
T a_entry = A(ai, aj);
T b_entry = B(bi, bj);
T mult = a_entry * b_entry;
C(p,q) += mult;
}
}
}
}
//-----------------------------------------------------------------------------
//* output operator
template<typename T>
string Matrix<T>::to_string() const
{
string s;
for (INDEX i=0; i<nRows(); i++)
{
if (i) s += '\n';
for (INDEX j=0; j<nCols(); j++)
{
if (j) s+= '\t';
s += ATC_Utility::to_string((*this)(i,j),myPrecision)+" ";
}
}
return s;
}
//-----------------------------------------------------------------------------
//* output operator that wraps the matrix in a nice labeled box
template<typename T>
void Matrix<T>::print(ostream &o, const string &name) const
{
o << "------- Begin "<<name<<" -----------------\n";
this->print(o);
o << "\n------- End "<<name<<" -------------------\n";
}
//-----------------------------------------------------------------------------
//* print operator, use cout by default
template<typename T>
void Matrix<T>::print() const
{
print(cout);
}
//-----------------------------------------------------------------------------
//* named print operator, use cout by default
template<typename T>
void Matrix<T>::print(const string &name) const
{
print(cout, name);
}
//-----------------------------------------------------------------------------
//* element by element division
template<typename T>
DenseMatrix<T> Matrix<T>::operator/ (const Matrix<T>& B) const
{
SSCK(*this, B, "Matrix<T>::Operator/");
DenseMatrix<T> R(*this);
R /= B;
return R;
}
//-----------------------------------------------------------------------------
//* element-wise raise to a power
template<typename T>
DenseMatrix<T> Matrix<T>::pow(int n) const
{
DenseMatrix<T> R(*this);
int sz=this->size(); for(INDEX i=0; i<sz; i++)
{
double val = R[i];
for (int k=1; k<n; k++) val *= R[i];
for (int k=n; k<1; k++) val /= R[i];
R[i] = val;
}
return R;
}
//-----------------------------------------------------------------------------
//* element-wise raise to a power
template<typename T>
DenseMatrix<T> Matrix<T>::pow(double n) const
{
DenseMatrix<T> R(*this);
int sz=this->size(); for(INDEX i=0; i<sz; i++)
{
double val = R[i];
R[i] = pow(val,n);
}
return R;
}
//-----------------------------------------------------------------------------
//* returns the transpose of this matrix (makes a copy)
template <typename T>
DenseMatrix<T> Matrix<T>::transpose() const
{
DenseMatrix<T> t(this->nCols(), this->nRows());
int szi = this->nRows();
int szj = this->nCols();
for (INDEX i = 0; i < szi; i++)
for (INDEX j = 0; j < szj; j++)
t(j,i) = (*this)(i,j);
return t;
}
//-----------------------------------------------------------------------------
//* returns the transpose of a matrix (makes a copy)
template <typename T>
DenseMatrix<T> transpose(const Matrix<T> &A)
{
return A.transpose();
}
//-----------------------------------------------------------------------------
//* Returns the sum of all matrix elements
template<typename T>
T Matrix<T>::sum() const
{
if (!size()) return T(0);
T v = (*this)[0];
for (INDEX i=1; i<this->size(); i++) v += (*this)[i];
return v;
}
//-----------------------------------------------------------------------------
//* Returns the standard deviation of the matrix
template<typename T>
T Matrix<T>::stdev() const
{
GCHK(this->size()<2, "Matrix::stdev() size must be > 1");
T mean = this->mean();
T diff = (*this)[0]-mean;
T stdev = diff*diff;
for (INDEX i=1; i<this->size(); i++)
{
diff = (*this)[i]-mean;
stdev += diff*diff;
}
return sqrt(stdev/T(this->size()-1));
}
//-----------------------------------------------------------------------------
//* Returns the maximum of the matrix
template<typename T>
T Matrix<T>::max() const
{
GCHK(!this->size(), "Matrix::max() size must be > 0");
T v = (*this)[0];
for (INDEX i=1; i<this->size(); i++) v = std::max(v, (*this)[i]);
return v;
}
//-----------------------------------------------------------------------------
//* Returns the minimum of the matrix
template<typename T>
T Matrix<T>::min() const
{
GCHK(!this->size(), "Matrix::min() size must be > 0");
T v = (*this)[0];
for (INDEX i=1; i<this->size(); i++) v = std::min(v, (*this)[i]);
return v;
}
//-----------------------------------------------------------------------------
//* Returns the maximum absolute value of the matrix
template<typename T>
T Matrix<T>::maxabs() const
{
GCHK(!this->size(), "Matrix::maxabs() size must be > 0");
T v = (*this)[0];
for (INDEX i=1; i<this->size(); i++) v = ATC_Utility::max_abs(v, (*this)[i]);
return v;
}
//-----------------------------------------------------------------------------
//* Returns the minimum absoute value of the matrix
template<typename T>
T Matrix<T>::minabs() const
{
GCHK(!this->size(), "Matrix::minabs() size must be > 0");
T v = (*this)[0];
for (INDEX i=1; i<this->size(); i++) v = ATC_Utility::min_abs(v, (*this)[i]);
return v;
}
//-----------------------------------------------------------------------------
//* returns the L2 norm of the matrix
template<typename T>
T Matrix<T>::norm() const
{
GCHK(!this->size(), "Matrix::norm() size must be > 0");
return sqrt(dot(*this));
}
//-----------------------------------------------------------------------------
//* returns the L2 norm of the matrix
template<typename T>
T Matrix<T>::norm_sq() const
{
GCHK(!this->size(), "Matrix::norm() size must be > 0");
return dot(*this);
}
//-----------------------------------------------------------------------------
//* returns the average of the matrix
template<typename T>
T Matrix<T>::mean() const
{
GCHK(!this->size(), "Matrix::mean() size must be > 0");
return sum()/T(this->size());
}
//-----------------------------------------------------------------------------
//* Returns the dot product of two vectors
template<typename T>
T Matrix<T>::dot(const Matrix<T>& r) const
{
SSCK(*this, r, "Matrix<T>::dot");
if (!this->size()) return T(0);
T v = r[0]*(*this)[0];
for (INDEX i=1; i<this->size(); i++) v += r[i]*(*this)[i];
return v;
}
//-----------------------------------------------------------------------------
// returns the sum of the matrix diagonal
//-----------------------------------------------------------------------------
template<typename T>
T Matrix<T>::trace() const
{
const INDEX N = std::min(nRows(),nCols());
if (!N) return T(0);
T r = (*this)(0,0);
for (INDEX i=0; i<N; i++)
r += (*this)(i,i);
return r;
}
//-----------------------------------------------------------------------------
//* Adds a matrix to this one
template<typename T>
Matrix<T>& Matrix<T>::operator+=(const Matrix &r)
{
SSCK(*this, r, "operator+= or operator +");
int sz=this->size(); for(INDEX i=0; i<sz; i++) (*this)[i]+=r[i];
return *this;
}
//-----------------------------------------------------------------------------
// subtracts a matrix from this one
//-----------------------------------------------------------------------------
template<typename T>
Matrix<T>& Matrix<T>::operator-=(const Matrix &r)
{
SSCK(*this, r, "operator-= or operator -");
int sz=this->size(); for(INDEX i=0; i<sz; i++) (*this)[i]-=r[i];
return *this;
}
//-----------------------------------------------------------------------------
// multiplies each element in this by the corresponding element in R
//-----------------------------------------------------------------------------
template<typename T>
Matrix<T>& Matrix<T>::operator*=(const Matrix<T>& R)
{
if ((R.nCols()==1) && (this->nCols()>1)) { // multiply every entry in a row by the same value
int szi = this->nRows();
int szj = this->nCols();
for (INDEX i = 0; i < szi; i++)
for (INDEX j = 0; j < szj; j++)
{
(*this)(i,j) *= R[i];
}
}
else if (((R.nCols()==R.size()) && (R.nRows()==R.size())) && !((this->nCols()==this->size()) && (this->nRows()==this->size()))){
int szi = this->nRows();
int szj = this->nCols();
for (INDEX i = 0; i < szi; i++)
for (INDEX j = 0; j < szj; j++)
{
(*this)(i,j) *= R[i];
}
}
else { // multiply each entry by a different value
int sz = this->size();
for (INDEX i = 0; i < sz; i++)
{
(*this)[i] *= R[i];
}
}
return *this;
}
//-----------------------------------------------------------------------------
// divides each element in this by the corresponding element in R
//-----------------------------------------------------------------------------
template<typename T>
Matrix<T>& Matrix<T>::operator/=(const Matrix<T>& R)
{
if ((R.nCols()==1) && (this->nCols()>1)) { // divide every entry in a row by the same value
int szi = this->nRows();
int szj = this->nCols();
for (INDEX i = 0; i < szi; i++)
for (INDEX j = 0; j < szj; j++)
{
(*this)(i,j) /= R[i];
}
}
else { // divide each entry by a different value
SSCK(*this, R, "operator/= or operator/");
int sz = this->size();
for(INDEX i = 0; i < sz; i++)
{
GCHK(fabs(R[i])==0,"Operator/: division by zero");
(*this)[i] /= R[i];
}
}
return *this;
}
//-----------------------------------------------------------------------------
// divides each element in this by the corresponding element in R unless zero
//-----------------------------------------------------------------------------
template<typename T>
Matrix<T>& Matrix<T>::divide_zero_safe(const Matrix<T>& R)
{
if ((R.nCols()==1) && (this->nCols()>1)) { // divide every entry in a row by the same value
int szi = this->nRows();
int szj = this->nCols();
for (INDEX i = 0; i < szi; i++)
for (INDEX j = 0; j < szj; j++)
{
if(fabs(R[i])!=0) {
(*this)(i,j) /= R[i];
}
}
}
else { // divide each entry by a different value
SSCK(*this, R, "operator/= or operator/");
int sz = this->size();
for(INDEX i = 0; i < sz; i++)
{
if(fabs(R[i])!=0) {
(*this)[i] /= R[i];
}
}
}
return *this;
}
//-----------------------------------------------------------------------------
// scales this matrix by a constant
//-----------------------------------------------------------------------------
template<typename T>
Matrix<T>& Matrix<T>::operator*=(const T v)
{
int sz=this->size(); for(INDEX i=0; i<sz; i++) (*this)[i]*=v;
return *this;
}
//-----------------------------------------------------------------------------
// adds a constant to this matrix
//-----------------------------------------------------------------------------
template<typename T>
Matrix<T>& Matrix<T>::operator+=(const T v)
{
int sz=this->size(); for(INDEX i=0; i<sz; i++) (*this)[i]+=v;
return *this;
}
//-----------------------------------------------------------------------------
// substracts a constant to this matrix
//-----------------------------------------------------------------------------
template<typename T>
Matrix<T>& Matrix<T>::operator-=(const T v)
{
int sz=this->size(); for(INDEX i=0; i<sz; i++) (*this)[i]-=v;
return *this;
}
//-----------------------------------------------------------------------------
//* scales this matrix by the inverse of a constant
template<typename T>
Matrix<T>& Matrix<T>::operator/=(T v)
{
return (*this)*=(1.0/v);
}
//----------------------------------------------------------------------------
// Assigns one matrix to another
//----------------------------------------------------------------------------
template<typename T>
Matrix<T>& Matrix<T>::operator=(const Matrix<T> &r)
{
this->_set_equal(r);
return *this;
}
//----------------------------------------------------------------------------
// general matrix assignment (for densely packed matrices)
//----------------------------------------------------------------------------
template<typename T>
void Matrix<T>::_set_equal(const Matrix<T> &r)
{
this->resize(r.nRows(), r.nCols());
const Matrix<T> *pr = &r;
if (const SparseMatrix<T> *ps = sparse_cast(pr))
copy_sparse_to_matrix(ps, *this);
else if (diag_cast(pr)) // r is Diagonal?
{
this->zero();
for (INDEX i=0; i<r.size(); i++) (*this)(i,i) = r[i];
}
else memcpy(this->ptr(), r.ptr(), r.size()*sizeof(T));
}
//-----------------------------------------------------------------------------
//* sets all elements to a constant
template<typename T>
inline Matrix<T>& Matrix<T>::operator=(const T &v)
{
set_all_elements_to(v);
return *this;
}
//-----------------------------------------------------------------------------
//* sets all elements to a constant
template<typename T>
void Matrix<T>::set_all_elements_to(const T &v)
{
int sz=this->size(); for(INDEX i=0; i<sz; i++) (*this)[i] = v;
}
//----------------------------------------------------------------------------
// adds a matrix scaled by factor s to this one.
//----------------------------------------------------------------------------
template <typename T>
void Matrix<T>::add_scaled(const Matrix<T> &A, const T& s)
{
SSCK(A, *this, "Matrix::add_scaled");
int sz=this->size(); for(INDEX i=0; i<sz; i++) (*this)[i] += A[i]*s;
}
//-----------------------------------------------------------------------------
//* writes a matlab command to the console
template<typename T>
void Matrix<T>::matlab(const string &s) const
{
this->matlab(cout, s);
}
//-----------------------------------------------------------------------------
//* Writes a matlab script defining the vector to the stream
template<typename T>
void Matrix<T>::matlab(ostream &o, const string &s) const
{
o << s <<"=zeros(" << nRows() << ","<<nCols()<<");\n";
int szi = this->nRows();
int szj = this->nCols();
for (INDEX i = 0; i < szi; i++)
for (INDEX j = 0; j < szj; j++)
o << s << "("<<i+1<<","<<j+1<<")=" << (*this)(i,j) << ";\n";
}
//-----------------------------------------------------------------------------
//* writes a mathematica command to the console
template<typename T>
void Matrix<T>::mathematica(const string &s) const
{
this->mathematica(cout, s);
}
//-----------------------------------------------------------------------------
//* Writes a mathematica script defining the vector to the stream
template<typename T>
void Matrix<T>::mathematica(ostream &o, const string &s) const
{
o << s <<" = { \n";
o.precision(15);
o << std::fixed;
for(INDEX i=0; i< nRows(); i++) {
o <<" { " << (*this)(i,0);
for(INDEX j=1; j< nCols(); j++) o << ", " << (*this)(i,j);
if (i+1 == nRows()) { o <<" } \n"; }
else { o <<" }, \n"; }
}
o << "};\n";
o << std::scientific;
}
//-----------------------------------------------------------------------------
//* sets all matrix elements to zero
template<typename T>
inline Matrix<T>& Matrix<T>::zero()
{
set_all_elements_to(T(0));
return *this;
}
//-----------------------------------------------------------------------------
//* sets to identity
template<typename T>
inline Matrix<T>& Matrix<T>::identity(int nrows)
{
if (nrows == 0) {
SQCK(*this, "DenseMatrix::inv(), matrix not square"); // check matrix is square
nrows = nRows();
}
reset(nrows,nrows);
for(INDEX i=0; i< nRows(); i++) (*this)(i,i) = 1;
return *this;
}
//-----------------------------------------------------------------------------
//* returns the total number of elements
template<typename T>
inline INDEX Matrix<T>::size() const
{
return nRows()*nCols();
}
//-----------------------------------------------------------------------------
//* returns true if (i,j) is within the range of the matrix
template<typename T>
inline bool Matrix<T>::in_range(INDEX i, INDEX j) const
{
return i<nRows() && j<nCols();
}
//-----------------------------------------------------------------------------
//* returns true if the matrix size is rs x cs
template<typename T>
inline bool Matrix<T>::is_size(INDEX rs, INDEX cs) const
{
return nRows()==rs && nCols()==cs;
}
//-----------------------------------------------------------------------------
//* returns true if the matrix is square and not empty
template<typename T>
inline bool Matrix<T>::is_square() const
{
return nRows()==nCols() && nRows();
}
//-----------------------------------------------------------------------------
//* returns true if Matrix, m, is the same size as this
template<typename T>
inline bool Matrix<T>::same_size(const Matrix<T> &m) const
{
return is_size(m.nRows(), m.nCols());
}
//-----------------------------------------------------------------------------
//* returns true if Matrix a and Matrix b are the same size
template<typename T>
inline bool Matrix<T>::same_size(const Matrix<T> &a, const Matrix<T> &b)
{
return a.same_size(b);
}
//-----------------------------------------------------------------------------
//* returns true if Matrix a rows = Matrix b cols
template<typename T>
inline bool Matrix<T>::cols_equals_rows(const Matrix<T> &a, const Matrix<T> &b)
{
return a.nCols() == b.nRows();
}
//-----------------------------------------------------------------------------
//* returns true if no value is outside of the range
template<typename T>
inline bool Matrix<T>::check_range(T min, T max) const
{
for (INDEX i = 0; i < this->nRows(); i++) {
for (INDEX j = 0; j < this->nCols(); j++) {
T val = (*this)(i,j);
if ( (val > max) || (val < min) ) return false;
}
}
return true;
}
//-----------------------------------------------------------------------------
//* Displays indexing error message and quits
template<typename T>
void ierror(const Matrix<T> &a, const char *FILE, int LINE, INDEX i, INDEX j)
{
cout << "Error: Matrix indexing failure ";
cout << "in file: " << FILE << ", line: "<< LINE <<"\n";
cout << "Tried accessing index (" << i << ", " << j <<")\n";
cout << "Matrix size was "<< a.nRows() << "x" << a.nCols() << "\n";
ERROR_FOR_BACKTRACE
exit(EXIT_FAILURE);
}
//-----------------------------------------------------------------------------
//* Displays custom message and indexing error and quits
template<typename T>
void ierror(const Matrix<T> &a, INDEX i, INDEX j, const string m)
{
cout << m << "\n";
cout << "Tried accessing index (" << i << ", " << j <<")\n";
cout << "Matrix size was "<< a.nRows() << "x" << a.nCols() << "\n";
ERROR_FOR_BACKTRACE
exit(EXIT_FAILURE);
}
//-----------------------------------------------------------------------------
//* Displays matrix compatibility error message
template<typename T>
void merror(const Matrix<T> &a, const Matrix<T> &b, const string m)
{
cout << "Error: " << m << "\n";
cout << "Matrix sizes were " << a.nRows() << "x" << a.nCols();
if (&a != &b) cout << ", and "<< b.nRows() << "x" << b.nCols();
cout << "\n";
if (a.size() < 100) a.print("Matrix");
ERROR_FOR_BACKTRACE
exit(EXIT_FAILURE);
}
//-----------------------------------------------------------------------------
//* returns upper or lower half of a partitioned matrix
//* A1 is the on-diagonal square matrix, A2 is the off-diagonal matrix
//* rowsIn is the rows to be placed in A1
//* rows is the map for A1, (rows,colsC) is the map for A2
template <typename T>
void Matrix<T>::row_partition(const set<int> & rowsIn,
set<int> & rows, set<int> & colsC,
DenseMatrix<T> & A1, DenseMatrix<T> & A2, bool complement) const
{
if (complement) {
for (INDEX i = 0; i < this->nRows(); i++) {
if (rowsIn.find(i) == rowsIn.end() ) rows.insert(i);
}
}
else rows = rowsIn;
// complement of set "rows" in set of this.cols is "cols"
for (INDEX i = 0; i < this->nCols(); i++) {
if (rows.find(i) == rows.end() ) colsC.insert(i);
}
// degenerate cases
if (int(rows.size()) == this->nCols()) {
A1 = (*this);
A2.reset(0,0);
return;
}
else if (rows.size() == 0) {
A1.reset(0,0);
A2 = (*this);
return;
}
// non-degenerate case
int nrows = rows.size();
int ncolsC = colsC.size();
A1.reset(nrows,nrows);
A2.reset(nrows,ncolsC);
set<int>::const_iterator itrI, itrJ;
INDEX i =0;
for (itrI = rows.begin(); itrI != rows.end(); itrI++) {
INDEX j = 0;
for (itrJ = rows.begin(); itrJ != rows.end(); itrJ++) {
A1(i,j) = (*this)(*itrI,*itrJ);
j++;
}
j = 0;
for (itrJ = colsC.begin(); itrJ != colsC.end(); itrJ++) {
A2(i,j) = (*this)(*itrI,*itrJ);
j++;
}
i++;
}
}
template <typename T>
set<int> Matrix<T>::row_partition(const set<int> & rows,
DenseMatrix<T> & A1, DenseMatrix<T> & A2) const
{
// complement of set "rows" in set of this.cols is "cols"
set<int> colsC;
for (INDEX i = 0; i < this->nCols(); i++) {
if (rows.find(i) == rows.end() ) colsC.insert(i);
}
// degenerate cases
if (int(rows.size()) == this->nCols()) {
A1 = (*this);
A2.reset(0,0);
return colsC;
}
else if (rows.size() == 0) {
A1.reset(0,0);
A2 = (*this);
return colsC;
}
// non-degenerate case
int nrows = rows.size();
int ncolsC = colsC.size();
A1.reset(nrows,nrows);
A2.reset(nrows,ncolsC);
set<int>::const_iterator itrI, itrJ;
INDEX i =0;
for (itrI = rows.begin(); itrI != rows.end(); itrI++) {
INDEX j = 0;
for (itrJ = rows.begin(); itrJ != rows.end(); itrJ++) {
A1(i,j) = (*this)(*itrI,*itrJ);
j++;
}
j = 0;
for (itrJ = colsC.begin(); itrJ != colsC.end(); itrJ++) {
A2(i,j) = (*this)(*itrI,*itrJ);
j++;
}
i++;
}
return colsC;
}
//-----------------------------------------------------------------------------
//* returns row & column mapped matrix
template <typename T>
void Matrix<T>::map(const set<int> & rows, const set<int> & cols,
DenseMatrix<T> & A ) const
{
if (rows.size() == 0 || cols.size() == 0 ) {
A.reset(0,0);
return;
}
int nrows = rows.size();
int ncols = cols.size();
A.reset(nrows,ncols);
set<int>::const_iterator itrI, itrJ;
INDEX i =0;
for (itrI = rows.begin(); itrI != rows.end(); itrI++) {
INDEX j = 0;
for (itrJ = cols.begin(); itrJ != cols.end(); itrJ++) {
A(i,j) = (*this)(*itrI,*itrJ);
j++;
}
i++;
}
}
//-----------------------------------------------------------------------------
//* inserts elements from a smaller matrix
template <typename T>
void Matrix<T>::insert(const set<int> & rows, const set<int> & cols,
const DenseMatrix<T> & A )
{
if (rows.size() == 0 || cols.size() == 0 ) return;
set<int>::const_iterator itrI, itrJ;
int i =0;
for (itrI = rows.begin(); itrI != rows.end(); itrI++) {
int j = 0;
for (itrJ = cols.begin(); itrJ != cols.end(); itrJ++) {
(*this)(*itrI,*itrJ) = A(i,j);
//cout << *itrI << " " << *itrJ << " : " << (*this)(*itrI,*itrJ) << "\n";
j++;
}
i++;
}
}
//-----------------------------------------------------------------------------
//* assemble elements from a smaller matrix
template <typename T>
void Matrix<T>::assemble(const set<int> & rows, const set<int> & cols,
const DenseMatrix<T> & A )
{
if (rows.size() == 0 || cols.size() == 0 ) return;
set<int>::const_iterator itrI, itrJ;
int i =0;
for (itrI = rows.begin(); itrI != rows.end(); itrI++) {
int j = 0;
for (itrJ = cols.begin(); itrJ != cols.end(); itrJ++) {
(*this)(*itrI,*itrJ) += A(i,j);
j++;
}
i++;
}
}
//-----------------------------------------------------------------------------
} // end namespace
#endif
c4science · Help