mat.h

#pragma once
#include "baseexception.h"
#include "basestring.h"
#include "symtensor.h"
#include "matrix.h"
#include "basememory.h" 
// mat.h
 
// Mat.H - Header file for the matrix & vector classes.
//
// History: 03-Jul-96 Dave Potyondy
// last mod. 11-Jul-96 DOP
// last mod. 22-Jul-96 DOP added scalar multiply
// last mod. 24-Jul-96 DOP inlined oper. () and arr_idx functions
// last mod. 09-Aug-96 DOP Added read/write member functions
// last mod. 03-May-11 JTD Adapted for use as FISH matrix primitive.
//
// =======================================================================
 
namespace itasca {
 /* Exception handler for Matrix Class
 */
 struct NegMatExcept : public Exception {
 NegMatExcept();
 };
 
 /* General matrix class for (m x n) matrices */
 class BASE_EXPORT Mat {
 public:
 using size_type = uint32;
 
 Mat(); //default constructor
 Mat(size_type m, size_type n); //constructor of (m x n) matrix
 
 Mat(const Mat &mtx); //copy constructor
 Mat(Mat &&mtx); // rvalue-reference copy constructor
 
 Mat(const SymTensor &t);
 Mat(const DVect2 &v);
 Mat(const DVect3 &v);
 template <unsigned SX, unsigned SY>
 Mat(const Matrix<double, SX, SY> &v);
 Mat(const DMatrix<2, 2> &v);
 Mat(const DMatrix<3, 3> &v);
 
 virtual ~Mat(); //destructor
 
 // full access arr(i,j)
 double &operator()(size_type i, size_type j) { return arr[arr_idx(i, j)]; }
 // read access arr(i.j)
 const double &operator()(size_type i, size_type j) const { return arr[arr_idx(i, j)]; }
 
 Mat &operator =(const Mat &mtx); //assignment
 Mat &operator= (Mat &&mtx);
 
 Mat operator +(const Mat &mtx) const; //matrix addition
 Mat operator -(const Mat &mtx) const; //matrix subtraction
 Mat operator *(const Mat &mtx) const; //matrix multiply
 Mat operator *(double scal) const; //scalar multiply
 Mat operator* (const DVect2 &v) const; // 2d vector multiply
 Mat operator* (const DVect3 &v) const; // 3d vector multiply
 void operator*=(double s);
 void operator/=(double s);
 void operator+=(const Mat &mtx);
 void operator-=(const Mat &mtx);
 
 //fill matrix with given value
 void fill(double val) { std::fill(arr, arr + len, val); }
 void zero() { fill(0.0); }
 void identity(); // set matrix to identity matrix
 Mat transpose() const; // return matrix transpose
 // multiply by given scalar
 void scalMult(double scal) { auto *e = arr + len; for (auto *v = arr; v < e; ++v) *v *= scal; }
 
 bool equals(const Mat &mtx) const; // are the matrices equal?
 bool exactEquals(const Mat &mtx) const;
 bool operator== (const Mat &mtx) const { return equals(mtx); }
 virtual bool symmetric() const; // is the matrix symmetric?
 virtual double maxNorm() const; // return infinity-norm: max abs(elem)
 UVect2 size() const { UVect2 v(msize, nsize); return v; }
 
 UVect2 blockSize() const { UVect2 v(blk_m, blk_n); return v; }
 void setBlockSize(size_type blk_msize, size_type blk_nsize) { blk_m = blk_msize; blk_n = blk_nsize; }
 void addBlock(const Mat &src, //source matrix
 size_type src_bi, size_type src_bj, //block index of source
 size_type dst_bi, size_type dst_bj); //block index of dest.
 virtual void addGenBlock(const Mat &src, //source matrix
 size_type src_i, size_type src_j, //u.l.-corner of source
 size_type dst_i, size_type dst_j); //u.l.-corner of dest.
 
 SymTensor toTensor() const;
 DVect2 toVect2() const;
 DVect3 toVect3() const;
 template <unsigned SX, unsigned SY>
 Matrix<double, SX, SY> toMatrix() const;
 double *data() const { return arr; }
 protected:
 double *arr = nullptr; //store in row-major order as C-array -- (1,2,3,4) - (1,2)
 size_type msize, nsize; //matrix is (msize x nsize) (3,4)
 size_type len; //number of stored entries
 size_type blk_m = 1, blk_n = 1; //block size (m x n)
 
 bool same(const double *v1, const double *v2) const;
 inline size_type arr_idx(size_type i, size_type j) const;
 };
 
 template <unsigned SX, unsigned SY>
 Mat::Mat(const Matrix<double, SX, SY> &m) : arr(nullptr) {
 msize = SX;
 nsize = SY;
 len = SX * SY;
 arr = NEWC(double[SX * SY]);
 for (uint32 i = 0; i < SX; ++i)
 for (uint32 j = 0; j < SY; ++j)
 operator()(i, j) = m(i, j);
 }
 
 template <unsigned SX, unsigned SY>
 Matrix<double, SX, SY> Mat::toMatrix() const {
 Matrix<double, SX, SY> ret;
 if (msize != SX || nsize != SY) {
 if (msize == 1 && SY == 1 && nsize == SX) { // Allow generalized vectors to switch
 for (uint32 i = 0; i < SX; ++i)
 ret(i, 0) = operator()(0, i);
 return ret;
 }
 if (nsize == 1 && SX == 1 && msize == SY) { // Allow 
 for (uint32 i = 0; i < SY; ++i)
 ret(0, i) = operator()(i, 0);
 return ret;
 }
 throw Exception("Matrix size does not match {},{}, is {},{} instead.", SX, SY, msize, nsize);
 }
 for (uint32 i = 0; i < SX; ++i)
 for (uint32 j = 0; j < SY; ++j)
 ret(i, j) = operator()(i, j);
 return ret;
 }
 
 // =========================================
 // Returns index into the array.
 // Assumes indexing starts at (0,0).
 Mat::size_type Mat::arr_idx(size_type i, size_type j) const {
#ifdef _DEBUG
 if ((i >= msize) || (j >= nsize))
 throw Exception("Matrix ({} x {}) : index ({},{}) is out of bounds.", msize, nsize, i, j);
#endif
 return nsize * (i)+j;
 }
 
}
 
namespace base {
 using itasca::Mat;
 BASE_EXPORT string ts(const itasca::Mat &m, int width = 0, char notation = '\0', int precision = -1, char fill = ' ');
}
// EOF Mat.H