Itasca C++ Interface
matrix.h
Go to the documentation of this file.
1 #pragma once
8 #ifdef __LINUX
9 #include <string.h>
10 #endif
11 
12 #include "symtensor.h"
13 #include <cassert>
14 
21 template <class T,unsigned SX,unsigned SY=SX>
22 class Matrix {
23 private:
25  class Column {
26  public:
28  constexpr Column(unsigned x,Matrix<T,SX,SY> &m) : x_(x), m_(m) { }
30  constexpr Column(const Column &r) : x_(r.x_), m_(r.m_) { }
32  constexpr const T &operator[](unsigned y) const { return m_.get(x_,y); }
34  T &operator[](unsigned y) { return m_.get(x_,y); }
35  private:
36  void operator=(const Column &r)=delete;
37  unsigned x_;
38  Matrix<T,SX,SY> &m_;
39  };
40  // Using template-meta programming to unroll innermost loop of matrix multiply
41  template <unsigned SZ,unsigned I,unsigned J,unsigned K>
42  class innerLoop {
43  public:
44  static double execute(const Matrix<T,SX,SY> &l,const Matrix<T,SY,SZ> &m) {
45  return l.t_[SX-I][SY-K]*m.t_[SY-K][SZ-J] + innerLoop<SZ,I,J,K-1>::execute(l,m);
46  }
47  };
48  template <unsigned SZ,unsigned I,unsigned J> // Termination specialization
49  class innerLoop<SZ,I,J,1> {
50  public:
51  static double execute(const Matrix<T,SX,SY> &l,const Matrix<T,SY,SZ> &m) {
52  return l.t_[SX-I][SY-1]*m.t_[SY-1][SZ-J];
53  }
54  };
55  template <unsigned I,unsigned K> // Specialization for column matrix (SZ=1)
56  class innerLoop<1,I,1,K> {
57  public:
58  static double execute(const Matrix<T,SX,SY> &l,const Matrix<T,SY,1> &m) {
59  return l.t_[SX-I][SY-K]*m.t_[SY-K] + innerLoop<1,I,1,K-1>::execute(l,m);
60  }
61  };
62  template <unsigned I> // termination specialization for column matrix specialization
63  class innerLoop<1,I,1,1> {
64  public:
65  static double execute(const Matrix<T,SX,SY> &l,const Matrix<T,SY,1> &m) {
66  return l.t_[SX-I][SY-1]*m.t_[SY-1];
67  }
68  };
69  // Using template meta-programming to unroll second loop of matrix multiply.
70  template <unsigned SZ,unsigned I,unsigned J>
71  class loop2Multiply {
72  public:
73  static void execute(const Matrix<T,SX,SY> &l,const Matrix<T,SY,SZ> &m,Matrix<T,SX,SZ> &ret) {
74  ret.t_[SX-I][SZ-J] = innerLoop<SZ,I,J,SY>::execute(l,m);
75  loop2Multiply<SZ,I,J-1>::execute(l,m,ret);
76  }
77  };
78  template <unsigned SZ,unsigned I> // Termination specialization
79  class loop2Multiply<SZ,I,1> {
80  public:
81  static void execute(const Matrix<T,SX,SY> &l,const Matrix<T,SY,SZ> &m,Matrix<T,SX,SZ> &ret) {
82  ret.t_[SX-I][SZ-1] = innerLoop<SZ,I,1,SY>::execute(l,m);
83  }
84  };
85  template <unsigned I> // Column matrix specialization (SZ=1)
86  class loop2Multiply<1,I,1> {
87  public:
88  static void execute(const Matrix<T,SX,SY> &l,const Matrix<T,SY,1> &m,Matrix<T,SX,1> &ret) {
89  ret.t_[SX-I] = innerLoop<1,I,1,SY>::execute(l,m);
90  }
91  };
92  // Using template meta-programming to unroll outermost loop of matrix multiply.
93  template <unsigned SZ,unsigned I>
94  class loopMultiply {
95  public:
96  static void execute(const Matrix<T,SX,SY> &l,const Matrix<T,SY,SZ> &m,Matrix<T,SX,SZ> &ret) {
97  loop2Multiply<SZ,I,SZ>::execute(l,m,ret);
98  loopMultiply<SZ,I-1>::execute(l,m,ret);
99  }
100  };
101  template <unsigned SZ> // Termination specialization
102  class loopMultiply<SZ,1> {
103  public:
104  static void execute(const Matrix<T,SX,SY> &l,const Matrix<T,SY,SZ> &m,Matrix<T,SX,SZ> &ret) {
105  loop2Multiply<SZ,1,SZ>::execute(l,m,ret);
106  }
107  };
108 
109 public:
111 #ifdef _DEBUG
112  Matrix() { set(initVal<T>()); }
113 #else
114 PUSHWARNING
115 VSWARNING(26495) // Initialization warning - on purpose
116  Matrix() {}
117 POPWARNING
118 #endif
120  explicit constexpr Matrix(const T &t) { set(t); }
121  explicit constexpr Matrix(const std::array<std::array<T,SY>,SX> &a) { for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) t_[i][j] = a[i][j]; }
122  constexpr Matrix(const std::initializer_list<std::initializer_list<T>> l);
124  constexpr Matrix(const Matrix<T,SX,SY> &m) { operator=(m); }
126  constexpr const Matrix<T,SX,SY> &operator=(const Matrix<T,SX,SY> &m) { for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) t_[i][j] = m.t_[i][j]; return *this; }
127 
129  constexpr const T & get(unsigned x,unsigned y) const { assert(x<SX); assert(y<SY); return t_[x][y]; }
131  constexpr T & get(unsigned x,unsigned y) { assert(x<SX); assert(y<SY); return t_[x][y]; }
133  constexpr const Column operator[](unsigned y) const { return Column(y,*this); }
135  constexpr Column operator[](unsigned y) { return Column(y,*this); }
137  constexpr const T & operator()(unsigned x,unsigned y) const { return get(x,y); }
139  constexpr T & operator()(unsigned x,unsigned y) { return get(x,y); }
140 
142  constexpr const Matrix<T,SX,SY> &operator+=(const Matrix<T,SX,SY> &m) { for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) t_[i][j] += m.t_[i][j]; return *this; }
144  constexpr const Matrix<T,SX,SY> &operator-=(const Matrix<T,SX,SY> &m) { for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) t_[i][j] -= m.t_[i][j]; return *this; }
146  constexpr const Matrix<T,SX,SY> &operator*=(const T &t) { for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) t_[i][j] *= t; return *this; }
148  constexpr const Matrix<T,SX,SY> &operator/=(const T &t) { for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) t_[i][j] /= t; return *this; }
149 
151  Matrix<T,SX,SY> operator+(const Matrix<T,SX,SY> &m) const { Matrix<T,SX,SY> ret; for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) ret.t_[i][j] = t_[i][j] + m.t_[i][j]; return ret; }
153  Matrix<T,SX,SY> operator-(const Matrix<T,SX,SY> &m) const { Matrix<T,SX,SY> ret; for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) ret.t_[i][j] = t_[i][j] - m.t_[i][j]; return ret; }
155  Matrix<T,SX,SY> operator*(const T &t) const { Matrix<T,SX,SY> ret; for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) ret.t_[i][j] = t_[i][j] * t; return ret; }
157  Matrix<T,SX,SY> operator/(const T &t) const { Matrix<T,SX,SY> ret; for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) ret.t_[i][j] = t_[i][j] / t; return ret; }
159  template <unsigned SZ>
161  Matrix<T,SX,SZ> ret;
162  // OK - we are using template expansion to do loop unrolling at compile time. Fun!
163  //for (unsigned i=0;i<SX;++i) { // ROWS OF FIRST
164  // for (unsigned j=0;j<SZ;++j) { // COLUMNS OF SECOND
165  // T &val = ret.get(i,j);
166  // val = 0;
167  // for (unsigned k=0;k<SY;++k) // Add up
168  // val += get(i,k)*m.get(k,j);
169  // }
170  //}
171  loopMultiply<SZ,SX>::execute(*this,m,ret);
172  return ret;
173  }
175  //template <class T>
176  //Matrix<T,SX,1> operator*(const Matrix<T,SY,1> &m) const {
177  // Matrix<T,SX,1> ret;
178  // // OK - we are using template expansion to do loop unrolling at compile time. Fun!
179  // //for (unsigned i=0;i<SX;++i) { // ROWS OF FIRST
180  // // T &val = ret.get(i,0);
181  // // val = 0.0;
182  // // for (unsigned k=0;k<SY;++k) // Add up
183  // // val += get(i,k)*m.get(k,0);
184  // //}
185  // loopMultiply<1,SX>::execute(*this,m,ret);
186  // return ret;
187  //}
188 
189  // Block operations
191  template <unsigned BX,unsigned BY,unsigned SX2,unsigned SY2>
192  void addBlock(const Matrix<T,SX2,SY2> &m,unsigned iSrc,unsigned jSrc,unsigned iDst,unsigned jDst) {
193  addGenBlock<BX,BY,SX2,SY2>(m,iSrc*BX,jSrc*BY,iDst*BX,jDst*BY);
194  }
196  template <unsigned BX,unsigned BY,unsigned SX2,unsigned SY2>
197  void addGenBlock(const Matrix<T,SX2,SY2> &m,unsigned iSrc,unsigned jSrc,unsigned iDst,unsigned jDst) {
198  for (unsigned i=iSrc,id=iDst;i<(iSrc+BX);i++,id++)
199  for (unsigned j=jSrc,jd=jDst;j<(jSrc+BY);j++,jd++)
200  get(id,jd) += m(i,j);
201  }
202 
204  T maxNorm() const { T ret(0); for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) ret = std::max(ret,std::abs(t_[i][j])); return ret; }
206  Matrix<T,SY,SX> transpose() const { Matrix<T,SY,SX> ret; for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) ret(j,i) = t_[i][j]; return ret; }
208  T trace() const { T tt = t_[0][0]; unsigned len = std::min(SX,SY); for (unsigned i=1; i<len; ++i) tt += t_[i][i]; return tt; }
210  constexpr void set(const T &t=T()) { for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SY;++j) t_[i][j] = t; }
211 
213  static constexpr Matrix<T,SX,SY> identity() { Matrix<T,SX,SY> ret(0.0); for (unsigned i=0;i<std::min(SX,SY);++i) ret(i,i) = 1.0; return ret; }
214 
215  T t_[SX][SY];
216 };
217 
218 template <class T,unsigned SX,unsigned SY>
219 constexpr Matrix<T,SX,SY>::Matrix(const std::initializer_list<std::initializer_list<T>> l) {
220  uint32 i=0;
221  for (auto &y : l) {
222  uint32 j = 0;
223  for (auto &x : y) {
224  t_[i][j] = x;
225  ++j;
226  }
227  ++i;
228  }
229 }
230 
231 
233 template <class T,unsigned SX>
234 class Matrix<T,SX,1> {
235 public:
237 #ifdef _DEBUG
238  Matrix() { set(initVal<T>()); }
239 #else
240  Matrix() { }
241 #endif
243  constexpr explicit Matrix(const T &t) { set(t); }
245  constexpr Matrix(const Matrix<T,SX,1> &m) { operator=(m); }
246  constexpr Matrix(const std::array<T,SX> &a) { for (uint32 i=0;i<SX;++i) t_[i] = a[i]; }
247  constexpr Matrix(const std::initializer_list<T> l) { uint32 i=0; for (auto &x : l) { t_[i] = x; ++i; } }
249  constexpr const Matrix<T,SX,1> &operator=(const Matrix<T,SX,1> &m) { for (uint32 i=0;i<SX;++i) t_[i] = m.t_[i]; return *this; }
250 
252  constexpr const T & get(unsigned x,[[maybe_unused]] unsigned y) const { assert(x<SX); assert(!y); return t_[x]; }
254  constexpr const T & get(unsigned x) const { assert(x<SX); return t_[x]; }
256  constexpr T & get(unsigned x,[[maybe_unused]] unsigned y) { assert(x<SX); assert(!y); return t_[x]; }
258  constexpr T & get(unsigned x) { assert(x<SX); return t_[x]; }
260  constexpr const T & operator()(unsigned x,unsigned y) const { return get(x,y); }
262  constexpr const T & operator()(unsigned x) const { return get(x); }
264  constexpr T & operator()(unsigned x,unsigned y) { return get(x,y); }
266  constexpr T & operator()(unsigned x) { return get(x); }
267 
269  constexpr const Matrix<T,SX,1> &operator+=(const Matrix<T,SX,1> &m) { for (unsigned i=0;i<SX;++i) t_[i] += m.t_[i]; return *this; }
271  constexpr const Matrix<T,SX,1> &operator-=(const Matrix<T,SX,1> &m) { for (unsigned i=0;i<SX;++i) t_[i] -= m.t_[i]; return *this; }
273  constexpr const Matrix<T,SX,1> &operator*=(const T &t) { for (unsigned i=0;i<SX;++i) t_[i] *= t; return *this; }
275  constexpr const Matrix<T,SX,1> &operator/=(const T &t) { for (unsigned i=0;i<SX;++i) t_[i] /= t; return *this; }
276 
278  Matrix<T,SX,1> operator+(const Matrix<T,SX,1> &m) const { Matrix<T,SX,1> ret; for (unsigned i=0;i<SX;++i) ret.t_[i] = t_[i] + m.t_[i]; return ret; }
280  Matrix<T,SX,1> operator-(const Matrix<T,SX,1> &m) const { Matrix<T,SX,1> ret; for (unsigned i=0;i<SX;++i) ret.t_[i] = t_[i] - m.t_[i]; return ret; }
282  Matrix<T,SX,1> operator*(const T &t) const { Matrix<T,SX,1> ret; for (unsigned i=0;i<SX;++i) ret.t_[i] = t_[i] * t; return ret; }
284  Matrix<T,SX,1> operator/(const T &t) const { Matrix<T,SX,1> ret; for (unsigned i=0;i<SX;++i) ret.t_[i] = t_[i] / t; return ret; }
286  template <unsigned SZ> Matrix<T,SX,SZ> operator*(const Matrix<T,1,SZ> &m) const {
287  Matrix<T,SX,SZ> ret;
288  for (unsigned i=0;i<SX;++i) // ROWS OF FIRST
289  for (unsigned j=0;j<SZ;++j) // COLUMNS OF SECOND
290  ret.get(i,j) = get(i,0)*m.get(0,j);
291  return ret;
292  }
293 
294  // Block operations
296  template <unsigned BX,unsigned BY,unsigned SX2,unsigned SY2> void addBlock(const Matrix<T,SX2,SY2> &m,unsigned iSrc,unsigned jSrc,unsigned iDst,unsigned jDst) {
297  addGenBlock<BX,BY,SX2,SY2>(m,iSrc*BX,jSrc*BY,iDst*BX,jDst*BY);
298  }
300  template <unsigned BX,unsigned SX2> void addBlock(const Matrix<T,SX2,1> &m,unsigned iSrc,unsigned iDst) {
301  addGenBlock<BX,SX2>(m,iSrc*BX,iDst*BX);
302  }
304  template <unsigned BX,unsigned BY,unsigned SX2,unsigned SY2> void addGenBlock(const Matrix<T,SX2,SY2> &m,unsigned iSrc,unsigned jSrc,unsigned iDst,unsigned jDst) {
305  for (unsigned i=iSrc,id=iDst;i<(iSrc+BX);i++,id++)
306  for (unsigned j=jSrc,jd=jDst;j<(jSrc+BY);j++,jd++)
307  get(id,jd) += m(i,j);
308  }
310  template <unsigned BX,unsigned SX2> void addGenBlock(const Matrix<T,SX2,1> &m,unsigned iSrc,unsigned iDst) {
311  for (unsigned i=iSrc,id=iDst;i<(iSrc+BX);i++,id++)
312  get(id) += m(i);
313  }
314 
316  T maxNorm() const { T ret(0); for (unsigned i=0;i<SX;++i) ret = std::max(ret,std::abs(t_[i])); return ret; }
318  Matrix<T,1,SX> transpose() const { Matrix<T,1,SX> ret; for (unsigned i=0;i<SX;++i) ret(0,i) = t_[i]; return ret; }
320  void set(const T &t=T()) { for (unsigned i=0;i<SX;++i) t_[i] = t; }
321 
323  static Matrix<T,SX,1> identity() { Matrix<T,SX,1> ret(0.0); ret.t_[0] = 1.0; return ret; }
324  T t_[SX];
325 };
326 
331 template <class T,unsigned SX>
332 class SymMatrix {
333 private:
335  class Column {
336  public:
338  Column(unsigned x,SymMatrix<T,SX> &m) : x_(x), m_(m) { }
340  Column(const Column &r) : x_(r.x_), m_(r.m_) { }
342  const T &operator[](unsigned y) const { return m_.get(x_,y); }
344  T &operator[](unsigned y) { return m_.get(x_,y); }
345  private:
346  void operator=(const Column &r); // DISALLOWED
347  unsigned x_;
348  SymMatrix<T,SX> &m_;
349  };
350  // Using template-meta programming to unroll innermost loop of matrix multiply
351  template <unsigned SZ,unsigned I,unsigned J,unsigned K>
352  class innerLoop {
353  public:
354  static double execute(const SymMatrix<T,SX> &l,const Matrix<T,SX,SZ> &m) {
355  return l.t_[index(SX-I,SX-K)]*m.t_[SX-K][SZ-J] + innerLoop<SZ,I,J,K-1>::execute(l,m);
356  }
357  };
358  template <unsigned SZ,unsigned I,unsigned J> // Termination specialization
359  class innerLoop<SZ,I,J,1> {
360  public:
361  static double execute(const SymMatrix<T,SX> &l,const Matrix<T,SX,SZ> &m) {
362  return l.t_[index(SX-I,SX-1)]*m.t_[SX-1][SZ-J];
363  }
364  };
365  template <unsigned I,unsigned K> // Specialization for column matrix (SZ=1)
366  class innerLoop<1,I,1,K> {
367  public:
368  static double execute(const SymMatrix<T,SX> &l,const Matrix<T,SX,1> &m) {
369  return l.t_[index(SX-I,SX-K)]*m.t_[SX-K] + innerLoop<1,I,1,K-1>::execute(l,m);
370  }
371  };
372  template <unsigned I> // termination specialization for column matrix specialization
373  class innerLoop<1,I,1,1> {
374  public:
375  static double execute(const SymMatrix<T,SX> &l,const Matrix<T,SX,1> &m) {
376  return l.t_[index(SX-I,SX-1)]*m.t_[SX-1];
377  }
378  };
379  // Using template meta-programming to unroll second loop of matrix multiply.
380  template <unsigned SZ,unsigned I,unsigned J>
381  class loop2Multiply {
382  public:
383  static void execute(const SymMatrix<T,SX> &l,const Matrix<T,SX,SZ> &m,Matrix<T,SX,SZ> &ret) {
384  ret.t_[SX-I][SZ-J] = innerLoop<SZ,I,J,SX>::execute(l,m);
385  loop2Multiply<SZ,I,J-1>::execute(l,m,ret);
386  }
387  };
388  template <unsigned SZ,unsigned I> // Termination specialization
389  class loop2Multiply<SZ,I,1> {
390  public:
391  static void execute(const SymMatrix<T,SX> &l,const Matrix<T,SX,SZ> &m,Matrix<T,SX,SZ> &ret) {
392  ret.t_[SX-I][SZ-1] = innerLoop<SZ,I,1,SX>::execute(l,m);
393  }
394  };
395  template <unsigned I> // Column matrix specialization (SZ=1)
396  class loop2Multiply<1,I,1> {
397  public:
398  static void execute(const SymMatrix<T,SX> &l,const Matrix<T,SX,1> &m,Matrix<T,SX,1> &ret) {
399  ret.t_[SX-I] = innerLoop<1,I,1,SX>::execute(l,m);
400  }
401  };
402  // Using template meta-programming to unroll outermost loop of matrix multiply.
403  template <unsigned SZ,unsigned I>
404  class loopMultiply {
405  public:
406  static void execute(const SymMatrix<T,SX> &l,const Matrix<T,SX,SZ> &m,Matrix<T,SX,SZ> &ret) {
407  loop2Multiply<SZ,I,SZ>::execute(l,m,ret);
408  loopMultiply<SZ,I-1>::execute(l,m,ret);
409  }
410  };
411  template <unsigned SZ> // Termination specialization
412  class loopMultiply<SZ,1> {
413  public:
414  static void execute(const SymMatrix<T,SX> &l,const Matrix<T,SX,SZ> &m,Matrix<T,SX,SZ> &ret) {
415  loop2Multiply<SZ,1,SZ>::execute(l,m,ret);
416  }
417  };
418 
419 
420  // Using template-meta programming to unroll innermost loop of matrix multiply
421  template <unsigned I,unsigned J,unsigned K>
422  class innerLoopS {
423  public:
424  static double execute(const SymMatrix<T,SX> &l,const SymMatrix<T,SX> &m) {
425  return l.t_[index(SX-I,SX-K)]*m.t_[index(SX-K,SX-J)] + innerLoopS<I,J,K-1>::execute(l,m);
426  }
427  };
428 
429  template <unsigned I,unsigned J> // Termination specialization
430  class innerLoopS<I,J,1> {
431  public:
432  static double execute(const SymMatrix<T,SX> &l,const SymMatrix<T,SX> &m) {
433  return l.t_[index(SX-I,SX-1)]*m.t_[index(SX-1,SX-J)];
434  // what is X here? seems undefined??
435  }
436  };
437 
438  // Using template meta-programming to unroll second loop of matrix multiply.
439  template <unsigned I,unsigned J>
440  class loop2MultiplyS {
441  public:
442  static void execute(const SymMatrix<T,SX> &l,const SymMatrix<T,SX> &m,Matrix<T,SX,SX> &ret) {
443  ret.t_[SX-I][SX-J] = innerLoopS<I,J,SX>::execute(l,m);
444  loop2MultiplyS<I,J-1>::execute(l,m,ret);
445  }
446  };
447  template <unsigned I> // Termination specialization
448  class loop2MultiplyS<I,1> {
449  public:
450  static void execute(const SymMatrix<T,SX> &l,const SymMatrix<T,SX> &m,Matrix<T,SX,SX> &ret) {
451  ret.t_[SX-I][SX-1] = innerLoopS<I,1,SX>::execute(l,m);
452  }
453  };
454  // Using template meta-programming to unroll outermost loop of matrix multiply.
455  template <unsigned I>
456  class loopMultiplyS {
457  public:
458  static void execute(const SymMatrix<T,SX> &l,const SymMatrix<T,SX> &m,Matrix<T,SX,SX> &ret) {
459  loop2MultiplyS<I,SX>::execute(l,m,ret);
460  loopMultiplyS<I-1>::execute(l,m,ret);
461  }
462  };
463 #ifndef __LINUX
464  template <> // Termination specialization
465  class loopMultiplyS<1> {
466  public:
467  static void execute(const SymMatrix<T,SX> &l,const SymMatrix<T,SX> &m,Matrix<T,SX,SX> &ret) {
468  loop2MultiplyS<1,SX>::execute(l,m,ret);
469  }
470  };
471 #endif
472 
473 public:
475 #ifdef _DEBUG
476  SymMatrix() { set(initVal<T>()); }
477 #else
478  SymMatrix() { }
479 #endif
481  explicit SymMatrix(const T &t) { set(t); }
483  SymMatrix(const SymMatrix<T,SX> &m) { memcpy(t_,m.t_,sizeof(T)*len_); }
485  const SymMatrix<T,SX> &operator=(const SymMatrix<T,SX> &m) { memcpy(t_,m.t_,sizeof(T)*len_); return *this; }
486 
488  const T & get(unsigned x,unsigned y) const { return t_[index(x,y)]; }
490  T & get(unsigned x,unsigned y) { return t_[index(x,y)]; }
492  const Column operator[](unsigned y) const { return Column(y,*this); }
494  Column operator[](unsigned y) { return Column(y,*this); }
496  const T & operator()(unsigned x,unsigned y) const { return get(x,y); }
498  T & operator()(unsigned x,unsigned y) { return get(x,y); }
499 
501  const SymMatrix<T,SX> &operator+=(const SymMatrix<T,SX> &m) { for (unsigned i=0;i<len_;++i) t_[i] += m.t_[i]; return *this; }
503  const SymMatrix<T,SX> &operator-=(const SymMatrix<T,SX> &m) { for (unsigned i=0;i<len_;++i) t_[i] -= m.t_[i]; return *this; }
505  const SymMatrix<T,SX> &operator*=(const T &t) { for (unsigned i=0;i<len_;++i) t_[i] *= t; return *this; }
507  const SymMatrix<T,SX> &operator/=(const T &t) { for (unsigned i=0;i<len_;++i) t_[i] /= t; return *this; }
508 
510  SymMatrix<T,SX> operator+(const SymMatrix<T,SX> &m) const { SymMatrix<T,SX> ret; for (unsigned i=0;i<len_;++i) ret.t_[i] = t_[i] + m.t_[i]; return ret; }
512  Matrix<T,SX,SX> operator+(const Matrix<T,SX,SX> &m) const { Matrix<T,SX,SX> ret; for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SX;++j) ret(i,j) = get(i,j) + m.get(i,j); return ret; }
514  SymMatrix<T,SX> operator-(const SymMatrix<T,SX> &m) const { SymMatrix<T,SX> ret; for (unsigned i=0;i<len_;++i) ret.t_[i] = t_[i] - m.t_[i]; return ret; }
516  Matrix<T, SX, SX> operator-(const Matrix<T, SX, SX>& m) const { Matrix<T, SX, SX> ret; for (unsigned i = 0; i < SX; ++i) for (unsigned j = 0; j < SX; ++j) ret(i, j) = get(i, j) - m.get(i, j); return ret; }
517 
519  SymMatrix<T,SX> operator*(const T &t) const { SymMatrix<T,SX> ret; for (unsigned i=0;i<len_;++i) ret.t_[i] = t_[i] * t; return ret; }
521  SymMatrix<T,SX> operator/(const T &t) const { SymMatrix<T,SX> ret; for (unsigned i=0;i<len_;++i) ret.t_[i] = t_[i] / t; return ret; }
524  Matrix<T,SX,SX> ret;
525  loopMultiplyS<SX>::execute(*this,m,ret);
526  //for (unsigned i=0;i<SX;++i) { // ROWS OF FIRST
527  // for (unsigned j=0;j<SX;++j) { // COLUMNS OF SECOND
528  // T &val = ret.get(i,j);
529  // val = 0.0;
530  // for (unsigned k=0;k<SX;++k) // Add up
531  // val += get(i,k)*m.get(k,j);
532  // }
533  //}
534  return ret;
535  }
537  template <unsigned SZ> Matrix<T,SX,SZ> operator*(const Matrix<T,SX,SZ> &m) const {
538  Matrix<T,SX,SZ> ret;
539  loopMultiply<SZ,SX>::execute(*this,m,ret);
540  //for (unsigned i=0;i<SX;++i) { // ROWS OF FIRST
541  // for (unsigned j=0;j<SY;++j) { // COLUMNS OF SECOND
542  // T &val = ret.get(i,j);
543  // val = 0.0;
544  // for (unsigned k=0;k<SX;++k) // Add up
545  // val += get(i,k)*m.get(k,j);
546  // }
547  //}
548  return ret;
549  }
551  //template <> Matrix<T,SX,1> operator*(const Matrix<T,SX,1> &m) const {
552  // Matrix<T,SX,1> ret;
553  // for (unsigned i=0;i<SX;++i) { // ROWS OF FIRST
554  // T &val = ret.get(i,0);
555  // val = 0.0;
556  // for (unsigned k=0;k<SX;++k) // Add up
557  // val += get(i,k)*m.get(k,0);
558  // }
559  // return ret;
560  //}
561 
562  // Block operations
564  template <unsigned BX,unsigned BY,unsigned SX2,unsigned SY2> void addBlock(const Matrix<T,SX2,SY2> &m,unsigned iSrc,unsigned jSrc,unsigned iDst,unsigned jDst) {
565  addGenBlock<BX,BY,SX2,SY2>(m,iSrc*BX,jSrc*BY,iDst*BX,jDst*BY);
566  }
568  template <unsigned BX,unsigned BY,unsigned SX2> void addBlock(const SymMatrix<T,SX2> &m,unsigned iSrc,unsigned jSrc,unsigned iDst,unsigned jDst) {
569  addGenBlock<BX,BY,SX2>(m,iSrc*BX,jSrc*BY,iDst*BX,jDst*BY);
570  }
572  template <unsigned BX,unsigned BY,unsigned SX2,unsigned SY2> void addGenBlock(const Matrix<T,SX2,SY2> &m,unsigned iSrc,unsigned jSrc,unsigned iDst,unsigned jDst) {
573  for (unsigned i=iSrc,id=iDst;i<(iSrc+BX);i++,id++) {
574  for (unsigned j=jSrc,jd=jDst;j<(jSrc+BY);j++,jd++) {
575  if (id<=jd)
576  get(id,jd) += m(i,j);
577  }
578  }
579  }
581  template <unsigned BX,unsigned BY,unsigned SX2> void addGenBlock(const SymMatrix<T,SX2> &m,unsigned iSrc,unsigned jSrc,unsigned iDst,unsigned jDst) {
582  for (unsigned i=iSrc,id=iDst;i<(iSrc+BX);i++,id++) {
583  for (unsigned j=jSrc,jd=jDst;j<(jSrc+BY);j++,jd++) {
584  if (id<=jd)
585  get(id,jd) += m(i,j);
586  }
587  }
588  }
589 
591  T maxNorm() const { T ret(0); for (unsigned i=0;i<len_;++i) ret = std::max(std::abs(t_[i]),ret); return ret; }
593  SymMatrix<T,SX> transpose() const { return *this; }
595  T trace() const { T tt = t_[0]; for (unsigned i=1; i<SX; ++i) tt += get(i,i); return tt; }
597  void set(const T &t=T()) { for (unsigned i=0;i<len_;++i) t_[i] = t; }
599  Matrix<T,SX,SX> toMatrix() const { Matrix<T,SX,SX> ret; for (unsigned i=0;i<SX;++i) for (unsigned j=0;j<SX;++j) ret(i,j)=get(i,j); return ret; }
600 
602  static SymMatrix<T,SX> identity() { SymMatrix<T,SX> ret(0.0); for (unsigned i=0;i<SX;++i) ret(i,i) = 1.0; return ret; }
605  SymMatrix<T,SX> ret;
606 #ifdef _DEBUG
607  double maxdif = 0.0;
608 #endif
609  for (unsigned i=0;i<SX;++i) {
610  ret(i,i) = m(i,i);
611  for (unsigned j=i+1;j<SX;++j) {
612  ret(i,j) = (m(i,j) + m(j,i)) * 0.5;
613 #ifdef _DEBUG
614  maxdif = std::max(std::abs(ret(i,j) - m(i,j)),maxdif);
615 #endif
616  }
617  }
618 #ifdef _DEBUG
619  assert(maxdif <= ret.maxNorm()*limits<float>::epsilon());
620 #endif
621  return ret;
622  }
623 
624  static constexpr unsigned len_ = (SX*SX + SX) / 2;
625  static constexpr unsigned indexConst_ = 1 + 2*SX;
626  static constexpr unsigned index(unsigned x,unsigned y) {
627  unsigned x2 = x > y ? y : x;
628  unsigned y2 = x > y ? x : y;
629  unsigned i = (((indexConst_-x2)*x2)/2) + (y2-x2);
630  return i;
631  }
632 
633  T t_[len_];
634 };
635 
636 // MOO NOTE!:: There are optimizations that can be made here!
639 // Using template-meta programming to unroll innermost loop of matrix multiply
640 template <class T,unsigned SX,unsigned SY,unsigned I,unsigned J,unsigned K>
641 class innerLoopMS {
642 public:
643  static double execute(const Matrix<T,SX,SY> &l,const SymMatrix<T,SY> &m) {
644  return l.t_[SX-I][SY-K]*m.t_[m.index(SY-K,SY-J)] + innerLoopMS<T,SX,SY,I,J,K-1>::execute(l,m);
645  }
646 };
647 template <class T,unsigned SX,unsigned SY,unsigned I,unsigned J> // Termination specialization
648 class innerLoopMS<T,SX,SY,I,J,1> {
649 public:
650  static double execute(const Matrix<T,SX,SY> &l,const SymMatrix<T,SY> &m) {
651  return l.t_[SX-I][SY-1]*m.t_[m.index(SY-1,SY-J)];
652  }
653 };
654 // Using template meta-programming to unroll second loop of matrix multiply.
655 template <class T,unsigned SX,unsigned SY,unsigned I,unsigned J>
657 public:
658  static void execute(const Matrix<T,SX,SY> &l,const SymMatrix<T,SY> &m,Matrix<T,SX,SY> &ret) {
659  ret.t_[SX-I][SY-J] = innerLoopMS<T,SX,SY,I,J,SY>::execute(l,m);
661  }
662 };
663 template <class T,unsigned SX,unsigned SY,unsigned I> // Termination specialization
664 class loop2MultiplyMS<T,SX,SY,I,1> {
665 public:
666  static void execute(const Matrix<T,SX,SY> &l,const SymMatrix<T,SY> &m,Matrix<T,SX,SY> &ret) {
667  ret.t_[SX-I][SY-1] = innerLoopMS<T,SX,SY,I,1,SY>::execute(l,m);
668  }
669 };
670 // Using template meta-programming to unroll outermost loop of matrix multiply.
671 template <class T,unsigned SX,unsigned SY,unsigned I>
673 public:
674  static void execute(const Matrix<T,SX,SY> &l,const SymMatrix<T,SY> &m,Matrix<T,SX,SY> &ret) {
677  }
678 };
679 template <class T,unsigned SX,unsigned SY> // Termination specialization
680 class loopMultiplyMS<T,SX,SY,1> {
681 public:
682  static void execute(const Matrix<T,SX,SY> &l,const SymMatrix<T,SY> &m,Matrix<T,SX,SY> &ret) {
684  }
685 };
686 
687 template <class T,unsigned SX,unsigned SY>
689  Matrix<T,SX,SY> ret;
691  //for (unsigned i=0;i<SX;++i) { // ROWS OF FIRST
692  // for (unsigned j=0;j<SY;++j) { // COLUMNS OF SECOND
693  // T &val = ret.get(i,j);
694  // val = 0.0;
695  // for (unsigned k=0;k<SY;++k) // Add up
696  // val += m.get(i,k)*s.get(k,j);
697  // }
698  //}
699  return ret;
700 }
701 
707 template <class T,unsigned S> class VMatrix : public Matrix<T,S,1> {
708 public:
709  using Matrix<T,S,1>::Matrix;
712  //constexpr VMatrix() : Matrix<T,S,1>() { }
714  //explicit constexpr VMatrix(const T &t) : Matrix<T,S,1>(t) { }
716  constexpr VMatrix(const VMatrix<T,S> &m) : Matrix<T,S,1>(m) { }
718  constexpr VMatrix(const Matrix<T,S,1> &m) : Matrix<T,S,1>(m) { }
720  //constexpr const VMatrix<T,S> &operator=(const VMatrix<T,S> &m) { Matrix<T,S,1>::operator=(m); return *this; }
722  //constexpr const VMatrix<T,S> &operator=(const Matrix<T,S,1> &m) { Matrix<T,S,1>::operator=(m); return *this; }
723 
725  //const T & get(unsigned x) const { return Matrix<T,S,1>::get(x,0); }
727  //T & get(unsigned x) { return Matrix<T,S,1>::get(x,0); }
729  constexpr const T &operator[](unsigned x) const { return this->get(x); }
731  constexpr T & operator[](unsigned x) { return this->get(x); }
733  //const T & operator()(unsigned x) const { return get(x); }
735  //T & operator()(unsigned x) { return get(x); }
736 };
737 
741 template <unsigned SX,unsigned SY=SX> class DMatrix : public Matrix<double,SX,SY> {
742 public:
746  constexpr DMatrix() : Matrix<double,SX,SY>() { }
747  //explicit constexpr DMatrix(const double t[SX][SY]) : Matrix<double,SX,SY>(t) { }
749  //explicit constexpr DMatrix(const double &t) : Matrix<double,SX,SY>(t) { }
751  constexpr DMatrix(const DMatrix<SX,SY> &m) : Matrix<double,SX,SY>(m) { }
753  constexpr DMatrix(const Matrix<double,SX,SY> &m) : Matrix<double,SX,SY>(m) { }
755  constexpr const DMatrix<SX,SY> &operator=(const DMatrix<SX,SY> &m) { Matrix<double,SX,SY>::operator=(m); return *this; }
757  //constexpr const DMatrix<SX,SY> &operator=(const Matrix<double,SX,SY> &m) { Matrix<double,SX,SY>::operator=(m); return *this; }
758 };
759 
763 template <unsigned SX> class DSymMatrix : public SymMatrix<double,SX> {
764 public:
766  DSymMatrix() : SymMatrix<double,SX>() { }
768  explicit DSymMatrix(const double &t) : SymMatrix<double,SX>(t) { }
770  DSymMatrix(const DSymMatrix<SX> &m) : SymMatrix<double,SX>(m) { }
772  DSymMatrix(const SymMatrix<double,SX> &m) : SymMatrix<double,SX>(m) { }
777 };
778 
781 template <unsigned S> class DVMatrix : public VMatrix<double,S> {
782 public:
785 
787  //DVMatrix() : VMatrix<double,S>() { }
789  //explicit DVMatrix(const double &d) : VMatrix<double,S>(d) { }
791  constexpr DVMatrix(const DVMatrix<S> &m) : VMatrix<double,S>(m) { }
793  constexpr DVMatrix(const Matrix<double,S,1> &m) : VMatrix<double,S>(m) { }
795  //const DVMatrix<S> &operator=(const DVMatrix<S> &m) { VMatrix<double,S>::operator=(m); return *this; }
797  //const DVMatrix<S> &operator=(const Matrix<double,S,1> &m) { VMatrix<double,S>::operator=(m); return *this; }
798 };
799 
800 template <class T,unsigned S>
801 double innerProduct(const VMatrix<T,S> &v1,const VMatrix<T,S> &v2) {
802  double d=0.0;
803  for (uint32 i=0;i<S;++i)
804  d += v1(i)*v2(i);
805  return d;
806 }
807 
808 template <class T,unsigned S>
810  Matrix<T,S,S> ret;
811  for (uint32 i=0;i<S;++i) {
812  for (uint32 j=0;j<S;++j)
813  ret(i,j) = v1(i) * v2(j);
814  }
815  return ret;
816 }
817 
821 template <class T>
822 Matrix<T,3,3> outerProduct(const Vector3<T> &v1,const Vector3<T> &v2) {
823  Matrix<T,3,3> ret;
824  ret.get(0,0) = v1.x() * v2.x();
825  ret.get(0,1) = v1.x() * v2.y();
826  ret.get(0,2) = v1.x() * v2.z();
827  ret.get(1,0) = v1.y() * v2.x();
828  ret.get(1,1) = v1.y() * v2.y();
829  ret.get(1,2) = v1.y() * v2.z();
830  ret.get(2,0) = v1.z() * v2.x();
831  ret.get(2,1) = v1.z() * v2.y();
832  ret.get(2,2) = v1.z() * v2.z();
833  return ret;
834 }
835 
839 template <class T>
840 Matrix<T,2,2> outerProduct(const Vector2<T> &v1,const Vector2<T> &v2) {
841  Matrix<T,2,2> ret;
842  ret.get(0,0) = v1.x() * v2.x();
843  ret.get(0,1) = v1.x() * v2.y();
844  ret.get(1,0) = v1.y() * v2.x();
845  ret.get(1,1) = v1.y() * v2.y();
846  return ret;
847 }
848 
853 template <class T>
854 constexpr T determinant(const Matrix<T,3,3> &mat) {
855  T t = mat.get(0,0) * (mat.get(1,1) * mat.get(2,2) - mat.get(2,1) * mat.get(1,2));
856  t -= mat.get(1,0) * (mat.get(0,1) * mat.get(2,2) - mat.get(2,1) * mat.get(0,2));
857  t += mat.get(2,0) * (mat.get(0,1) * mat.get(1,2) - mat.get(1,1) * mat.get(0,2));
858  return t;
859 }
860 
864 template <class T>
865 constexpr T determinant(const Matrix<T,2,2> &mat) {
866  return mat.get(0,0) * mat.get(1,1) - mat.get(0,1) * mat.get(1,0);
867 }
868 
873 template <class T>
874 constexpr Matrix<T,3,3> inverse(const Matrix<T,3,3> &mat) {
875  // First calculate the adjoint
876  const T a11 = mat(0,0), a21 = mat(1,0), a31 = mat(2,0);
877  const T a12 = mat(0,1), a22 = mat(1,1), a32 = mat(2,1);
878  const T a13 = mat(0,2), a23 = mat(1,2), a33 = mat(2,2);
879 
880  Matrix<T,3,3> rm = {{ a22*a33 - a32*a23,-a12*a33 + a32*a13, a12*a23 - a22*a13},
881  {-a21*a33 + a31*a23, a11*a33 - a31*a13,-a11*a23 + a21*a13},
882  { a21*a32 - a31*a22,-a11*a32 + a31*a12, a11*a22 - a21*a12}};
883  //Matrix<T,3,3> rm;
884  //rm.get(0,0) = a22*a33 - a32*a23;
885  //rm.get(1,0) = -a21*a33 + a31*a23;
886  //rm.get(2,0) = a21*a32 - a31*a22;
887  //rm.get(0,1) = -a12*a33 + a32*a13;
888  //rm.get(1,1) = a11*a33 - a31*a13;
889  //rm.get(2,1) = -a11*a32 + a31*a12;
890  //rm.get(0,2) = a12*a23 - a22*a13;
891  //rm.get(1,2) = -a11*a23 + a21*a13;
892  //rm.get(2,2) = a11*a22 - a21*a12;
893  // Divide by the determinant
894  return rm / determinant(mat);
895 }
896 
901 template <class T>
902 constexpr Matrix<T,2,2> inverse(const Matrix<T,2,2> &mat) {
903  Matrix<T,2,2> rm;
904  rm.get(0,0) = mat.get(1,1);
905  rm.get(1,1) = mat.get(0,0);
906  rm.get(0,1) = mat.get(0,1)*(-1);
907  rm.get(1,0) = mat.get(1,0)*(-1);
908  //Divide by the determinant
909  return rm / determinant(mat);
910 }
911 
912 // Returns both inverse and determinant (Before inverse).
913 template <class T>
914 constexpr std::tuple<Matrix<T,3,3>,double> inverseDet(const Matrix<T,3,3> &mat) {
915  // First calculate the adjoint
916  const T &a11 = mat(0,0), a21 = mat(1,0), a31 = mat(2,0);
917  const T &a12 = mat(0,1), a22 = mat(1,1), a32 = mat(2,1);
918  const T &a13 = mat(0,2), a23 = mat(1,2), a33 = mat(2,2);
919 
920  Matrix<T,3,3> rm;
921  rm.get(0,0) = a22*a33 - a32*a23;
922  rm.get(1,0) = -a21*a33 + a31*a23;
923  rm.get(2,0) = a21*a32 - a31*a22;
924  rm.get(0,1) = -a12*a33 + a32*a13;
925  rm.get(1,1) = a11*a33 - a31*a13;
926  rm.get(2,1) = -a11*a32 + a31*a12;
927  rm.get(0,2) = a12*a23 - a22*a13;
928  rm.get(1,2) = -a11*a23 + a21*a13;
929  rm.get(2,2) = a11*a22 - a21*a12;
930  // Divide by the determinant
931  double det = determinant(mat);
932  return {rm / det,det};
933 }
934 
937 template <class T>
938 constexpr VMatrix<T,2> toMatrix(const Vector2<T> &v) {
939  VMatrix<T,2> ret;
940  ret.get(0) = v.x();
941  ret.get(1) = v.y();
942  return ret;
943 }
944 
947 template <class T>
948 constexpr VMatrix<T,3> toMatrix(const Vector3<T> &v) {
949  VMatrix<T,3> ret;
950  ret.get(0) = v.x();
951  ret.get(1) = v.y();
952  ret.get(2) = v.z();
953  return ret;
954 }
955 
957 template <class T,unsigned SX>
958 constexpr VMatrix<T,SX> toVMatrix(const Vector2<T> &v,unsigned start=0) {
959  VMatrix<T,SX> ret(0.0);
960  ret.get(start) = v.x();
961  ret.get(start+1) = v.y();
962  return ret;
963 }
964 
967 template <unsigned SX> constexpr DVMatrix<SX> toDVMatrix(const DVect2 &v,unsigned start=0) { return toVMatrix<double,SX>(v,start); }
968 
970 template <class T,unsigned SX>
971 constexpr VMatrix<T,SX> toVMatrix(const Vector3<T> &v,unsigned start=0) {
972  VMatrix<T,SX> ret(0.0);
973  ret.get(start) = v.x();
974  ret.get(start+1) = v.y();
975  ret.get(start+2) = v.z();
976  return ret;
977 }
978 
981 template <unsigned SX>
982 DVMatrix<SX> constexpr toDVMatrix(const DVect3 &v,unsigned start=0) { return toVMatrix<double,SX>(v,start); }
983 
986 template <class T>
987 constexpr Vector2<T> operator*(const Matrix<T,2,2> &m,const Vector2<T> &v) {
988  Vector2<T> ret(v.x()*m.get(0,0) + v.y()*m.get(0,1),
989  v.x()*m.get(1,0) + v.y()*m.get(1,1));
990  return ret;
991 }
992 
995 template <class T>
996 constexpr Vector2<T> operator*(const SymMatrix<T,2> &m,const Vector2<T> &v) {
997  Vector2<T> ret(v.x()*m.get(0,0) + v.y()*m.get(0,1),
998  v.x()*m.get(1,0) + v.y()*m.get(1,1));
999  return ret;
1000 }
1001 
1004 template <class T>
1005 constexpr Vector3<T> operator*(const Matrix<T,3,3> &m,const Vector3<T> &v) {
1006  Vector3<T> ret(v.x()*m.get(0,0) + v.y()*m.get(0,1) + v.z()*m.get(0,2),
1007  v.x()*m.get(1,0) + v.y()*m.get(1,1) + v.z()*m.get(1,2),
1008  v.x()*m.get(2,0) + v.y()*m.get(2,1) + v.z()*m.get(2,2));
1009  return ret;
1010 }
1011 
1014 template <class T>
1015 Vector3<T> operator*(const SymMatrix<T,3> &m,const Vector3<T> &v) {
1016  Vector3<T> ret(v.x()*m.get(0,0) + v.y()*m.get(0,1) + v.z()*m.get(0,2),
1017  v.x()*m.get(1,0) + v.y()*m.get(1,1) + v.z()*m.get(1,2),
1018  v.x()*m.get(2,0) + v.y()*m.get(2,1) + v.z()*m.get(2,2));
1019  return ret;
1020 }
1021 
1024 template <unsigned SX>
1025 DSymMatrix<SX> toSymMatrix(const SymTensor &s) {
1026  static_assert(SX <= 3);
1027  static_assert(SX >= 2);
1028  DSymMatrix<SX> ret;
1029  ret.get(0,0) = s.s11();
1030  ret.get(0,1) = ret.get(1,0) = s.s12();
1031  ret.get(1,1) = s.s22();
1032  if constexpr (SX > 2) {
1033  ret.get(0,2) = ret.get(2,0) = s.s13();
1034  ret.get(1,2) = ret.get(2,1) = s.s23();
1035  ret.get(2,2) = s.s33();
1036  }
1037  return ret;
1038 }
1039 
1042 BASE_EXPORT DMatrix<3,3> toMatrix(const SymTensor &s);
1043 
1046 BASE_EXPORT DMatrix<2,2> toMatrix2(const SymTensor &s);
1047 
1050 BASE_EXPORT DSymMatrix<3> toSymMatrix(const SymTensor &s);
1051 
1055 inline SymTensor toSymTensor(const DMatrix<3,3> &m) {
1056  SymTensor ret(m.get(0,0),
1057  m.get(1,1),
1058  m.get(2,2),
1059  (m.get(0,1)+m.get(1,0))*0.5,
1060  (m.get(0,2)+m.get(2,0))*0.5,
1061  (m.get(1,2)+m.get(2,1))*0.5);
1062  return ret;
1063 }
1064 
1065 inline SymTensor toSymTensor(const DMatrix<2,2> &m) {
1066  SymTensor ret(m.get(0,0),
1067  m.get(1,1),
1068  0.0,
1069  (m.get(0,1)+m.get(1,0))*0.5,
1070  0.0,
1071  0.0);
1072  return ret;
1073 }
1074 
1077 inline SymTensor toSymTensor(const DSymMatrix<3> &m) {
1078  SymTensor ret(m.get(0,0),m.get(1,1),m.get(1,2),
1079  m.get(0,1),m.get(0,2),m.get(1,2));
1080  return ret;
1081 }
1082 
1084 template <class T,unsigned SX>
1085 Vector2<T> toVector2(const VMatrix<T,SX> &m,unsigned start=0) { Vector2<T> ret(m(start),m(start+1)); return ret; }
1086 
1088 template <class T,unsigned SX>
1089 Vector3<T> toVector3(const VMatrix<T,SX> &m,unsigned start=0) { Vector3<T> ret(m(start),m(start+1),m(start+2)); return ret; }
1090 
1092 template <class T>
1093 Vector2<T> toVector(const VMatrix<T,2> &m) { Vector2<T> ret(m(0),m(1)); return ret; }
1094 
1096 template <class T>
1097 Vector3<T> toVector(const VMatrix<T,3> &m) { Vector3<T> ret(m(0),m(1),m(2)); return ret; }
1098 
1100 template <class T,unsigned SY>
1101 Vector2<T> columnToVector(const Matrix<T,2,SY> &m,unsigned col) { Vector2<T> ret(m(0,col),m(1,col)); return ret; }
1103 template <class T,unsigned SY>
1104 Vector3<T> columnToVector(const Matrix<T,3,SY> &m,unsigned col) { Vector3<T> ret(m(0,col),m(1,col),m(2,col)); return ret; }
1106 template <class T,unsigned SX>
1107 Vector2<T> rowToVector(const Matrix<T,SX,2> &m,unsigned row) { Vector2<T> ret(m(row,0),m(row,1)); return ret; }
1109 template <class T,unsigned SX>
1110 Vector3<T> rowToVector(const Matrix<T,SX,3> &m,unsigned row) { Vector3<T> ret(m(row,0),m(row,1),m(row,2)); return ret; }
1111 
1112 template <class T,unsigned SY>
1113 constexpr void vectorToColumn(Matrix<T,2,SY> &m,const DVect2 &v,unsigned col) { m(0,col) = v.x(); m(1,col) = v.y(); }
1114 
1115 template <class T,unsigned SY>
1116 constexpr void vectorToColumn(Matrix<T,3,SY> &m,const DVect3 &v,unsigned col) { m(0,col) = v.x(); m(1,col) = v.y(); m(2,col) = v.z(); }
1117 
1118 template <class T,unsigned SX,unsigned SY>
1119 constexpr void vectorToColumn(Matrix<T,SX,SY> &m,const VMatrix<T,SX> &v,unsigned col) { for (uint32 i=0;i<SX;++i) m(i,col) = v(i); }
1120 
1121 template <class T,unsigned SX>
1122 void vectorToRow(Matrix<T,SX,2> &m,const DVect2 &v,unsigned row) { m(row,0) = v.x(); m(row,1) = v.y(); }
1123 
1124 template <class T,unsigned SX>
1125 void vectorToRow(Matrix<T,SX,3> &m,const DVect3 &v,unsigned row) { m(row,0) = v.x(); m(row,1) = v.y(); m(row,2) = v.z(); }
1126 
1127 template <class T,unsigned SX,unsigned SY>
1128 constexpr void vectorToRow(Matrix<T,SX,SY> &m,const VMatrix<T,SY> &v,unsigned row) { for (uint32 j=0;j<SY;++j) m(row,j) = v(j); }
1129 
1130 template <class T,unsigned SX,unsigned SY>
1131 VMatrix<T,SX> column(const Matrix<T,SX,SY> &m,unsigned c) {
1132  VMatrix<T,SX> ret;
1133  for (unsigned i=0;i<SX;++i)
1134  ret(i) = m(i,c);
1135  return ret;
1136 }
1138 // EoF
DMatrix is a Matrix that defaults to type double...
Definition: matrix.h:741
constexpr DMatrix()
Default constructor, no data initialization.
Definition: matrix.h:746
constexpr DMatrix(const DMatrix< SX, SY > &m)
Copy constructor.
Definition: matrix.h:751
constexpr DMatrix(const Matrix< double, SX, SY > &m)
Copy constructor, for Matrix.
Definition: matrix.h:753
constexpr const DMatrix< SX, SY > & operator=(const DMatrix< SX, SY > &m)
Equality operator.
Definition: matrix.h:755
DSymMatrix is a SymMatrix that defaults to type double...
Definition: matrix.h:763
DSymMatrix(const double &t)
Explicit contructor, initializes all elements to t.
Definition: matrix.h:768
const DSymMatrix< SX > & operator=(const DSymMatrix< SX > &m)
Equality operator.
Definition: matrix.h:774
DSymMatrix(const SymMatrix< double, SX > &m)
Copy constructor, for SymMatrix.
Definition: matrix.h:772
DSymMatrix(const DSymMatrix< SX > &m)
Copy constructor.
Definition: matrix.h:770
DSymMatrix()
Default constructor, no data initialization.
Definition: matrix.h:766
const DSymMatrix< SX > & operator=(const SymMatrix< double, SX > &m)
Equality operator, for Matrix.
Definition: matrix.h:776
DVMatrix is a double version of VMatrix.
Definition: matrix.h:781
constexpr DVMatrix(const DVMatrix< S > &m)
Copy constructor.
Definition: matrix.h:791
constexpr DVMatrix(const Matrix< double, S, 1 > &m)
Copy constructor, for Matrix class.
Definition: matrix.h:793
A specialization of the Matrix class for the case when SY=1.
Definition: matrix.h:234
constexpr Matrix(const Matrix< T, SX, 1 > &m)
Copy constructor.
Definition: matrix.h:245
constexpr const Matrix< T, SX, 1 > & operator*=(const T &t)
In-place multiplication by a scalar.
Definition: matrix.h:273
Matrix< T, SX, 1 > operator*(const T &t) const
Binary scalar multiplication operator.
Definition: matrix.h:282
constexpr const T & get(unsigned x) const
Retrieve constant value at row x. Bounds checking is done in a debug compile.
Definition: matrix.h:254
constexpr T & operator()(unsigned x)
() operator access to get(x)
Definition: matrix.h:266
Matrix< T, SX, 1 > operator/(const T &t) const
Binary scalar division operator.
Definition: matrix.h:284
void addBlock(const Matrix< T, SX2, 1 > &m, unsigned iSrc, unsigned iDst)
Adds a block to this matrix from matrix m using the block size BX from block indice iSrc to indice iD...
Definition: matrix.h:300
void addGenBlock(const Matrix< T, SX2, SY2 > &m, unsigned iSrc, unsigned jSrc, unsigned iDst, unsigned jDst)
Adds a block to this matrix from matrix m using the block size (BX,BY), from normal not block indices...
Definition: matrix.h:304
void set(const T &t=T())
Sets all matrix elemnts to T.
Definition: matrix.h:320
constexpr const Matrix< T, SX, 1 > & operator/=(const T &t)
In-place division by a scalar.
Definition: matrix.h:275
constexpr const T & operator()(unsigned x, unsigned y) const
() operator access to const get(x,y)
Definition: matrix.h:260
constexpr const Matrix< T, SX, 1 > & operator=(const Matrix< T, SX, 1 > &m)
Equality operator.
Definition: matrix.h:249
Matrix< T, SX, 1 > operator+(const Matrix< T, SX, 1 > &m) const
Binary addition operator.
Definition: matrix.h:278
Matrix< T, SX, SZ > operator*(const Matrix< T, 1, SZ > &m) const
Binary matrix multiplication operator – simplest naive approach.
Definition: matrix.h:286
static Matrix< T, SX, 1 > identity()
Returns an identity matrix (or as close as you can get if not diagonal).
Definition: matrix.h:323
T maxNorm() const
Returns the infinity norm of the matrix, or the maximum absolute magnitude of any element.
Definition: matrix.h:316
constexpr const T & operator()(unsigned x) const
() operator access to const get(x)
Definition: matrix.h:262
constexpr const T & get(unsigned x,[[maybe_unused]] unsigned y) const
Retrieve constant value at row x column y. Bounds checking is done in a debug compile.
Definition: matrix.h:252
constexpr T & get(unsigned x,[[maybe_unused]] unsigned y)
Retrieve value at row x column y. Bounds checking is done in a debug compile.
Definition: matrix.h:256
constexpr const Matrix< T, SX, 1 > & operator+=(const Matrix< T, SX, 1 > &m)
In-place matrix addition.
Definition: matrix.h:269
constexpr const Matrix< T, SX, 1 > & operator-=(const Matrix< T, SX, 1 > &m)
In-place matrix subtraction.
Definition: matrix.h:271
Matrix< T, SX, 1 > operator-(const Matrix< T, SX, 1 > &m) const
Binary subtraction operator.
Definition: matrix.h:280
Matrix< T, 1, SX > transpose() const
returns the transposed matrix of this matrix
Definition: matrix.h:318
constexpr Matrix(const T &t)
Explicit constructor, initializes all elements to value t.
Definition: matrix.h:243
Matrix()
Default constructor, does nothing and no initialization.
Definition: matrix.h:240
void addGenBlock(const Matrix< T, SX2, 1 > &m, unsigned iSrc, unsigned iDst)
Adds a block to this matrix from matrix m using the block size BX from normal non block indice iSrc t...
Definition: matrix.h:310
constexpr T & operator()(unsigned x, unsigned y)
() operator access to get(x,y)
Definition: matrix.h:264
void addBlock(const Matrix< T, SX2, SY2 > &m, unsigned iSrc, unsigned jSrc, unsigned iDst, unsigned jDst)
Adds a block to this matrix from matrix m using the block size (BX,BY), from block indices (iSrc,...
Definition: matrix.h:296
constexpr T & get(unsigned x)
Retrieve value at row x. Bounds checking is done in a debug compile.
Definition: matrix.h:258
A template-based matrix class, size fixed at compile time. Defaults to symmetric sized matrix.
Definition: matrix.h:22
constexpr T & operator()(unsigned x, unsigned y)
() operator access to get(x,y)
Definition: matrix.h:139
void addGenBlock(const Matrix< T, SX2, SY2 > &m, unsigned iSrc, unsigned jSrc, unsigned iDst, unsigned jDst)
Adds a block to this matrix from matrix m using the block size (BX,BY), from normal not block indices...
Definition: matrix.h:197
static constexpr Matrix< T, SX, SY > identity()
Returns an identity matrix (or as close as you can get if not diagonal).
Definition: matrix.h:213
T maxNorm() const
Returns the infinity norm of the matrix, or the maximum absolute magnitude of any element.
Definition: matrix.h:204
Matrix< T, SX, SZ > operator*(const Matrix< T, SY, SZ > &m) const
Binary matrix multiplication operator – simplest naive approach.
Definition: matrix.h:160
constexpr const T & operator()(unsigned x, unsigned y) const
() operator access to get(x,y)
Definition: matrix.h:137
constexpr Matrix(const Matrix< T, SX, SY > &m)
Copy constructor.
Definition: matrix.h:124
constexpr POPWARNING Matrix(const T &t)
Explicit constructor, initializes all elements to value t.
Definition: matrix.h:120
constexpr const Matrix< T, SX, SY > & operator-=(const Matrix< T, SX, SY > &m)
In-place matrix subtraction.
Definition: matrix.h:144
constexpr const Matrix< T, SX, SY > & operator=(const Matrix< T, SX, SY > &m)
Equality operator.
Definition: matrix.h:126
constexpr const Matrix< T, SX, SY > & operator+=(const Matrix< T, SX, SY > &m)
In-place matrix addition.
Definition: matrix.h:142
constexpr const Matrix< T, SX, SY > & operator/=(const T &t)
In-place division by a scalar.
Definition: matrix.h:148
void addBlock(const Matrix< T, SX2, SY2 > &m, unsigned iSrc, unsigned jSrc, unsigned iDst, unsigned jDst)
Adds a block to this matrix from matrix m using the block size (BX,BY), from block indices (iSrc,...
Definition: matrix.h:192
T trace() const
Return the trace of the matrix or the sum of the diagonal components. Works also if the matrix is not...
Definition: matrix.h:208
constexpr const Column operator[](unsigned y) const
Array access operator returns a Column helper class, which has it's own [] operator to access a colum...
Definition: matrix.h:133
Matrix< T, SY, SX > transpose() const
Return the transpose of the matrix.
Definition: matrix.h:206
constexpr T & get(unsigned x, unsigned y)
Retrieve value at row x column y. Bounds checking is done in a debug compile.
Definition: matrix.h:131
constexpr const Matrix< T, SX, SY > & operator*=(const T &t)
In-place multiplication by a scalar.
Definition: matrix.h:146
constexpr Column operator[](unsigned y)
Array access operator returns a Column helper class, which has it's own [] operator to access a colum...
Definition: matrix.h:135
constexpr void set(const T &t=T())
Set all entries in the matrix to t.
Definition: matrix.h:210
Matrix< T, SX, SY > operator-(const Matrix< T, SX, SY > &m) const
Binary subtraction operator.
Definition: matrix.h:153
Matrix< T, SX, SY > operator/(const T &t) const
Binary scalar division operator.
Definition: matrix.h:157
PUSHWARNING VSWARNING(26495) Matrix()
Default constructor, does nothing and no initialization.
Definition: matrix.h:115
Matrix< T, SX, SY > operator*(const T &t) const
Binary scalar multiplication operator.
Definition: matrix.h:155
constexpr const T & get(unsigned x, unsigned y) const
Retrieve value at row x column y. Bounds checking is done in a debug compile.
Definition: matrix.h:129
Matrix< T, SX, SY > operator+(const Matrix< T, SX, SY > &m) const
Binary addition operator.
Definition: matrix.h:151
A template-based symmetric matrix class, size fixed at compile time. This primitive template-based ma...
Definition: matrix.h:332
const SymMatrix< T, SX > & operator/=(const T &t)
In-place division by a scalar.
Definition: matrix.h:507
SymMatrix(const T &t)
Explicit constructor, initializes all elements to value t.
Definition: matrix.h:481
static SymMatrix< T, SX > fromMatrix(const Matrix< T, SX, SX > &m)
Assign from full matrix.
Definition: matrix.h:604
SymMatrix< T, SX > operator*(const T &t) const
Binary scalar multiplication operator for a symetric matrix.
Definition: matrix.h:519
Matrix< T, SX, SZ > operator*(const Matrix< T, SX, SZ > &m) const
Binary matrix multiplication operator – simplest naive approach.
Definition: matrix.h:537
const SymMatrix< T, SX > & operator*=(const T &t)
In-place multiplication by a scalar.
Definition: matrix.h:505
SymMatrix(const SymMatrix< T, SX > &m)
Copy constructor.
Definition: matrix.h:483
T & operator()(unsigned x, unsigned y)
() operator access to get(x,y)
Definition: matrix.h:498
const Column operator[](unsigned y) const
Array access operator returns a Column helper class, which has it's own [] operator to access a colum...
Definition: matrix.h:492
void addBlock(const SymMatrix< T, SX2 > &m, unsigned iSrc, unsigned jSrc, unsigned iDst, unsigned jDst)
Adds a block to this matrix from matrix m using the block size (BX,BY), from block indices (iSrc,...
Definition: matrix.h:568
SymMatrix< T, SX > transpose() const
Return the transpose of the matrix.
Definition: matrix.h:593
void addGenBlock(const SymMatrix< T, SX2 > &m, unsigned iSrc, unsigned jSrc, unsigned iDst, unsigned jDst)
Adds a block to this matrix from matrix m using the block size (BX,BY), from normal not block indices...
Definition: matrix.h:581
const SymMatrix< T, SX > & operator+=(const SymMatrix< T, SX > &m)
In-place matrix addition.
Definition: matrix.h:501
void addBlock(const Matrix< T, SX2, SY2 > &m, unsigned iSrc, unsigned jSrc, unsigned iDst, unsigned jDst)
Adds a block to this matrix from matrix m using the block size (BX,BY), from block indices (iSrc,...
Definition: matrix.h:564
T & get(unsigned x, unsigned y)
Retrieve value at row x column y. Bounds checking is done in a debug compile.
Definition: matrix.h:490
SymMatrix< T, SX > operator-(const SymMatrix< T, SX > &m) const
Binary subtraction operator for a symetric matrix.
Definition: matrix.h:514
const SymMatrix< T, SX > & operator-=(const SymMatrix< T, SX > &m)
In-place matrix subtraction.
Definition: matrix.h:503
T trace() const
Return the trace of the matrix or the sum of the diagonal components. Works also if the matrix is not...
Definition: matrix.h:595
SymMatrix< T, SX > operator+(const SymMatrix< T, SX > &m) const
Binary addition operator for a symetric matris.
Definition: matrix.h:510
const SymMatrix< T, SX > & operator=(const SymMatrix< T, SX > &m)
Equality operator.
Definition: matrix.h:485
const T & operator()(unsigned x, unsigned y) const
() operator access to get(x,y)
Definition: matrix.h:496
void addGenBlock(const Matrix< T, SX2, SY2 > &m, unsigned iSrc, unsigned jSrc, unsigned iDst, unsigned jDst)
Adds a block to this matrix from matrix m using the block size (BX,BY), from normal not block indices...
Definition: matrix.h:572
T maxNorm() const
Returns the infinity norm of the matrix, or the maximum absolute magnitude of any element.
Definition: matrix.h:591
Matrix< T, SX, SX > operator+(const Matrix< T, SX, SX > &m) const
Binary addition operator.
Definition: matrix.h:512
void set(const T &t=T())
Sets all matrix elements to t.
Definition: matrix.h:597
static SymMatrix< T, SX > identity()
Returns an identity matrix (or as close as you can get if not diagonal).
Definition: matrix.h:602
Matrix< T, SX, SX > toMatrix() const
Returns its copy.
Definition: matrix.h:599
const T & get(unsigned x, unsigned y) const
Retrieve value at row x column y. Bounds checking is done in a debug compile.
Definition: matrix.h:488
Column operator[](unsigned y)
Array access operator returns a Column helper class, which has it's own [] operator to access a colum...
Definition: matrix.h:494
Matrix< T, SX, SX > operator*(const SymMatrix< T, SX > &m) const
Binary matrix multiplication operator – simplest naive approach.
Definition: matrix.h:523
SymMatrix< T, SX > operator/(const T &t) const
Binary scalar division operator for a symetric matrix.
Definition: matrix.h:521
SymMatrix()
Default constructor, does nothing and no initialization.
Definition: matrix.h:478
Matrix< T, SX, SX > operator-(const Matrix< T, SX, SX > &m) const
Binary subtraction operator.
Definition: matrix.h:516
A 1-Dimensional version of Matrix, to represent a vector.
Definition: matrix.h:707
constexpr VMatrix(const Matrix< T, S, 1 > &m)
Copy contructor, works on Matrix if SY is 1.
Definition: matrix.h:718
constexpr T & operator[](unsigned x)
1D version of array operator, which currently unfortunately eliminates [x][0] syntax on VMatrix (may ...
Definition: matrix.h:731
constexpr const T & operator[](unsigned x) const
1D version of array operator, which currently unfortunately eliminates [x][0] syntax on VMatrix (may ...
Definition: matrix.h:729
constexpr VMatrix(const VMatrix< T, S > &m)
Copy constructor.
Definition: matrix.h:716
Definition: matrix.h:641
debug checked shorthand for std::numeric_limits<T>::
Definition: limit.h:25
Definition: matrix.h:656
Definition: matrix.h:672
constexpr T determinant(const Matrix< T, 3, 3 > &mat)
Returns the determinant of a 3X3 Matrix.
Definition: matrix.h:854
constexpr VMatrix< T, SX > toVMatrix(const Vector2< T > &v, unsigned start=0)
Converts a Vector2 into a VMatrix of arbitrary size, at an arbitrary starting index.
Definition: matrix.h:958
BASE_EXPORT DMatrix< 2, 2 > toMatrix2(const SymTensor &s)
Definition: matrix.cpp:16
Vector2< T > toVector(const VMatrix< T, 2 > &m)
Converts a VMatrix to a Vector3, using three elements starting at index start.
Definition: matrix.h:1093
constexpr Matrix< T, 3, 3 > inverse(const Matrix< T, 3, 3 > &mat)
Returns the inverse of a 3X3 Matrix.
Definition: matrix.h:874
constexpr Vector3< T > operator*(const Matrix< T, 3, 3 > &m, const Vector3< T > &v)
Definition: matrix.h:1005
Matrix< T, 3, 3 > outerProduct(const Vector3< T > &v1, const Vector3< T > &v2)
Creates a 3X3 Matrix from the outer product of two Vector3 types.
Definition: matrix.h:822
constexpr VMatrix< T, 2 > toMatrix(const Vector2< T > &v)
Definition: matrix.h:938
DSymMatrix< SX > toSymMatrix(const SymTensor &s)
Definition: matrix.h:1025
constexpr DVMatrix< SX > toDVMatrix(const DVect2 &v, unsigned start=0)
Definition: matrix.h:967
constexpr VMatrix< T, 3 > toMatrix(const Vector3< T > &v)
Definition: matrix.h:948
Vector2< T > toVector2(const VMatrix< T, SX > &m, unsigned start=0)
Converts a VMatrix to a Vector2, using two elements starting at index start.
Definition: matrix.h:1085
SymTensor toSymTensor(const DSymMatrix< 3 > &m)
Definition: matrix.h:1077
SymTensor toSymTensor(const DMatrix< 3, 3 > &m)
Definition: matrix.h:1055
constexpr Vector2< T > operator*(const Matrix< T, 2, 2 > &m, const Vector2< T > &v)
Definition: matrix.h:987
constexpr Vector2< T > operator*(const SymMatrix< T, 2 > &m, const Vector2< T > &v)
Definition: matrix.h:996
Vector2< T > rowToVector(const Matrix< T, SX, 2 > &m, unsigned row)
returns in a Vector2<t> the row row from matrix Matrix<T,SX,2>
Definition: matrix.h:1107
constexpr T determinant(const Matrix< T, 2, 2 > &mat)
Returns the determinant of a 2X2 Matrix.
Definition: matrix.h:865
#define BASE_EXPORT
Definition: basedef.h:24
Vector3< T > toVector3(const VMatrix< T, SX > &m, unsigned start=0)
Converts a VMatrix to a Vector3, using three elements starting at index start.
Definition: matrix.h:1089
Vector3< T > operator*(const SymMatrix< T, 3 > &m, const Vector3< T > &v)
Definition: matrix.h:1015
Vector2< T > columnToVector(const Matrix< T, 2, SY > &m, unsigned col)
returns in a Vector2<t> the column col from matrix Matrix<T,SX,2>
Definition: matrix.h:1101
Matrix< T, 2, 2 > outerProduct(const Vector2< T > &v1, const Vector2< T > &v2)
Creates a 2X2 Matrix from the outer product of two Vector2 types.
Definition: matrix.h:840
constexpr Matrix< T, 2, 2 > inverse(const Matrix< T, 2, 2 > &mat)
Returns the inverse of a 2X2 Matrix.
Definition: matrix.h:902
A Symmetric 2nd order tensor.