Commit c84228af authored by Kirill Terekhov's avatar Kirill Terekhov

Some features for dense matrices

Made all algorithms for matrices abstract from matrix implementation.

Added few algorithms:
Solve - find A*X=B with A and B matrices
PseudoSolve - find X=A^+*B using singular value decomposition.
Kronecker product
operator / for B^{-1}*A
operator *= for in-place multiplication

Added ability to specify storage type for Matrix, allows for
Matrix<Storage::real,Storage::real_array> or
Matrix<Storage::var,Storage::var_array> to directly manipulate data
stored in the mesh.

Class SubMatrix with full matrix functionality to address and
manipulate some submatrix in original matrix.

Fully commented in doxygen inmost_dense.h.

Silence solver initializer when database file provided as “”.
parent d887be83
......@@ -6,24 +6,11 @@
#endif
#include <iomanip>
// Matrix with n rows and m columns
// __m__
// | |
// n| |
// |_____|
//
// todo:
// 1. expression templates for operations
// (???) how to for multiplication?
// 2. (ok) template matrix type for AD variables
// 3. template container type for data storage.
// 4. option for wrapper container around provided data storage. (to perform matrix operations with existing data)
// 5. class Subset for fortran-like access to matrix.
namespace INMOST
{
/// Structure that selects desired class, depending on the operation.
template<class A, class B> struct Promote;
template<> struct Promote<INMOST_DATA_REAL_TYPE, INMOST_DATA_REAL_TYPE> {typedef INMOST_DATA_REAL_TYPE type;};
#if defined(USE_AUTODIFF)
......@@ -39,18 +26,24 @@ namespace INMOST
__INLINE INMOST_DATA_REAL_TYPE get_value(INMOST_DATA_REAL_TYPE x) {return x;}
#endif
template<typename Var, typename Storage = array<Var> >
class SubMatrix;
template<typename Var, typename Storage = array<Var> >
class Matrix;
/// Abstract class for a matrix,
/// used to abstract away all the data storage and access
/// and provide common implementation of the algorithms.
template<typename Var>
class Matrix
class AbstractMatrix
{
public:
typedef unsigned enumerator;
protected:
array<Var> space;
enumerator n, m;
private:
/// Sign function for SVD algorithm.
static Var sign_func(const Var & a, const Var & b) {return (b >= 0.0 ? fabs(a) : -fabs(a));}
/// Max function for SVD algorithm.
static INMOST_DATA_REAL_TYPE max_func(INMOST_DATA_REAL_TYPE x, INMOST_DATA_REAL_TYPE y) { return x > y ? x : y; }
/// Function for QR rotation in SVD algorithm.
static Var pythag(const Var & a, const Var & b)
{
Var at = fabs(a), bt = fabs(b), ct, result;
......@@ -60,94 +53,98 @@ namespace INMOST
return result;
}
public:
bool CheckNans()
{
for(enumerator k = 0; k < n*m; ++k)
if( check_nans(space[k]) ) return true;
return false;
}
void RemoveRow(enumerator row)
{
for(enumerator k = row+1; k < n; ++k)
{
for(enumerator l = 0; l < m; ++l)
(*this)(k-1,l) = (*this)(k,l);
}
space.resize((n-1)*m);
--n;
}
void RemoveRows(enumerator first, enumerator last)
typedef unsigned enumerator;
/// Construct empty matrix.
AbstractMatrix() {}
/// Construct from matrix of the same type.
/// @param other Another matrix of the same type.
AbstractMatrix(const AbstractMatrix & b)
{
enumerator shift = last-first;
for(enumerator k = last+1; k < n; ++k)
{
for(enumerator l = 0; l < m; ++l)
(*this)(k-shift-1,l) = (*this)(k,l);
}
space.resize((n-shift)*m);
n-=shift;
Resize(b.Rows(),b.Cols());
for(enumerator i = 0; i < b.Rows(); ++i)
for(enumerator j = 0; j < b.Cols(); ++j)
(*this)(i,j) = b(i,j);
return *this;
}
void RemoveColumn(enumerator col)
/// Construct from matrix of another type.
/// @param other Another matrix of different type.
template<typename typeB>
AbstractMatrix(const AbstractMatrix<typeB> & b)
{
Matrix<Var> tmp(n,m-1);
for(enumerator k = 0; k < n; ++k)
{
for(enumerator l = 0; l < col; ++l)
tmp(k,l) = (*this)(k,l);
for(enumerator l = col+1; l < m; ++l)
tmp(k,l-1) = (*this)(k,l);
}
this->Swap(tmp);
Resize(b.Rows(),b.Cols());
for(enumerator i = 0; i < b.Rows(); ++i)
for(enumerator j = 0; j < b.Cols(); ++j)
assign((*this)(i,j),b(i,j));
return *this;
}
void RemoveColumns(enumerator first, enumerator last)
/// Assign matrix of the same type.
/// @param other Another matrix of the same type.
/// @return Reference to matrix.
AbstractMatrix & operator =(AbstractMatrix const & other)
{
enumerator shift = last-first;
Matrix<Var> tmp(n,m-shift);
for(enumerator k = 0; k < n; ++k)
{
for(enumerator l = 0; l < first; ++l)
tmp(k,l) = (*this)(k,l);
for(enumerator l = last+1; l < m; ++l)
tmp(k,l-shift-1) = (*this)(k,l);
}
this->Swap(tmp);
Resize(other.Rows(),other.Cols());
for(enumerator i = 0; i < other.Rows(); ++i)
for(enumerator j = 0; j < other.Cols(); ++j)
assign((*this)(i,j),other(i,j));
return *this;
}
void RemoveSubset(enumerator firstrow, enumerator lastrow, enumerator firstcol, enumerator lastcol)
/// Assign matrix of another type.
/// @param other Another matrix of different type.
/// @return Reference to matrix.
template<typename typeB>
AbstractMatrix & operator =(AbstractMatrix<typeB> const & other)
{
enumerator shiftrow = lastrow-firstrow;
enumerator shiftcol = lastcol-firstcol;
Matrix<Var> tmp(n-shiftrow, m-shiftcol);
for(enumerator k = 0; k < firstrow; ++k)
{
for(enumerator l = 0; l < firstcol; ++l)
tmp(k,l) = (*this)(k,l);
for(enumerator l = lastcol+1; l < m; ++l)
tmp(k,l-shiftcol-1) = (*this)(k,l);
}
for(enumerator k = lastrow+1; k < n; ++k)
{
for(enumerator l = 0; l < firstcol; ++l)
tmp(k-shiftrow-1,l) = (*this)(k,l);
for(enumerator l = lastcol+1; l < m; ++l)
tmp(k-shiftrow-1,l-shiftcol-1) = (*this)(k,l);
}
this->Swap(tmp);
Resize(other.Rows(),other.Cols());
for(enumerator i = 0; i < other.Rows(); ++i)
for(enumerator j = 0; j < other.Cols(); ++j)
assign((*this)(i,j),other(i,j));
return *this;
}
void Swap(Matrix & b)
/// Obtain number of rows.
/// @return Number of rows.
virtual enumerator Rows() const = 0;
/// Obtain number of columns.
/// @return Number of columns.
virtual enumerator Cols() const = 0;
/// Access element of the matrix by row and column indices.
/// @param i Column index.
/// @param j Row index.
/// @return Reference to element.
virtual Var & operator () (enumerator i, enumerator j) = 0;
/// Access element of the matrix by row and column indices
/// without right to change the element.
/// @param i Column index.
/// @param j Row index.
/// @return Reference to constant element.
virtual const Var & operator () (enumerator i, enumerator j) const = 0;
/// Resize the matrix into different size.
/// @param nrows New number of rows.
/// @param ncols New number of columns.
virtual void Resize(enumerator rows, enumerator cols) = 0;
/// Check all matrix entries for nans.
/// Also checks derivatives for matrices of variables.
bool CheckNans()
{
space.swap(b.space);
std::swap(n,b.n);
std::swap(m,b.m);
for(enumerator i = 0; i < Rows(); ++i)
for(enumerator j = 0; j < Cols(); ++j)
if( check_nans((*this)(i,j)) ) return true;
return false;
}
/// Singular value decomposition.
/// Reconstruct matrix: A = U*Sigma*V.Transpose().
/// Source http://www.public.iastate.edu/~dicook/JSS/paper/code/svd.c
/// With few corrections for robustness.
/// @param U Left unitary matrix, U^T U = I.
/// @param Sigma Diagonal matrix with singular values.
/// @param V Right unitary matrix, not transposed.
/// @param order_singular_values
bool SVD(Matrix & U, Matrix & Sigma, Matrix & V, bool order_singular_values = true)
/// \warning Somehow result differ in auto-differential mode.
/// \todo Test implementation for auto-differentiation.
bool SVD(AbstractMatrix<Var> & U, AbstractMatrix<Var> & Sigma, AbstractMatrix<Var> & V, bool order_singular_values = true) const
{
int flag, i, its, j, jj, k, l, nm;
int n = Rows();
int m = Cols();
Var c, f, h, s, x, y, z;
Var g = 0.0, scale = 0.0;
INMOST_DATA_REAL_TYPE anorm = 0.0;
......@@ -169,7 +166,6 @@ namespace INMOST
Sigma.Resize(m,m);
Sigma.Zero();
V.Resize(m,m);
std::swap(n,m); //this how original algorithm takes it
// Householder reduction to bidiagonal form
for (i = 0; i < (int)n; i++)
......@@ -434,22 +430,411 @@ namespace INMOST
std::swap(n,m);
return true;
}
Matrix() : space(0),n(0),m(0) {}
Matrix(Var * pspace, enumerator pn, enumerator pm) : space(pspace,pspace+pn*pm), n(pn), m(pm) {}
/// Set all the elements of the matrix to zero.
void Zero()
{
for(enumerator i = 0; i < Rows(); ++i)
for(enumerator j = 0; j < Cols(); ++j)
(*this)(i,j) = 0.0;
}
/// Transpose the current matrix.
/// @return Transposed matrix.
Matrix<Var> Transpose() const;
///Exchange contents of two matrices.
virtual void Swap(AbstractMatrix<Var> & b);
/// Unary minus. Change sign of each element of the matrix.
Matrix<Var> operator-() const;
/// Subtract a matrix.
/// @param other The matrix to be subtracted.
/// @return Difference of current matrix with another matrix.
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type> operator-(const AbstractMatrix<typeB> & other) const;
/// Subtract a matrix and store result in the current.
/// @param other The matrix to be subtracted.
/// @return Reference to the current matrix.
template<typename typeB>
AbstractMatrix & operator-=(const AbstractMatrix<typeB> & other);
/// Add a matrix.
/// @param other The matrix to be added.
/// @return Sum of current matrix with another matrix.
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type> operator+(const AbstractMatrix<typeB> & other) const;
/// Add a matrix and store result in the current.
/// @param other The matrix to be added.
/// @return Reference to the current matrix.
template<typename typeB>
AbstractMatrix & operator+=(const AbstractMatrix<typeB> & other);
/// Multiply the matrix by another matrix.
/// @param other Matrix to be multiplied from right.
/// @return Matrix multipled by another matrix.
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type> operator*(const AbstractMatrix<typeB> & other) const;
/// Multiply matrix with another matrix in-place.
/// @param B Another matrix to the right in multiplication.
/// @return Reference to current matrix.
template<typename typeB>
AbstractMatrix & operator*=(const AbstractMatrix<typeB> & B);
/// Multiply the matrix by a coefficient.
/// @param coef Coefficient.
/// @return Matrix multiplied by the coefficient.
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type> operator*(typeB coef) const;
/// Multiply the matrix by the coefficient of the same type and store the result.
/// @param coef Coefficient.
/// @return Reference to the current matrix.
template<typename typeB>
AbstractMatrix & operator*=(typeB coef);
/// Divide the matrix by a coefficient of a different type.
/// @param coef Coefficient.
/// @return Matrix divided by the coefficient.
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type>
operator/(typeB coef) const;
/// Divide the matrix by the coefficient of the same type and store the result.
/// @param coef Coefficient.
/// @return Reference to the current matrix.
template<typename typeB>
AbstractMatrix &
operator/=(typeB coef);
/// Performs B^{-1}*A, multiplication by inverse matrix from left.
/// @param other Matrix to be inverted and multiplied from left.
/// @return Multiplication of current matrix by inverse of another
/// matrix from left and boolean.
/// If boolean is true, then the matrix was inverted successfully.
/// \todo (test) Use Solve here.
template<typename typeB>
std::pair<Matrix<typename Promote<Var,typeB>::type>,bool>
operator/(const AbstractMatrix<typeB> & other) const
{
return other.Solve(*this);
}
/// Kronecker product, latex symbol \otimes.
/// @param other Matrix on the right of the Kronecker product.
/// @return Result of the Kronecker product of current and another matrix.
template<typename typeB>
Matrix<typename Promote<Var,typeB>::type>
Kronecker(const AbstractMatrix<typeB> & other) const;
/// Inverts matrix using Crout-LU decomposition with full pivoting for
/// maximum element. Works well for non-singular matrices, for singular
/// matrices look see Matrix::PseudoInvert.
/// @param print_fail Verbose output for singular matrices.
/// @return A pair of inverse matrix and boolean. If boolean is true,
/// then the matrix was inverted successfully.
/// @see Matrix::PseudoInvert.
/// \todo (test) Activate and test implementation with Solve.
std::pair<Matrix<Var>,bool> Invert(bool print_fail = false) const;
/// Finds X in A*X=B, where A and B are matrices.
/// Converts system into A^T*A*X=A^T*B.
/// Inverts matrix A^T*A using Crout-LU decomposition with full pivoting for
/// maximum element. Works well for non-singular matrices, for singular
/// matrices see Matrix::PseudoInvert.
/// @param print_fail Verbose output for singular matrices.
/// @return A pair of inverse matrix and boolean. If boolean is true,
/// then the matrix was inverted successfully.
/// @see Matrix::PseudoInvert.
/// \warning Number of rows in matrices A and B should match.
/// \todo
/// 1. Test implementation.
/// 2. Maximum product transversal + block pivoting instead of pivoting
/// by maximum element.
template<typename typeB>
std::pair<Matrix<typename Promote<Var,typeB>::type>,bool>
Solve(const AbstractMatrix<typeB> & B, bool print_fail = false) const;
/// Calculate sum of the diagonal elements of the matrix.
/// @return Trace of the matrix.
Var Trace() const
{
assert(Cols() == Rows());
Var ret = 0.0;
for(enumerator i = 0; i < Rows(); ++i) ret += (*this)(i,i);
return ret;
}
/// Output matrix to screen.
/// Does not output derivatices.
/// @param threshold Elements smaller then the number are considered to be zero.
void Print(INMOST_DATA_REAL_TYPE threshold = 1.0e-10) const
{
for(enumerator k = 0; k < Rows(); ++k)
{
for(enumerator l = 0; l < Cols(); ++l)
{
if( fabs(get_value((*this)(k,l))) > threshold )
std::cout << std::setw(10) << get_value((*this)(k,l));
else
std::cout << std::setw(10) << 0;
std::cout << " ";
}
std::cout << std::endl;
}
}
/// Check if the matrix is symmetric.
/// @return Returns true if the matrix is symmetric, otherwise false.
bool isSymmetric() const
{
if( Rows() != Cols() ) return false;
for(enumerator k = 0; k < Rows(); ++k)
{
for(enumerator l = k+1; l < Rows(); ++l)
if( fabs((*this)(k,l)-(*this)(l,k)) > 1.0e-7 )
return false;
}
return true;
}
/// Computes dot product by summing up multiplication of entries with the
/// same indices in the current and the provided matrix.
/// @param other Provided matrix.
/// @return Dot product of two matrices.
template<typename typeB>
typename Promote<Var,typeB>::type
DotProduct(const AbstractMatrix<typeB> & other) const
{
assert(Cols() == other.Cols());
assert(Rows() == other.Rows());
typename Promote<Var,typeB>::type ret = 0.0;
for(enumerator i = 0; i < Cols(); ++i)
for(enumerator j = 0; j < Rows(); ++j)
ret += ((*this)(i,j))*other(i,j);
return ret;
}
/// Computes dot product by summing up multiplication of entries with the
/// same indices in the current and the provided matrix.
/// @param other Provided matrix.
/// @return Dot product of two matrices.
template<typename typeB>
typename Promote<Var,typeB>::type operator ^(const AbstractMatrix<typeB> & other) const
{
return DotProduct(other);
}
/// Computes frobenious norm of the matrix.
/// @return Frobenius norm of the matrix.
Var FrobeniusNorm() const
{
Var ret = 0;
for(enumerator i = 0; i < Rows(); ++i)
for(enumerator j = 0; j < Cols(); ++j)
ret += (*this)(i,j)*(*this)(i,j);
return sqrt(ret);
}
/// Calculates Moore-Penrose pseudo-inverse of the matrix.
/// @param tol Thershold for singular values. Singular values smaller
/// then threshold are considered to be zero.
/// @param print_fail Verbose output if singular value decomposition
/// algorithm has failed.
/// @return A pair of pseudo-inverse matrix and boolean. If boolean is true,
/// then the matrix was inverted successfully.
std::pair<Matrix<Var>,bool> PseudoInvert(INMOST_DATA_REAL_TYPE tol = 0, bool print_fail = false) const;
/// Solves the system of equations of the form A*X=B, with A and B matrices.
/// Uses Moore-Penrose pseudo-inverse of the matrix A and calculates X = A^+*B.
/// @param B Matrix at the right hand side.
/// @param tol Thershold for singular values. Singular values smaller
/// then threshold are considered to be zero.
/// @param print_fail Verbose output if singular value decomposition
/// algorithm has failed.
/// @return A pair of the solution matrix X and boolean. If boolean is true,
/// then the matrix was inverted successfully.
template<typename typeB>
std::pair<Matrix<typename Promote<Var,typeB>::type>,bool>
PseudoSolve(const AbstractMatrix<typeB> & B, INMOST_DATA_REAL_TYPE tol = 0, bool print_fail = false) const;
/// Extract submatrix of a matrix.
/// Let A = {a_ij}, i in [0,n), j in [0,m) be original matrix.
/// Then the method returns B = {a_ij}, i in [ibeg,iend),
/// and j in [jbeg,jend).
/// @param ibeg Starting row in the original matrix.
/// @param iend Last row (excluded) in the original matrix.
/// @param jbeg Starting column in the original matrix.
/// @param jend Last column (excluded) in the original matrix.
/// @return Submatrix of the original matrix.
Matrix<Var> ExtractSubMatrix(enumerator ibeg, enumerator iend, enumerator jbeg, enumerator jend) const;
/// Change representation of the matrix into matrix of another size.
/// Useful to change representation from matrix into vector and back.
/// Replaces original number of columns and rows with a new one.
/// @return Matrix with same entries and provided number of rows and columns.
Matrix<Var> Repack(enumerator rows, enumerator cols) const;
};
/// Class for linear algebra operations on dense matrices.
/// Matrix with n rows and m columns.
///
/// __m__
/// | |
/// n| |
/// |_____|
///
/// \todo:
/// 1. expression templates for operations
/// (???) how to for multiplication?
/// efficient multiplication would require all the
/// matrix elements to be precomputed.
/// consider number 5 instead.
/// 2. (ok) template matrix type for AD variables
/// 3. (ok,test) template container type for data storage.
/// 4. (ok,test) option for wrapper container around provided data storage.
/// (to perform matrix operations with existing data)
/// 5. consider multi-threaded stack to get space for
/// matrices for local operations and returns.
/// 6. class SubMatrix for fortran-like access to matrix.
/// 7. Uniform implementation of algorithms for Matrix and Submatrix.
/// to achieve: make abdstract class with abstract element access
/// operator, make matrix and submatrix ancestors of that class
template<typename Var, typename storage_type>
class Matrix : public AbstractMatrix<Var>
{
public:
typedef unsigned enumerator; //< Integer type for indexes.
protected:
//array<Var> space; //< Array of row-wise stored elements.
storage_type space; //< Array of row-wise stored elements.
enumerator n; //< Number of rows.
enumerator m; //< Number of columns.
public:
/// Erase single row.
/// @param row Position of row, should be from 0 to Matrix::Rows()-1.
void RemoveRow(enumerator row)
{
assert(row < n);
for(enumerator k = row+1; k < n; ++k)
{
for(enumerator l = 0; l < m; ++l)
(*this)(k-1,l) = (*this)(k,l);
}
space.resize((n-1)*m);
--n;
}
/// Erase multiple rows. Assumes first <= last.
/// If first == last, then no row is deleted.
/// @param first Position of the first row, should be from 0 to Matrix::Rows()-1.
/// @param last Position behind the last row, should be from 0 to Matrix::Rows().
/// \todo check
void RemoveRows(enumerator first, enumerator last)
{
assert(first < n);
assert(last <= n);
assert(first <= last);
enumerator shift = last-first;
for(enumerator k = last; k < n; ++k)
{
for(enumerator l = 0; l < m; ++l)
(*this)(k-shift,l) = (*this)(k,l);
}
space.resize((n-shift)*m);
n-=shift;
}
/// Erase single column.
/// @param row Position of column, should be from 0 to Matrix::Cols()-1.
void RemoveColumn(enumerator col)
{
assert(col < m);
Matrix<Var> tmp(n,m-1);
for(enumerator k = 0; k < n; ++k)
{
for(enumerator l = 0; l < col; ++l)
tmp(k,l) = (*this)(k,l);
for(enumerator l = col+1; l < m; ++l)
tmp(k,l-1) = (*this)(k,l);
}
this->Swap(tmp);
}
/// Erase multiple columns. Assumes first <= last.
/// If first == last, then no column is deleted.
/// @param first Position of the first column, should be from 0 to Matrix::Cols()-1.
/// @param last Position behind the last column, should be from 0 to Matrix::Cols().
/// \todo check
void RemoveColumns(enumerator first, enumerator last)
{
assert(first < m);
assert(last <= m);
assert(first <= last);
enumerator shift = last-first;
Matrix<Var> tmp(n,m-shift);
for(enumerator k = 0; k < n; ++k)
{
for(enumerator l = 0; l < first; ++l)
tmp(k,l) = (*this)(k,l);
for(enumerator l = last; l < m; ++l)
tmp(k,l-shift) = (*this)(k,l);
}
this->Swap(tmp);
}
/// Erase part of the matrix.
/// @param firstrow Position of the first row, should be from 0 to Matrix::Rows()-1.
/// @param lastrow Position behind the last row, should be from 0 to Matrix::Rows().
/// @param firstcol Position of the first column, should be from 0 to Matrix::Cols()-1.
/// @param lastcol Position behind the last column, should be from 0 to Matrix::Cols().
void RemoveSubset(enumerator firstrow, enumerator lastrow, enumerator firstcol, enumerator lastcol)
{
enumerator shiftrow = lastrow-firstrow;
enumerator shiftcol = lastcol-firstcol;
Matrix<Var> tmp(n-shiftrow, m-shiftcol);
for(enumerator k = 0; k < firstrow; ++k)
{
for(enumerator l = 0; l < firstcol; ++l)
tmp(k,l) = (*this)(k,l);
for(enumerator l = lastcol; l < m; ++l)
tmp(k,l-shiftcol) = (*this)(k,l);
}
for(enumerator k = lastrow; k < n; ++k)
{
for(enumerator l = 0; l < firstcol; ++l)
tmp(k-shiftrow,l) = (*this)(k,l);
for(enumerator l = lastcol; l < m; ++l)
tmp(k-shiftrow,l-shiftcol) = (*this)(k,l);
}
this->Swap(tmp);
}
///Exchange contents of two matrices.
void Swap(AbstractMatrix<Var> & b)
{
if( dynamic_cast<Matrix<Var,storage_type> *>(&b) != NULL )
{
Matrix<Var,storage_type> & bb = dynamic_cast<Matrix<Var,storage_type> &>(b);
space.swap(bb.space);
std::swap(n,bb.n);
std::swap(m,bb.m);
}
else AbstractMatrix<Var>::Swap(b);
}
/// Construct empty matrix.
Matrix() : space(), n(0), m(0) {}
/// Construct the matrix from provided array and sizes.
/// @param pspace Array of elements of the matrix, stored in row-wise format.
/// @param pn Number of rows.
/// @param pm Number of columns.
Matrix(const Var * pspace, enumerator pn, enumerator pm) : space(pspace,pspace+pn*pm), n(pn), m(pm) {}
/// Construct the matrix with the provided storage.
/// Could be used to wrap existing array.
/// @param pspace Storage of elements of the matrix, stored in row-wise format.
/// @param pn Number of rows.
/// @param pm Number of columns.
/// \todo Do we need reference for pspace or just pspace?
Matrix(storage_type pspace, enumerator pn, enumerator pm) : space(pspace), n(pn), m(pm) {}
/// Construct a matrix with provided sizes.
/// @param pn Number of rows.
/// @param pm Number of columns.
/// \warning The matrix does not necessery have zero entries.
Matrix(enumerator pn, enumerator pm) : space(pn*pm), n(pn), m(pm) {}
/// Copy matrix.
/// @param other Another matrix of the same type.
Matrix(const Matrix & other) : space(other.n*other.m), n(other.n), m(other.m)
{
for(enumerator i = 0; i < n*m; ++i)
space[i] = other.space[i];
}
/// Construct matrix from matrix of different type.
/// Uses assign function declared in inmost_expression.h.
/// Copies derivative information if possible.
/// @param other Another matrix of different type.
template<typename typeB>
Matrix(const Matrix<typeB> & other) : space(other.Cols()*other.Rows()), n(other.Rows()), m(other.Cols())
Matrix(const AbstractMatrix<typeB> & other) : space(other.Cols()*other.Rows()), n(other.Rows()), m(other.Cols())
{
for(enumerator i = 0; i < n; ++i)
for(enumerator j = 0; j < m; ++j)
(*this)(i,j) = get_value(other(i,j));
assign((*this)(i,j),other(i,j));
}
/// Delete matrix.
~Matrix() {}
/// Resize the matrix into different size.
/// @param nrows New number of rows.
/// @param ncols New number of columns.
void Resize(enumerator nrows, enumerator mcols)
{
if( space.size() != mcols*nrows )
......@@ -457,19 +842,27 @@ namespace INMOST
n = nrows;
m = mcols;
}
/// Assign matrix of the same type.
/// @param other Another matrix of the same type.
/// @return Reference to matrix.
Matrix & operator =(Matrix const & other)
{
if( n*m != other.n*other.m ) space.resize(other.n*other.m);
if( n*m != other.n*other.m )
space.resize(other.n*other.m);
for(enumerator i = 0; i < other.n*other.m; ++i)
space[i] = other.space[i];
n = other.n;
m = other.m;
return *this;
}
/// Assign matrix of another type.
/// @param other Another matrix of different type.
/// @return Reference to matrix.
template<typename typeB>
Matrix & operator =(Matrix<typeB> const & other)
Matrix & operator =(AbstractMatrix<typeB> const & other)
{
if( Cols()*Rows() != other.Cols()*other.Rows() ) space.resize(other.Cols()*other.Rows());
if( Cols()*Rows() != other.Cols()*other.Rows() )
space.resize(other.Cols()*other.Rows());
for(enumerator i = 0; i < other.Rows(); ++i)
for(enumerator j = 0; j < other.Cols(); ++j)
assign((*this)(i,j),other(i,j));
......@@ -477,8 +870,10 @@ namespace INMOST
m = other.Cols();
return *this;
}
// i is in [0,n] - row index
// j is in [0,m] - column index
/// Access element of the matrix by row and column indices.
/// @param i Column index.
/// @param j Row index.
/// @return Reference to element.
Var & operator()(enumerator i, enumerator j)
{