Itasca C++ Interface
mat.h
1 #pragma once
2 #include "baseexception.h"
3 #include "basestring.h"
4 #include "symtensor.h"
5 #include "matrix.h"
6 #include "basememory.h"
7 // mat.h
8 
9 // Mat.H - Header file for the matrix & vector classes.
10 //
11 // History: 03-Jul-96 Dave Potyondy
12 // last mod. 11-Jul-96 DOP
13 // last mod. 22-Jul-96 DOP added scalar multiply
14 // last mod. 24-Jul-96 DOP inlined oper. () and arr_idx functions
15 // last mod. 09-Aug-96 DOP Added read/write member functions
16 // last mod. 03-May-11 JTD Adapted for use as FISH matrix primitive.
17 //
18 // =======================================================================
19 
20 namespace itasca {
21  /* Exception handler for Matrix Class
22  */
23  struct NegMatExcept : public Exception {
24  NegMatExcept();
25  };
26 
27  /* General matrix class for (m x n) matrices */
28  class BASE_EXPORT Mat {
29  public:
30  using size_type = uint32;
31 
32  Mat(); //default constructor
33  Mat(size_type m, size_type n); //constructor of (m x n) matrix
34 
35  Mat(const Mat &mtx); //copy constructor
36  Mat(Mat &&mtx); // rvalue-reference copy constructor
37 
38  Mat(const SymTensor &t);
39  Mat(const DVect2 &v);
40  Mat(const DVect3 &v);
41  template <unsigned SX, unsigned SY>
42  Mat(const Matrix<double, SX, SY> &v);
43  Mat(const DMatrix<2, 2> &v);
44  Mat(const DMatrix<3, 3> &v);
45 
46  virtual ~Mat(); //destructor
47 
48  // full access arr(i,j)
49  double &operator()(size_type i, size_type j) { return arr[arr_idx(i, j)]; }
50  // read access arr(i.j)
51  const double &operator()(size_type i, size_type j) const { return arr[arr_idx(i, j)]; }
52 
53  Mat &operator =(const Mat &mtx); //assignment
54  Mat &operator= (Mat &&mtx);
55 
56  Mat operator +(const Mat &mtx) const; //matrix addition
57  Mat operator -(const Mat &mtx) const; //matrix subtraction
58  Mat operator *(const Mat &mtx) const; //matrix multiply
59  Mat operator *(double scal) const; //scalar multiply
60  Mat operator* (const DVect2 &v) const; // 2d vector multiply
61  Mat operator* (const DVect3 &v) const; // 3d vector multiply
62  void operator*=(double s);
63  void operator/=(double s);
64  void operator+=(const Mat &mtx);
65  void operator-=(const Mat &mtx);
66 
67  //fill matrix with given value
68  void fill(double val) { std::fill(arr, arr + len, val); }
69  void zero() { fill(0.0); }
70  void identity(); // set matrix to identity matrix
71  Mat transpose() const; // return matrix transpose
72  // multiply by given scalar
73  void scalMult(double scal) { auto *e = arr + len; for (auto *v = arr; v < e; ++v) *v *= scal; }
74 
75  bool equals(const Mat &mtx) const; // are the matrices equal?
76  bool exactEquals(const Mat &mtx) const;
77  bool operator== (const Mat &mtx) const { return equals(mtx); }
78  virtual bool symmetric() const; // is the matrix symmetric?
79  virtual double maxNorm() const; // return infinity-norm: max abs(elem)
80  UVect2 size() const { UVect2 v(msize, nsize); return v; }
81 
82  UVect2 blockSize() const { UVect2 v(blk_m, blk_n); return v; }
83  void setBlockSize(size_type blk_msize, size_type blk_nsize) { blk_m = blk_msize; blk_n = blk_nsize; }
84  void addBlock(const Mat &src, //source matrix
85  size_type src_bi, size_type src_bj, //block index of source
86  size_type dst_bi, size_type dst_bj); //block index of dest.
87  virtual void addGenBlock(const Mat &src, //source matrix
88  size_type src_i, size_type src_j, //u.l.-corner of source
89  size_type dst_i, size_type dst_j); //u.l.-corner of dest.
90 
91  SymTensor toTensor() const;
92  DVect2 toVect2() const;
93  DVect3 toVect3() const;
94  template <unsigned SX, unsigned SY>
95  Matrix<double, SX, SY> toMatrix() const;
96  double *data() const { return arr; }
97  protected:
98  double *arr = nullptr; //store in row-major order as C-array -- (1,2,3,4) - (1,2)
99  size_type msize, nsize; //matrix is (msize x nsize) (3,4)
100  size_type len; //number of stored entries
101  size_type blk_m = 1, blk_n = 1; //block size (m x n)
102 
103  bool same(const double *v1, const double *v2) const;
104  inline size_type arr_idx(size_type i, size_type j) const;
105  };
106 
107  template <unsigned SX, unsigned SY>
108  Mat::Mat(const Matrix<double, SX, SY> &m) : arr(nullptr) {
109  msize = SX;
110  nsize = SY;
111  len = SX * SY;
112  arr = NEWC(double[SX * SY]);
113  for (uint32 i = 0; i < SX; ++i)
114  for (uint32 j = 0; j < SY; ++j)
115  operator()(i, j) = m(i, j);
116  }
117 
118  template <unsigned SX, unsigned SY>
119  Matrix<double, SX, SY> Mat::toMatrix() const {
121  if (msize != SX || nsize != SY) {
122  if (msize == 1 && SY == 1 && nsize == SX) { // Allow generalized vectors to switch
123  for (uint32 i = 0; i < SX; ++i)
124  ret(i, 0) = operator()(0, i);
125  return ret;
126  }
127  if (nsize == 1 && SX == 1 && msize == SY) { // Allow
128  for (uint32 i = 0; i < SY; ++i)
129  ret(0, i) = operator()(i, 0);
130  return ret;
131  }
132  throw Exception("Matrix size does not match {},{}, is {},{} instead.", SX, SY, msize, nsize);
133  }
134  for (uint32 i = 0; i < SX; ++i)
135  for (uint32 j = 0; j < SY; ++j)
136  ret(i, j) = operator()(i, j);
137  return ret;
138  }
139 
140  // =========================================
141  // Returns index into the array.
142  // Assumes indexing starts at (0,0).
143  Mat::size_type Mat::arr_idx(size_type i, size_type j) const {
144 #ifdef _DEBUG
145  if ((i >= msize) || (j >= nsize))
146  throw Exception("Matrix ({} x {}) : index ({},{}) is out of bounds.", msize, nsize, i, j);
147 #endif
148  return nsize * (i)+j;
149  }
150 
151 }
152 
153 namespace base {
154  BASE_EXPORT string ts(const itasca::Mat &m, int width = 0, char notation = '\0', int precision = -1, char fill = ' ');
155 }
156 // EOF Mat.H
Comment point for memory allocation in all modules.
QString helper functions, plus some additions.
DMatrix is a Matrix that defaults to type double...
Definition: matrix.h:741
Base exception class for all Itasca code.
Definition: baseexception.h:10
A template-based matrix class, size fixed at compile time. Defaults to symmetric sized matrix.
Definition: matrix.h:22
Definition: mat.h:28
constexpr const Vector2< T > & toVect2(const Vector2< T > &v)
Conversion between vectors of different dimension.
Definition: vect.h:339
constexpr Vector3< T > toVect3(const Vector2< T > &v, const T &t=0)
Conversion between vectors of different dimension.
Definition: vect.h:341
#define BASE_EXPORT
Definition: basedef.h:24
A template-based matrix class, size fixed at compile time.
namespace Itasca
Definition: basememory.cpp:10
Definition: mat.h:23
A Symmetric 2nd order tensor.