MorphoGraphX
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Matrix.hpp
Go to the documentation of this file.
1 #ifndef MATRIX_HPP
2 #define MATRIX_HPP
3 
10 #include <Config.hpp>
11 
12 #include <cuda/CudaGlobal.hpp>
13 #include <StaticAssert.hpp>
14 #include <Util.hpp>
15 #include <Vector.hpp>
16 
17 #include <cmath>
18 #ifndef COMPILE_CUDA
19 # include <QTextStream>
20 #endif
21 
22 namespace mgx {
23 namespace util {
24 
25 enum MatrixLayout { C_STYLE, GL_STYLE };
26 
34 template <size_t nRows, size_t nCols, typename T = double> class Matrix {
35 private:
36  Vector<nCols, T> rows[nRows];
37 
38 public:
39  typedef T value_type;
40  typedef T& reference_type;
41  typedef const T& const_reference_type;
42  typedef T* pointer_type;
43  typedef const T* const_pointer_type;
44  typedef T* iterator;
45  typedef const T* const_iterator;
46 
47  static const size_t numcols = nCols;
48  static const size_t numrows = nRows;
49 
53  CU_HOST_DEVICE
54  Matrix(void)
55  {
56  for(size_t i = 0; i < nRows; i++)
57  for(size_t j = 0; j < nCols; j++)
58  rows[i][j] = 0;
59  }
60 
66  template <typename T1> CU_HOST_DEVICE Matrix(const Matrix<nRows, nCols, T1>& mat)
67  {
68  for(size_t i = 0; i < nRows; i++)
69  for(size_t j = 0; j < nCols; j++)
70  rows[i][j] = T(mat(i, j));
71  }
72 
78  template <typename T1> CU_HOST_DEVICE Matrix(const Vector<nCols, T1>* vecs)
79  {
80  for(size_t i = 0; i < nRows; i++)
81  rows[i] = vecs[i];
82  }
83 
93  template <typename T1> CU_HOST_DEVICE Matrix(const T1* values, MatrixLayout layout = C_STYLE)
94  {
95  switch(layout) {
96  case C_STYLE: {
97  for(size_t i = 0; i < nRows; i++) {
98  rows[i] = Vector<nCols, T>(values + (i * nCols));
99  }
100  } break;
101  case GL_STYLE: {
102  for(size_t j = 0; j < nCols; j++)
103  for(size_t i = 0; i < nRows; i++) {
104  rows[i][j] = T(values[i + j * nRows]);
105  }
106  } break;
107  }
108  }
109 
119  CU_HOST_DEVICE
120  Matrix(const T* values, MatrixLayout layout = C_STYLE)
121  {
122  switch(layout) {
123  case C_STYLE: {
124  for(size_t i = 0; i < nRows; i++) {
125  rows[i] = Vector<nCols, T>(values + (i * nCols));
126  }
127  } break;
128  case GL_STYLE: {
129  for(size_t i = 0; i < nRows; i++)
130  for(size_t j = 0; j < nCols; j++) {
131  rows[i][j] = values[i + j * nRows];
132  }
133  } break;
134  }
135  }
136 
142  CU_HOST_DEVICE
143  Matrix(const T& value)
144  {
145  for(size_t i = 0; i < nRows; i++)
146  for(size_t j = 0; j < nCols; j++)
147  rows[i][j] = (i == j) ? value : 0;
148  }
149 
155  CU_HOST_DEVICE
157  return Vector<2, size_t>(nRows, nCols);
158  }
159 
163  CU_HOST_DEVICE
164  static size_t nbRows() {
165  return nRows;
166  }
167 
171  CU_HOST_DEVICE
172  static size_t nbColumns() {
173  return nCols;
174  }
175 
181  CU_HOST_DEVICE
182  const T* c_data() const {
183  return rows[0].c_data();
184  }
185 
191  CU_HOST_DEVICE
192  T* data() {
193  return rows[0].data();
194  }
195 
199  CU_HOST_DEVICE
200  Matrix operator-(void) const
201  {
202  Matrix ans;
203 
204  for(size_t i = 0; i < nRows; i++)
205  ans.rows[i] = -rows[i];
206 
207  return ans;
208  }
209 
213  CU_HOST_DEVICE
214  Matrix operator+(const Matrix& mat) const
215  {
216  Matrix ans;
217 
218  for(size_t i = 0; i < nRows; i++)
219  ans.rows[i] = rows[i] + mat.rows[i];
220 
224  return ans;
225  }
226 
227  CU_HOST_DEVICE
228  Matrix operator-(const Matrix& mat) const
229  {
230  Matrix ans;
231 
232  for(size_t i = 0; i < nRows; i++)
233  ans.rows[i] = rows[i] - mat.rows[i];
234 
235  return ans;
236  }
237 
241  CU_HOST_DEVICE
242  Matrix operator*(const T& scalar) const
243  {
244  Matrix ans;
245 
246  for(size_t i = 0; i < nRows; i++)
247  ans[i] = rows[i] * scalar;
248 
249  return ans;
250  }
251 
255  CU_HOST_DEVICE
256  Matrix operator/(const T& scalar) const
257  {
258  Matrix ans;
259 
260  for(size_t i = 0; i < nRows; i++)
261  ans[i] = rows[i] / scalar;
262 
263  return ans;
264  }
265 
269  CU_HOST_DEVICE
270  friend Matrix operator*(const T& scalar, const Matrix& mat)
271  {
272  Matrix ans;
273 
274  for(size_t i = 0; i < nRows; i++)
275  ans[i] = scalar * mat.rows[i];
276 
277  return ans;
278  }
279 
283  CU_HOST_DEVICE
285  {
286  Vector<nRows, T> ans;
287  for(size_t i = 0; i < nRows; ++i) {
288  T value = 0;
289  for(size_t j = 0; j < nCols; ++j)
290  value += rows[i][j] * vec[j];
291  ans[i] = value;
292  }
293  return ans;
294  }
295 
296  CU_HOST_DEVICE
297  Matrix& operator=(const Matrix& mat)
298  {
299  for(size_t i = 0; i < nRows; i++)
300  rows[i] = mat.rows[i];
301 
302  return (*this);
303  }
304 
305  CU_HOST_DEVICE
306  Matrix& operator+=(const Matrix& mat) {
307  return ((*this) = (*this) + mat);
308  }
309 
310  CU_HOST_DEVICE
311  Matrix& operator-=(const Matrix& mat) {
312  return ((*this) = (*this) - mat);
313  }
314 
315  CU_HOST_DEVICE
316  Matrix& operator*=(const T& scalar) {
317  return ((*this) = (*this) * scalar);
318  }
319 
320  CU_HOST_DEVICE
321  Matrix& operator/=(const T& scalar) {
322  return ((*this) = (*this) / scalar);
323  }
324 
325  CU_HOST_DEVICE
326  Matrix& operator*=(const Matrix& mat) {
327  return ((*this) = (*this) * mat);
328  }
329 
330  bool operator==(const Matrix& mat) const
331  {
332  for(size_t i = 0; i < nRows; i++)
333  if(rows[i] != mat.rows[i])
334  return false;
335 
336  return true;
337  }
338 
339  CU_HOST_DEVICE
340  bool operator!=(const Matrix& mat) const {
341  return (!((*this) == mat));
342  }
343 
344 #ifndef COMPILE_CUDA
345  friend QTextStream& operator<<(QTextStream& out, const Matrix& mat)
346  {
347  for(size_t i = 0; i < nRows; i++) {
348  out << mat.rows[i];
349  if(i != (nRows - 1))
350  out << " ";
351  }
352 
353  return out;
354  }
355 
356  friend QTextStream& operator>>(QTextStream& in, Matrix& mat)
357  {
358  for(size_t i = 0; i < nRows && !in.atEnd(); i++)
359  in >> mat.rows[i];
360  return in;
361  }
362 #endif
363 
364  CU_HOST_DEVICE
365  friend std::ostream& operator<<(std::ostream& out, const Matrix& mat)
366  {
367  for(size_t i = 0; i < nRows; i++) {
368  out << mat.rows[i];
369  if(i != (nRows - 1))
370  out << " ";
371  }
372 
373  return out;
374  }
375 
376  CU_HOST_DEVICE
377  friend std::istream& operator>>(std::istream& in, Matrix& mat)
378  {
379  for(size_t i = 0; i < nRows && in; i++)
380  in >> mat.rows[i];
381  return in;
382  }
383 
389  CU_HOST_DEVICE
391  return rows[idx];
392  }
393 
399  CU_HOST_DEVICE
400  Vector<nCols, T> operator[](size_t idx) const {
401  return rows[idx];
402  }
403 
407  CU_HOST_DEVICE
408  T& operator()(size_t i, size_t j) {
409  return rows[i][j];
410  }
411 
415  CU_HOST_DEVICE
416  T operator()(size_t i, size_t j) const {
417  return rows[i][j];
418  }
419 
423  CU_HOST_DEVICE
424  static Matrix identity()
425  {
426  Matrix mat(1);
427  return mat;
428  }
429 
433  CU_HOST_DEVICE
434  Matrix& zero(void)
435  {
436  for(size_t i = 0; i < nRows; i++)
437  for(size_t j = 0; j < nCols; j++)
438  rows[i][j] = 0.0;
439  return (*this);
440  }
441 
447  CU_HOST_DEVICE
448  Matrix& operator=(const T& value)
449  {
450  for(size_t i = 0; i < nRows; ++i) {
451  for(size_t j = 0; j < nCols; ++j) {
452  if(i == j)
453  rows[i][j] = value;
454  else
455  rows[i][j] = 0;
456  }
457  }
458  return *this;
459  }
460 
464  CU_HOST_DEVICE
466  {
468  for(size_t i = 0; i < nRows; ++i)
469  for(size_t j = 0; j < nCols; ++j)
470  t[i][j] = rows[j][i];
471  return t;
472  }
473 
480  CU_HOST_DEVICE
481  static Matrix<3, 3, T> rotation(const Vector<3, T>& direction, T angle)
482  {
483  T ca = std::cos(angle);
484  T sa = std::sin(angle);
485  Matrix<3, 3, T> r;
486  T x = direction.x();
487  T y = direction.y();
488  T z = direction.z();
489  r[0].set(ca + (1 - ca) * x * x, (1 - ca) * x * y - sa * z, (1 - ca) * z * x + sa * y);
490  r[1].set((1 - ca) * y * x + sa * z, ca + (1 - ca) * y * y, (1 - ca) * z * y - sa * x);
491  r[2].set((1 - ca) * x * z - sa * y, (1 - ca) * y * z + sa * x, ca + (1 - ca) * z * z);
492  return r;
493  }
494 
501  CU_HOST_DEVICE
502  static Matrix<4, 4, T> rotation(const Vector<4, T>& direction, T angle)
503  {
504  T ca = std::cos(angle);
505  T sa = std::sin(angle);
506  Matrix<4, 4, T> r;
507  T x = direction.x() / direction.t();
508  T y = direction.y() / direction.t();
509  T z = direction.z() / direction.t();
510  r[0].set(ca + (1 - ca) * x * x, (1 - ca) * x * y - sa * z, (1 - ca) * z * x + sa * y, 0);
511  r[1].set((1 - ca) * y * x + sa * z, ca + (1 - ca) * y * y, (1 - ca) * z * y - sa * x, 0);
512  r[2].set((1 - ca) * x * z - sa * y, (1 - ca) * y * z + sa * x, ca + (1 - ca) * z * z, 0);
513  r[3].set(0, 0, 0, 1);
514  return r;
515  }
516 
520  CU_HOST_DEVICE
521  T trace() const
522  {
523  T acc = 0;
524  for(size_t i = 0; i < min(nRows, nCols); ++i) {
525  acc += rows[i][i];
526  }
527  return acc;
528  }
529 
530  CU_HOST_DEVICE
531  void fillArray(T* array, MatrixLayout layout = C_STYLE)
532  {
533  switch(layout) {
534  case C_STYLE:
535  memcpy(array, &rows[0][0], sizeof(T) * nRows * nCols);
536  break;
537  case GL_STYLE: {
538  size_t k = 0;
539  for(size_t c = 0; c < nCols; ++c)
540  for(size_t r = 0; r < nRows; ++r, ++k)
541  array[k] = rows[r][c];
542  } break;
543  }
544  }
545 
549  CU_HOST_DEVICE
551  {
552  STATIC_ASSERT(nRows == nCols);
553  Vector<nRows, T> res;
554  for(size_t i = 0; i < nRows; ++i)
555  res[i] = (*this)(i, i);
556  return res;
557  }
558 };
559 
564 template <size_t nRows, size_t nCols, typename T>
565 CU_HOST_DEVICE Vector<nCols, T> operator*(const Vector<nCols, T>& vec, const Matrix<nRows, nCols, T>& mat);
566 
571 template <size_t nRows, size_t nSize, size_t nCols, typename T>
572 CU_HOST_DEVICE Matrix<nRows, nCols, T> operator*(const Matrix<nRows, nSize, T>& mat1,
573  const Matrix<nSize, nCols, T>& mat2);
574 
579 template <typename T> T CU_HOST_DEVICE det(const Matrix<1, 1, T>& mat);
584 template <typename T> T CU_HOST_DEVICE det(const Matrix<2, 2, T>& mat);
589 template <typename T> T CU_HOST_DEVICE det(const Matrix<3, 3, T>& mat);
596 template <size_t nRows, typename T> T CU_HOST_DEVICE det(const Matrix<nRows, nRows, T>& mat);
597 
598 template <size_t nRows, typename T>
599 CU_HOST_DEVICE T matrix_minor(const Matrix<nRows, nRows, T>& mat, size_t i, size_t j);
600 
605 template <size_t nRows, typename T> CU_HOST_DEVICE T cofactor(const Matrix<nRows, nRows, T>& mat, size_t i, size_t j);
606 
612 template <typename T> CU_HOST_DEVICE Matrix<1, 1, T> inverse(const Matrix<1, 1, T>& mat);
618 template <typename T> CU_HOST_DEVICE Matrix<2, 2, T> inverse(const Matrix<2, 2, T>& mat);
624 template <typename T> CU_HOST_DEVICE Matrix<3, 3, T> inverse(const Matrix<3, 3, T>& mat);
630 template <typename T> CU_HOST_DEVICE Matrix<4, 4, T> inverse(const Matrix<4, 4, T>& mat);
631 
639 template <size_t nRows, typename T> CU_HOST_DEVICE Matrix<nRows, nRows, T> inverse(const Matrix<nRows, nRows, T>& mat);
640 
646 template <size_t nRows, size_t nCols, typename T>
647 CU_HOST_DEVICE Matrix<nCols, nRows, T> transpose(const Matrix<nRows, nCols, T>& mat);
648 
657 template <size_t nRows, size_t nCols, typename T> CU_HOST_DEVICE T norm(const Matrix<nRows, nCols, T>& mat);
658 
666 template <size_t nRows, size_t nCols, typename T> CU_HOST_DEVICE T normsq(const Matrix<nRows, nCols, T>& mat);
667 
673 template <size_t nRows, size_t nCols, typename T>
674 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(const T& (*fct)(const T &), const Matrix<nRows, nCols, T>& m)
675 {
677  for(size_t i = 0; i < nRows; ++i) {
678  const Vector<nCols, T>& mrow = m[i];
679  Vector<nCols, T>& rrow = result[i];
680  for(size_t j = 0; j < nCols; ++j) {
681  rrow[j] = (*fct)(mrow[j]);
682  }
683  }
684  return result;
685 }
686 
692 template <size_t nRows, size_t nCols, typename T>
693 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(T (*fct)(T), const Matrix<nRows, nCols, T>& m)
694 {
696  for(size_t i = 0; i < nRows; ++i) {
697  const Vector<nCols, T>& mrow = m[i];
698  Vector<nCols, T>& rrow = result[i];
699  for(size_t j = 0; j < nCols; ++j) {
700  rrow[j] = (*fct)(mrow[j]);
701  }
702  }
703  return result;
704 }
705 
711 template <size_t nRows, size_t nCols, typename T>
712 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(T (*fct)(const T&), const Matrix<nRows, nCols, T>& m)
713 {
715  for(size_t i = 0; i < nRows; ++i) {
716  const Vector<nCols, T>& mrow = m[i];
717  Vector<nCols, T>& rrow = result[i];
718  for(size_t j = 0; j < nCols; ++j) {
719  rrow[j] = (*fct)(mrow[j]);
720  }
721  }
722  return result;
723 }
724 
730 template <size_t nRows, size_t nCols, typename T, typename T1>
731 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(const T& (*fct)(const T1 &), const Matrix<nRows, nCols, T1>& m)
732 {
734  for(size_t i = 0; i < nRows; ++i) {
735  const Vector<nCols, T1>& mrow = m[i];
736  Vector<nCols, T>& rrow = result[i];
737  for(size_t j = 0; j < nCols; ++j) {
738  rrow[j] = (*fct)(mrow[j]);
739  }
740  }
741  return result;
742 }
743 
749 template <size_t nRows, size_t nCols, typename T, typename T1>
750 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(T (*fct)(T1), const Matrix<nRows, nCols, T1>& m)
751 {
753  for(size_t i = 0; i < nRows; ++i) {
754  const Vector<nCols, T1>& mrow = m[i];
755  Vector<nCols, T>& rrow = result[i];
756  for(size_t j = 0; j < nCols; ++j) {
757  rrow[j] = (*fct)(mrow[j]);
758  }
759  }
760  return result;
761 }
762 
768 template <size_t nRows, size_t nCols, typename T, typename T1>
769 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(T (*fct)(const T1&), const Matrix<nRows, nCols, T1>& m)
770 {
772  for(size_t i = 0; i < nRows; ++i) {
773  const Vector<nCols, T1>& mrow = m[i];
774  Vector<nCols, T>& rrow = result[i];
775  for(size_t j = 0; j < nCols; ++j) {
776  rrow[j] = (*fct)(mrow[j]);
777  }
778  }
779  return result;
780 }
781 
787 template <size_t nRows, size_t nCols, typename T>
788 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(T (*fct)(T, T), const Matrix<nRows, nCols, T>& m1,
789  const Matrix<nRows, nCols, T>& m2)
790 {
792  for(size_t i = 0; i < nRows; ++i) {
793  const Vector<nCols, T>& mrow1 = m1[i];
794  const Vector<nCols, T>& mrow2 = m2[i];
795  Vector<nCols, T>& rrow = result[i];
796  for(size_t j = 0; j < nCols; ++j) {
797  rrow[j] = (*fct)(mrow1[j], mrow2[j]);
798  }
799  }
800  return result;
801 }
802 
808 template <size_t nRows, size_t nCols, typename T>
809 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(T (*fct)(const T&, const T&), const Matrix<nRows, nCols, T>& m1,
810  const Matrix<nRows, nCols, T>& m2)
811 {
813  for(size_t i = 0; i < nRows; ++i) {
814  const Vector<nCols, T>& mrow1 = m1[i];
815  const Vector<nCols, T>& mrow2 = m2[i];
816  Vector<nCols, T>& rrow = result[i];
817  for(size_t j = 0; j < nCols; ++j) {
818  rrow[j] = (*fct)(mrow1[j], mrow2[j]);
819  }
820  }
821  return result;
822 }
823 
829 template <size_t nRows, size_t nCols, typename T>
830 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(const T& (*fct)(const T &, const T &), const Matrix<nRows, nCols, T>& m1,
831  const Matrix<nRows, nCols, T>& m2)
832 {
834  for(size_t i = 0; i < nRows; ++i) {
835  const Vector<nCols, T>& mrow1 = m1[i];
836  const Vector<nCols, T>& mrow2 = m2[i];
837  Vector<nCols, T>& rrow = result[i];
838  for(size_t j = 0; j < nCols; ++j) {
839  rrow[j] = (*fct)(mrow1[j], mrow2[j]);
840  }
841  }
842  return result;
843 }
844 
850 template <size_t nRows, size_t nCols, typename T, typename T1, typename T2>
851 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(T (*fct)(T1, T2), const Matrix<nRows, nCols, T1>& m1,
852  const Matrix<nRows, nCols, T2>& m2)
853 {
855  for(size_t i = 0; i < nRows; ++i) {
856  const Vector<nCols, T1>& mrow1 = m1[i];
857  const Vector<nCols, T2>& mrow2 = m2[i];
858  Vector<nCols, T>& rrow = result[i];
859  for(size_t j = 0; j < nCols; ++j) {
860  rrow[j] = (*fct)(mrow1[j], mrow2[j]);
861  }
862  }
863  return result;
864 }
865 
871 template <size_t nRows, size_t nCols, typename T, typename T1, typename T2>
872 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(T (*fct)(const T1&, const T2&), const Matrix<nRows, nCols, T1>& m1,
873  const Matrix<nRows, nCols, T2>& m2)
874 {
876  for(size_t i = 0; i < nRows; ++i) {
877  const Vector<nCols, T1>& mrow1 = m1[i];
878  const Vector<nCols, T2>& mrow2 = m2[i];
879  Vector<nCols, T>& rrow = result[i];
880  for(size_t j = 0; j < nCols; ++j) {
881  rrow[j] = (*fct)(mrow1[j], mrow2[j]);
882  }
883  }
884  return result;
885 }
886 
892 template <size_t nRows, size_t nCols, typename T, typename T1, typename T2>
893 CU_HOST_DEVICE Matrix<nRows, nCols, T> map(const T& (*fct)(const T1 &, const T2 &), const Matrix<nRows, nCols, T1>& m1,
894  const Matrix<nRows, nCols, T2>& m2)
895 {
897  for(size_t i = 0; i < nRows; ++i) {
898  const Vector<nCols, T1>& mrow1 = m1[i];
899  const Vector<nCols, T2>& mrow2 = m2[i];
900  Vector<nCols, T>& rrow = result[i];
901  for(size_t j = 0; j < nCols; ++j) {
902  rrow[j] = (*fct)(mrow1[j], mrow2[j]);
903  }
904  }
905  return result;
906 }
907 
908 template <size_t nRows, size_t nCols, typename T>
909 CU_HOST_DEVICE Vector<nCols, T> operator*(const Vector<nCols, T>& vec, const Matrix<nRows, nCols, T>& mat)
910 {
911  Vector<nCols, T> ans;
912  for(size_t i = 0; i < nCols; ++i) {
913  T value = 0;
914  for(size_t j = 0; j < nRows; ++j)
915  value += mat(j, i) * vec[j];
916  ans[i] = value;
917  }
918  return ans;
919 }
920 
921 template <size_t nRows, size_t intSize, size_t nCols, typename T>
922 CU_HOST_DEVICE Matrix<nRows, nCols, T> operator*(const Matrix<nRows, intSize, T>& mat1,
924 {
926  for(size_t i = 0; i < nRows; i++)
927  for(size_t j = 0; j < nCols; j++) {
928  T acc = 0;
929  for(size_t k = 0; k < intSize; k++)
930  acc += mat1(i, k) * mat2(k, j);
931  ans(i, j) = acc;
932  }
933  return ans;
934 }
935 
936 template <typename T> CU_HOST_DEVICE T det(const Matrix<1, 1, T>& mat) {
937  return mat(0, 0);
938 }
939 
940 template <typename T> CU_HOST_DEVICE T det(const Matrix<2, 2, T>& mat)
941 {
942  return mat(0, 0) * mat(1, 1) - mat(0, 1) * mat(1, 0);
943 }
944 
945 template <typename T> CU_HOST_DEVICE T det(const Matrix<3, 3, T>& mat)
946 {
947  return mat(0, 0) * mat(1, 1) * mat(2, 2) + mat(0, 1) * mat(1, 2) * mat(2, 0) + mat(0, 2) * mat(1, 0) * mat(2, 1)
948  - mat(0, 0) * mat(1, 2) * mat(2, 1) - mat(0, 1) * mat(1, 0) * mat(2, 2) - mat(0, 2) * mat(1, 1) * mat(2, 0);
949 }
950 
951 template <size_t nRows, typename T> CU_HOST_DEVICE T det(const Matrix<nRows, nRows, T>& mat)
952 {
953  STATIC_ASSERT(nRows > 3);
954  T acc = 0;
955  for(size_t i = 0; i < nRows; i++) {
956  acc += mat(i, 0) * cofactor(mat, i, 0);
957  }
958  return acc;
959 }
960 
961 template <size_t nRows, typename T>
962 CU_HOST_DEVICE T matrix_minor(const Matrix<nRows, nRows, T>& mat, size_t i, size_t j)
963 {
964  STATIC_ASSERT(nRows > 0);
965  Matrix<nRows - 1, nRows - 1, T> ans;
966  for(size_t i1 = 0, i2 = 0; i1 < nRows; i1++, i2++) {
967  if(i1 == i) {
968  i2--;
969  } else {
970  for(size_t j1 = 0, j2 = 0; j1 < nRows; j1++, j2++) {
971  if(j1 == j) {
972  j2--;
973  } else {
974  ans(i2, j2) = mat(i1, j1);
975  }
976  }
977  }
978  }
979  return det(ans);
980 }
981 
982 template <size_t nRows, typename T> CU_HOST_DEVICE T cofactor(const Matrix<nRows, nRows, T>& mat, size_t i, size_t j)
983 {
984  T inv = 1;
985  if((i + j) % 2) {
986  inv = -1;
987  }
988  return inv * matrix_minor(mat, i, j);
989 }
990 
991 template <typename T> CU_HOST_DEVICE Matrix<1, 1, T> inverse(const Matrix<1, 1, T>& mat)
992 {
993  Matrix<1, 1, T> ans;
994  ans[0][0] = 1 / mat(0, 0);
995  return ans;
996 }
997 
998 template <typename T> CU_HOST_DEVICE Matrix<2, 2, T> inverse(const Matrix<2, 2, T>& mat)
999 {
1000  Matrix<2, 2, T> ans;
1001  T d;
1002  d = det(mat);
1003  if(d == 0)
1004  return ans;
1005  ans(0, 0) = mat(1, 1) / d;
1006  ans(0, 1) = mat(0, 1) / -d;
1007  ans(1, 0) = mat(1, 0) / -d;
1008  ans(1, 1) = mat(0, 0) / d;
1009  return ans;
1010 }
1011 
1012 template <typename T> CU_HOST_DEVICE Matrix<3, 3, T> inverse(const Matrix<3, 3, T>& mat)
1013 {
1014  Matrix<3, 3, T> ans;
1015  T d;
1016  d = det(mat);
1017  if(d == 0)
1018  return ans;
1019  ans(0, 0) = (mat(1, 1) * mat(2, 2) - mat(1, 2) * mat(2, 1)) / d;
1020  ans(1, 1) = (mat(0, 0) * mat(2, 2) - mat(0, 2) * mat(2, 0)) / d;
1021  ans(2, 2) = (mat(1, 1) * mat(0, 0) - mat(1, 0) * mat(0, 1)) / d;
1022  ans(1, 0) = (mat(1, 2) * mat(2, 0) - mat(1, 0) * mat(2, 2)) / d;
1023  ans(0, 1) = (mat(2, 1) * mat(0, 2) - mat(0, 1) * mat(2, 2)) / d;
1024  ans(2, 0) = (mat(1, 0) * mat(2, 1) - mat(1, 1) * mat(2, 0)) / d;
1025  ans(0, 2) = (mat(0, 1) * mat(1, 2) - mat(1, 1) * mat(0, 2)) / d;
1026  ans(1, 2) = (mat(0, 2) * mat(1, 0) - mat(0, 0) * mat(1, 2)) / d;
1027  ans(2, 1) = (mat(2, 0) * mat(0, 1) - mat(0, 0) * mat(2, 1)) / d;
1028  return ans;
1029 }
1030 
1031 template <typename T> Matrix<4, 4, T> inverse(const Matrix<4, 4, T>& m)
1032 {
1033  Matrix<4, 4, T> inv;
1034 
1035  inv(0, 0) = (m(1, 1) * m(2, 2) * m(3, 3) - m(1, 1) * m(3, 2) * m(2, 3) - m(1, 2) * m(2, 1) * m(3, 3)
1036  + m(1, 2) * m(3, 1) * m(2, 3) + m(1, 3) * m(2, 1) * m(3, 2) - m(1, 3) * m(3, 1) * m(2, 2));
1037 
1038  inv(0, 1) = (-m(0, 1) * m(2, 2) * m(3, 3) + m(0, 1) * m(3, 2) * m(2, 3) + m(0, 2) * m(2, 1) * m(3, 3)
1039  - m(0, 2) * m(3, 1) * m(2, 3) - m(0, 3) * m(2, 1) * m(3, 2) + m(0, 3) * m(3, 1) * m(2, 2));
1040 
1041  inv(0, 2) = (m(0, 1) * m(1, 2) * m(3, 3) - m(0, 1) * m(3, 2) * m(1, 3) - m(0, 2) * m(1, 1) * m(3, 3)
1042  + m(0, 2) * m(3, 1) * m(1, 3) + m(0, 3) * m(1, 1) * m(3, 2) - m(0, 3) * m(3, 1) * m(1, 2));
1043 
1044  inv(0, 3) = (-m(0, 1) * m(1, 2) * m(2, 3) + m(0, 1) * m(2, 2) * m(1, 3) + m(0, 2) * m(1, 1) * m(2, 3)
1045  - m(0, 2) * m(2, 1) * m(1, 3) - m(0, 3) * m(1, 1) * m(2, 2) + m(0, 3) * m(2, 1) * m(1, 2));
1046 
1047  inv(1, 0) = (-m(1, 0) * m(2, 2) * m(3, 3) + m(1, 0) * m(3, 2) * m(2, 3) + m(1, 2) * m(2, 0) * m(3, 3)
1048  - m(1, 2) * m(3, 0) * m(2, 3) - m(1, 3) * m(2, 0) * m(3, 2) + m(1, 3) * m(3, 0) * m(2, 2));
1049 
1050  inv(1, 1) = (m(0, 0) * m(2, 2) * m(3, 3) - m(0, 0) * m(3, 2) * m(2, 3) - m(0, 2) * m(2, 0) * m(3, 3)
1051  + m(0, 2) * m(3, 0) * m(2, 3) + m(0, 3) * m(2, 0) * m(3, 2) - m(0, 3) * m(3, 0) * m(2, 2));
1052 
1053  inv(1, 2) = (-m(0, 0) * m(1, 2) * m(3, 3) + m(0, 0) * m(3, 2) * m(1, 3) + m(0, 2) * m(1, 0) * m(3, 3)
1054  - m(0, 2) * m(3, 0) * m(1, 3) - m(0, 3) * m(1, 0) * m(3, 2) + m(0, 3) * m(3, 0) * m(1, 2));
1055 
1056  inv(1, 3) = (m(0, 0) * m(1, 2) * m(2, 3) - m(0, 0) * m(2, 2) * m(1, 3) - m(0, 2) * m(1, 0) * m(2, 3)
1057  + m(0, 2) * m(2, 0) * m(1, 3) + m(0, 3) * m(1, 0) * m(2, 2) - m(0, 3) * m(2, 0) * m(1, 2));
1058 
1059  inv(2, 0) = (m(1, 0) * m(2, 1) * m(3, 3) - m(1, 0) * m(3, 1) * m(2, 3) - m(1, 1) * m(2, 0) * m(3, 3)
1060  + m(1, 1) * m(3, 0) * m(2, 3) + m(1, 3) * m(2, 0) * m(3, 1) - m(1, 3) * m(3, 0) * m(2, 1));
1061 
1062  inv(2, 1) = (-m(0, 0) * m(2, 1) * m(3, 3) + m(0, 0) * m(3, 1) * m(2, 3) + m(0, 1) * m(2, 0) * m(3, 3)
1063  - m(0, 1) * m(3, 0) * m(2, 3) - m(0, 3) * m(2, 0) * m(3, 1) + m(0, 3) * m(3, 0) * m(2, 1));
1064 
1065  inv(2, 2) = (m(0, 0) * m(1, 1) * m(3, 3) - m(0, 0) * m(3, 1) * m(1, 3) - m(0, 1) * m(1, 0) * m(3, 3)
1066  + m(0, 1) * m(3, 0) * m(1, 3) + m(0, 3) * m(1, 0) * m(3, 1) - m(0, 3) * m(3, 0) * m(1, 1));
1067 
1068  inv(2, 3) = (-m(0, 0) * m(1, 1) * m(2, 3) + m(0, 0) * m(2, 1) * m(1, 3) + m(0, 1) * m(1, 0) * m(2, 3)
1069  - m(0, 1) * m(2, 0) * m(1, 3) - m(0, 3) * m(1, 0) * m(2, 1) + m(0, 3) * m(2, 0) * m(1, 1));
1070 
1071  inv(3, 0) = (-m(1, 0) * m(2, 1) * m(3, 2) + m(1, 0) * m(3, 1) * m(2, 2) + m(1, 1) * m(2, 0) * m(3, 2)
1072  - m(1, 1) * m(3, 0) * m(2, 2) - m(1, 2) * m(2, 0) * m(3, 1) + m(1, 2) * m(3, 0) * m(2, 1));
1073 
1074  inv(3, 1) = (m(0, 0) * m(2, 1) * m(3, 2) - m(0, 0) * m(3, 1) * m(2, 2) - m(0, 1) * m(2, 0) * m(3, 2)
1075  + m(0, 1) * m(3, 0) * m(2, 2) + m(0, 2) * m(2, 0) * m(3, 1) - m(0, 2) * m(3, 0) * m(2, 1));
1076 
1077  inv(3, 2) = (-m(0, 0) * m(1, 1) * m(3, 2) + m(0, 0) * m(3, 1) * m(1, 2) + m(0, 1) * m(1, 0) * m(3, 2)
1078  - m(0, 1) * m(3, 0) * m(1, 2) - m(0, 2) * m(1, 0) * m(3, 1) + m(0, 2) * m(3, 0) * m(1, 1));
1079 
1080  inv(3, 3) = (m(0, 0) * m(1, 1) * m(2, 2) - m(0, 0) * m(2, 1) * m(1, 2) - m(0, 1) * m(1, 0) * m(2, 2)
1081  + m(0, 1) * m(2, 0) * m(1, 2) + m(0, 2) * m(1, 0) * m(2, 1) - m(0, 2) * m(2, 0) * m(1, 1));
1082 
1083  T det = m(0, 0) * inv(0, 0) + m(1, 0) * inv(0, 1) + m(2, 0) * inv(0, 2) + m(3, 0) * inv(0, 3);
1084 
1085  T inv_det = 1.0 / det;
1086 
1087  return inv_det * inv;
1088 }
1089 
1090 template <size_t nRows, typename T> CU_HOST_DEVICE Matrix<nRows, nRows, T> inverse(const Matrix<nRows, nRows, T>& mat)
1091 {
1093  T d = det(mat);
1094  if(d == 0)
1095  return ans;
1096  for(size_t i = 0; i < nRows; ++i)
1097  for(size_t j = 0; j < nRows; ++j) {
1098  ans(i, j) = cofactor(mat, j, i) / d;
1099  }
1100  return ans;
1101 }
1102 
1103 template <size_t nRows, size_t nCols, typename T>
1104 CU_HOST_DEVICE Matrix<nCols, nRows, T> transpose(const Matrix<nRows, nCols, T>& mat)
1105 {
1107  for(size_t i = 0; i < nRows; i++)
1108  for(size_t j = 0; j < nCols; j++)
1109  ans(j, i) = mat(i, j);
1110  return ans;
1111 }
1112 
1113 template <size_t nRows, size_t nCols, typename T> CU_HOST_DEVICE T normsq(const Matrix<nRows, nCols, T>& mat)
1114 {
1115  T acc = 0;
1116  for(size_t i = 0; i < nRows; i++)
1117  for(size_t j = 0; j < nCols; j++) {
1118  const T& v = mat(i, j);
1119  acc += v * v;
1120  }
1121  return acc;
1122 }
1123 
1124 template <size_t nRows, size_t nCols, typename T> CU_HOST_DEVICE T norm(const Matrix<nRows, nCols, T>& mat)
1125 {
1126  return std::sqrt(normsq(mat));
1127 }
1128 
1129 }}
1130 #endif
CU_HOST_DEVICE Matrix & zero(void)
Set the matrix to all zero.
Definition: Matrix.hpp:434
CU_HOST_DEVICE void t(const T &v)
Short access to the fourth element.
Definition: Vector.hpp:678
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(T(*fct)(const T &), const Matrix< nRows, nCols, T > &m)
Apply a unary function to each element of the matrix.
Definition: Matrix.hpp:712
CU_HOST_DEVICE Matrix(const T1 *values, MatrixLayout layout=C_STYLE)
Fill in the matrix with the values.
Definition: Matrix.hpp:93
CU_HOST_DEVICE const T * c_data() const
Returns a constant raw pointer on the data.
Definition: Vector.hpp:262
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(T(*fct)(const T1 &, const T2 &), const Matrix< nRows, nCols, T1 > &m1, const Matrix< nRows, nCols, T2 > &m2)
Apply a binary function to each element of the matrix.
Definition: Matrix.hpp:872
CU_HOST_DEVICE friend Matrix operator*(const T &scalar, const Matrix &mat)
Matrix-scalar multiplication.
Definition: Matrix.hpp:270
CU_HOST_DEVICE T operator()(size_t i, size_t j) const
Return the value at row i, column j.
Definition: Matrix.hpp:416
CU_HOST_DEVICE Matrix operator-(void) const
Matrix subtraction.
Definition: Matrix.hpp:200
CU_HOST_DEVICE Matrix & operator=(const T &value)
Set the matrix to a diagonal matrix.
Definition: Matrix.hpp:448
static CU_HOST_DEVICE Matrix< 4, 4, T > rotation(const Vector< 4, T > &direction, T angle)
Creates the 4x4 matrix corresponding to a rotation.
Definition: Matrix.hpp:502
CU_HOST_DEVICE T trace() const
Trace of the matrix.
Definition: Matrix.hpp:521
CU_HOST_DEVICE void x(const T &v)
Short access to the first element.
Definition: Vector.hpp:651
CU_HOST_DEVICE const T * c_data() const
Returns a constant raw pointer on the data.
Definition: Matrix.hpp:182
CU_HOST_DEVICE Vector< nCols, T > & operator[](size_t idx)
Returns the nth row.
Definition: Matrix.hpp:390
static CU_HOST_DEVICE Matrix< 3, 3, T > rotation(const Vector< 3, T > &direction, T angle)
Creates the 3x3 matrix corresponding to a rotation.
Definition: Matrix.hpp:481
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(const T &(*fct)(const T &, const T &), const Matrix< nRows, nCols, T > &m1, const Matrix< nRows, nCols, T > &m2)
Apply a binary function to each element of the matrix.
Definition: Matrix.hpp:830
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(T(*fct)(const T1 &), const Matrix< nRows, nCols, T1 > &m)
Apply a unary function to each element of the matrix.
Definition: Matrix.hpp:769
static CU_HOST_DEVICE size_t nbRows()
Returns the number of rows of the matrix.
Definition: Matrix.hpp:164
CU_HOST_DEVICE Matrix(const Vector< nCols, T1 > *vecs)
Fill the matrix with the array of vectors.
Definition: Matrix.hpp:78
CU_HOST_DEVICE Matrix(const Matrix< nRows, nCols, T1 > &mat)
Copy a matrix.
Definition: Matrix.hpp:66
static CU_HOST_DEVICE Vector< 2, size_t > size()
Returns the size of the matrix.
Definition: Matrix.hpp:156
CU_HOST_DEVICE T * data()
Returns a raw pointer on the data.
Definition: Matrix.hpp:192
static CU_HOST_DEVICE size_t nbColumns()
Returns the number of columns of the matrix.
Definition: Matrix.hpp:172
CU_HOST_DEVICE Matrix(const T &value)
Create a diagonal matrix.
Definition: Matrix.hpp:143
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(T(*fct)(T1), const Matrix< nRows, nCols, T1 > &m)
Apply a unary function to each element of the matrix.
Definition: Matrix.hpp:750
CU_HOST_DEVICE Matrix(const T *values, MatrixLayout layout=C_STYLE)
Fill in the matrix with the values.
Definition: Matrix.hpp:120
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(T(*fct)(T, T), const Matrix< nRows, nCols, T > &m1, const Matrix< nRows, nCols, T > &m2)
Apply a binary function to each element of the matrix.
Definition: Matrix.hpp:788
CU_HOST_DEVICE Matrix operator*(const T &scalar) const
Matrix-scalar multiplication.
Definition: Matrix.hpp:242
CU_HOST_DEVICE Matrix operator+(const Matrix &mat) const
Matrix addition.
Definition: Matrix.hpp:214
CU_HOST_DEVICE T * data()
Returns a raw pointer on the data.
Definition: Vector.hpp:224
CU_HOST_DEVICE void y(const T &v)
Short access to the second element.
Definition: Vector.hpp:660
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(const T &(*fct)(const T1 &), const Matrix< nRows, nCols, T1 > &m)
Apply a unary function to each element of the matrix.
Definition: Matrix.hpp:731
CU_HOST_DEVICE Vector< nRows, T > diag() const
Return the diagonal vector, if the matrix is square.
Definition: Matrix.hpp:550
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(const T &(*fct)(const T1 &, const T2 &), const Matrix< nRows, nCols, T1 > &m1, const Matrix< nRows, nCols, T2 > &m2)
Apply a binary function to each element of the matrix.
Definition: Matrix.hpp:893
CU_HOST_DEVICE Matrix< nCols, nRows, T > operator~()
Transpose the matrix.
Definition: Matrix.hpp:465
#define STATIC_ASSERT(B)
Assertion that works at compile time.
Definition: StaticAssert.hpp:47
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(const T &(*fct)(const T &), const Matrix< nRows, nCols, T > &m)
Apply a unary function to each element of the matrix.
Definition: Matrix.hpp:674
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(T(*fct)(T), const Matrix< nRows, nCols, T > &m)
Apply a unary function to each element of the matrix.
Definition: Matrix.hpp:693
static CU_HOST_DEVICE Matrix identity()
Returns an identity matrix.
Definition: Matrix.hpp:424
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(T(*fct)(T1, T2), const Matrix< nRows, nCols, T1 > &m1, const Matrix< nRows, nCols, T2 > &m2)
Apply a binary function to each element of the matrix.
Definition: Matrix.hpp:851
Define the STATIC_ASSERT macro.
Class representing a fixed-size matrix.
Definition: Matrix.hpp:34
Common definitions and utilities This file is shared by cuda, do not include headers that nvcc can't ...
CU_HOST_DEVICE T & operator()(size_t i, size_t j)
Return the value at row i, column j.
Definition: Matrix.hpp:408
CU_HOST_DEVICE void z(const T &v)
Short access to the third element.
Definition: Vector.hpp:669
CU_HOST_DEVICE Matrix operator/(const T &scalar) const
Matrix-scalar division.
Definition: Matrix.hpp:256
CU_HOST_DEVICE Matrix(void)
Create a matrix filled with 0s.
Definition: Matrix.hpp:54
CU_HOST_DEVICE Vector< nCols, T > operator[](size_t idx) const
Returns the nth row.
Definition: Matrix.hpp:400
CU_HOST_DEVICE Matrix< nRows, nCols, T > map(T(*fct)(const T &, const T &), const Matrix< nRows, nCols, T > &m1, const Matrix< nRows, nCols, T > &m2)
Apply a binary function to each element of the matrix.
Definition: Matrix.hpp:809
CU_HOST_DEVICE Vector< nRows, T > operator*(const Vector< nCols, T > &vec) const
Matrix*Column Vector.
Definition: Matrix.hpp:284
Defines the Vector class template This file is shared by cuda, do not include headers that nvcc can't...