www.mooseframework.org
ColumnMajorMatrix.h
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #pragma once
11 
12 #include "DenseMatrix.h"
13 #include "Moose.h" // using namespace libMesh
14 #include "MooseError.h"
15 #include "DualReal.h"
16 #include "MooseTypes.h"
17 
18 #include "libmesh/type_tensor.h"
19 #include "libmesh/dense_vector.h"
20 
21 #include "metaphysicl/raw_type.h"
22 
23 // C++ includes
24 #include <iomanip>
25 
31 template <typename T>
33 {
34 public:
39  explicit ColumnMajorMatrixTempl(const unsigned int rows = Moose::dim,
40  const unsigned int cols = Moose::dim);
41 
46 
50  explicit ColumnMajorMatrixTempl(const TypeTensor<T> & tensor);
51 
52  explicit ColumnMajorMatrixTempl(const DenseMatrix<T> & rhs);
53 
54  explicit ColumnMajorMatrixTempl(const DenseVector<T> & rhs);
55 
59  ColumnMajorMatrixTempl(const TypeVector<T> & col1,
60  const TypeVector<T> & col2,
61  const TypeVector<T> & col3);
62 
67  unsigned int numEntries() const;
68 
73  void reshape(const unsigned int rows, const unsigned int cols);
74 
79  T & operator()(const unsigned int i, const unsigned int j = 0);
80 
87  T operator()(const unsigned int i, const unsigned int j = 0) const;
88 
92  void print();
93 
97  void print_scientific(std::ostream & os);
98 
102  void fill(TypeTensor<T> & tensor);
103 
107  void fill(DenseMatrix<T> & rhs);
108 
112  void fill(DenseVector<T> & rhs);
113 
119 
125 
131 
135  void setDiag(T value);
136 
140  void addDiag(T value);
141 
145  T tr() const;
146 
150  void zero();
151 
155  void identity();
156 
160  T doubleContraction(const ColumnMajorMatrixTempl<T> & rhs) const;
161 
165  T norm();
166 
170  unsigned int n() const;
171 
175  unsigned int m() const;
176 
180  void eigen(ColumnMajorMatrixTempl<T> & eval, ColumnMajorMatrixTempl<T> & evec) const;
181 
185  void eigenNonsym(ColumnMajorMatrixTempl<T> & eval_real,
186  ColumnMajorMatrixTempl<T> & eval_img,
187  ColumnMajorMatrixTempl<T> & evec_right,
188  ColumnMajorMatrixTempl<T> & eve_left) const;
189 
193  void exp(ColumnMajorMatrixTempl<T> & z) const;
194 
198  void inverse(ColumnMajorMatrixTempl<T> & invA) const;
199 
203  T * rawData();
204  const T * rawData() const;
205 
210 
215  ColumnMajorMatrixTempl<T> & operator=(const TypeTensor<T> & rhs);
216 
221  ColumnMajorMatrixTempl<T> & operator=(const DenseMatrix<T> & rhs);
222 
227  ColumnMajorMatrixTempl<T> & operator=(const DenseVector<T> & rhs);
228 
233  template <typename T2>
235 
240 
244  ColumnMajorMatrixTempl<T> operator*(T scalar) const;
245 
249  ColumnMajorMatrixTempl<T> operator*(const TypeVector<T> & rhs) const;
250 
251  // /**
252  // * Matrix Vector Multiplication of the TypeTensor Product. Note that the
253  // * Tensor type is treated as a single dimension Vector for this operation
254  // */
255  // ColumnMajorMatrixTempl operator*(const TypeTensor<T> & rhs) const;
256 
261 
266 
271 
279 
283  ColumnMajorMatrixTempl<T> & operator+=(const TypeTensor<T> & rhs);
284 
292 
296  ColumnMajorMatrixTempl<T> operator+(T scalar) const;
297 
302 
307 
312 
316  void checkSquareness() const;
317 
321  void checkShapeEquality(const ColumnMajorMatrixTempl<T> & rhs) const;
322 
326  bool operator==(const ColumnMajorMatrixTempl<T> & rhs) const;
327  bool operator!=(const ColumnMajorMatrixTempl<T> & rhs) const;
328 
329 protected:
330  unsigned int _n_rows, _n_cols, _n_entries;
331  std::vector<T> _values;
332 
333  template <typename T2>
334  friend void dataStore(std::ostream &, ColumnMajorMatrixTempl<T2> &, void *);
335  template <typename T2>
336  friend void dataLoad(std::istream &, ColumnMajorMatrixTempl<T2> &, void *);
337 };
338 
339 template <typename T>
340 inline unsigned int
342 {
343  return _n_entries;
344 }
345 
346 template <typename T>
347 inline void
348 ColumnMajorMatrixTempl<T>::reshape(unsigned int rows, unsigned int cols)
349 {
350  if (cols * rows == _n_entries)
351  {
352  _n_rows = rows;
353  _n_cols = cols;
354  }
355  else
356  {
357  _n_rows = rows;
358  _n_cols = cols;
359  _n_entries = _n_rows * _n_cols;
360  _values.resize(_n_entries);
361  }
362 }
363 
364 template <typename T>
365 inline T &
366 ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j)
367 {
368  if ((i * j) >= _n_entries)
369  mooseError("Reference outside of ColumnMajorMatrix bounds!");
370 
371  // Row major indexing!
372  return _values[(j * _n_rows) + i];
373 }
374 
375 template <typename T>
376 inline T
377 ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j) const
378 {
379  if ((i * j) >= _n_entries)
380  mooseError("Reference outside of ColumnMajorMatrix bounds!");
381 
382  // Row major indexing!
383  return _values[(j * _n_rows) + i];
384 }
385 
386 template <typename T>
387 inline void
389 {
390  ColumnMajorMatrixTempl<T> & s = (*this);
391 
392  for (unsigned int i = 0; i < _n_rows; i++)
393  {
394  for (unsigned int j = 0; j < _n_cols; j++)
395  Moose::out << std::setw(15) << s(i, j) << " ";
396 
397  Moose::out << std::endl;
398  }
399 }
400 
401 template <typename T>
402 inline void
404 {
405  ColumnMajorMatrixTempl<T> & s = (*this);
406 
407  for (unsigned int i = 0; i < _n_rows; i++)
408  {
409  for (unsigned int j = 0; j < _n_cols; j++)
410  os << std::setw(15) << std::scientific << std::setprecision(8) << s(i, j) << " ";
411 
412  os << std::endl;
413  }
414 }
415 
416 template <typename T>
417 inline void
418 ColumnMajorMatrixTempl<T>::fill(TypeTensor<T> & tensor)
419 {
420  if (Moose::dim * Moose::dim != _n_entries)
421  mooseError(
422  "Cannot fill tensor! The ColumnMajorMatrix doesn't have the same number of entries!");
423 
424  for (const auto j : make_range(Moose::dim))
425  for (const auto i : make_range(Moose::dim))
426  tensor(i, j) = _values[j * Moose::dim + i];
427 }
428 
429 template <typename T>
430 inline void
431 ColumnMajorMatrixTempl<T>::fill(DenseMatrix<T> & rhs)
432 {
433  if (rhs.n() * rhs.m() != _n_entries)
434  mooseError(
435  "Cannot fill dense matrix! The ColumnMajorMatrix doesn't have the same number of entries!");
436 
437  for (unsigned int j = 0, index = 0; j < rhs.m(); ++j)
438  for (unsigned int i = 0; i < rhs.n(); ++i, ++index)
439  rhs(i, j) = _values[index];
440 }
441 
442 template <typename T>
443 inline void
444 ColumnMajorMatrixTempl<T>::fill(DenseVector<T> & rhs)
445 {
446  if (_n_rows != rhs.size() || _n_cols != 1)
447  mooseError("ColumnMajorMatrix and DenseVector must be the same shape for a fill!");
448 
449  for (unsigned int i = 0; i < _n_rows; ++i)
450  rhs(i) = (*this)(i);
451 }
452 
453 template <typename T>
456 {
457  const ColumnMajorMatrixTempl<T> & s = (*this);
458 
459  ColumnMajorMatrixTempl<T> ret_matrix(_n_cols, _n_rows);
460 
461  for (unsigned int i = 0; i < _n_rows; i++)
462  for (unsigned int j = 0; j < _n_cols; j++)
463  ret_matrix(j, i) = s(i, j);
464 
465  return ret_matrix;
466 }
467 
468 template <typename T>
471 {
472  ColumnMajorMatrixTempl<T> & s = (*this);
473 
474  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols), I(_n_rows, _n_cols);
475 
476  I.identity();
477 
478  for (unsigned int i = 0; i < _n_rows; i++)
479  for (unsigned int j = 0; j < _n_cols; j++)
480  ret_matrix(i, j) = s(i, j) - I(i, j) * (s.tr() / 3.0);
481 
482  return ret_matrix;
483 }
484 
485 template <typename T>
486 inline void
488 {
489  this->checkSquareness();
490 
491  for (unsigned int i = 0; i < _n_rows; i++)
492  (*this)(i, i) = value;
493 }
494 
495 template <typename T>
496 inline void
498 {
499  this->checkSquareness();
500 
501  for (unsigned int i = 0; i < _n_rows; i++)
502  (*this)(i, i) += value;
503 }
504 
505 template <typename T>
506 inline T
508 {
509  this->checkSquareness();
510 
511  T trace = 0;
512 
513  for (unsigned int i = 0; i < _n_rows; i++)
514  trace += (*this)(i, i);
515 
516  return trace;
517 }
518 
519 template <typename T>
520 inline void
522 {
523  for (unsigned int i = 0; i < _n_entries; i++)
524  _values[i] = 0;
525 }
526 
527 template <typename T>
528 inline void
530 {
531  this->checkSquareness();
532 
533  zero();
534 
535  for (unsigned int i = 0; i < _n_rows; i++)
536  (*this)(i, i) = 1;
537 }
538 
539 template <typename T>
540 inline T
542 {
543  this->checkShapeEquality(rhs);
544 
545  T value = 0;
546 
547  for (unsigned int j = 0; j < _n_cols; j++)
548  for (unsigned int i = 0; i < _n_rows; i++)
549  value += (*this)(i, j) * rhs(i, j);
550 
551  return value;
552 }
553 
554 template <typename T>
555 inline unsigned int
557 {
558  return _n_rows;
559 }
560 
561 template <typename T>
562 inline unsigned int
564 {
565  return _n_cols;
566 }
567 
568 template <typename T>
569 inline T *
571 {
572  return &_values[0];
573 }
574 
575 template <typename T>
576 inline const T *
578 {
579  return &_values[0];
580 }
581 
582 template <typename T>
584 ColumnMajorMatrixTempl<T>::operator=(const TypeTensor<T> & rhs)
585 {
586  // Resize the tensor if necessary
587  if ((Moose::dim * Moose::dim) != _n_entries)
588  {
589  _n_entries = Moose::dim * Moose::dim;
590  _values.resize(_n_entries);
591  }
592 
593  // Make sure the shape is correct
594  reshape(Moose::dim, Moose::dim);
595 
596  ColumnMajorMatrixTempl<T> & s = (*this);
597 
598  // Copy the values
599  for (unsigned int j = 0; j < _n_cols; j++)
600  for (unsigned int i = 0; i < _n_cols; i++)
601  s(i, j) = rhs(i, j);
602 
603  return *this;
604 }
605 
606 template <typename T>
607 template <typename T2>
610 {
611  this->reshape(rhs.m(), rhs.n());
612 
613  for (MooseIndex(rhs.m()) i = 0; i < rhs.m(); ++i)
614  for (MooseIndex(rhs.n()) j = 0; j < rhs.n(); ++j)
615  (*this)(i, j) = rhs(i, j);
616 
617  return *this;
618 }
619 
620 template <typename T>
623 {
624  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
625 
626  for (unsigned int i = 0; i < _n_entries; i++)
627  ret_matrix._values[i] = _values[i] * scalar;
628 
629  return ret_matrix;
630 }
631 
632 template <typename T>
634 ColumnMajorMatrixTempl<T>::operator*(const TypeVector<T> & rhs) const
635 {
636  if (_n_cols != Moose::dim)
637  mooseError("Cannot perform matvec operation! The column dimension of "
638  "the ColumnMajorMatrix does not match the TypeVector!");
639 
640  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
641 
642  for (unsigned int i = 0; i < _n_rows; ++i)
643  for (unsigned int j = 0; j < _n_cols; ++j)
644  ret_matrix._values[i] += (*this)(i, j) * rhs(j);
645 
646  return ret_matrix;
647 }
648 
649 // template <typename T>
650 // inline ColumnMajorMatrixTempl<T>
651 // ColumnMajorMatrixTempl<T>::operator*(const TypeTensor<T> & rhs) const
652 // {
653 // mooseAssert(_n_cols == LIBMESH_DIM*LIBMESH_DIM, "Cannot perform matvec operation! The
654 // ColumnMajorMatrixTempl<T> doesn't have the correct shape!");
655 
656 // ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
657 
658 // for (unsigned int i=0; i<_n_rows; ++i)
659 // for (unsigned int j=0; j<_n_cols; ++j)
660 // // Treat the Tensor as a column major column vector
661 // ret_matrix._values[i] += (*this)(i, j) * rhs(j%3, j/3);
662 
663 // return ret_matrix;
664 // }
665 
666 template <typename T>
669 {
670  if (_n_cols != rhs._n_rows)
671  mooseError(
672  "Cannot perform matrix multiply! The shapes of the two operands are not compatible!");
673 
674  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, rhs._n_cols);
675 
676  for (unsigned int i = 0; i < ret_matrix._n_rows; ++i)
677  for (unsigned int j = 0; j < ret_matrix._n_cols; ++j)
678  for (unsigned int k = 0; k < _n_cols; ++k)
679  ret_matrix(i, j) += (*this)(i, k) * rhs(k, j);
680 
681  return ret_matrix;
682 }
683 
684 template <typename T>
687 {
688  this->checkShapeEquality(rhs);
689 
690  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
691 
692  for (unsigned int i = 0; i < _n_entries; i++)
693  ret_matrix._values[i] = _values[i] + rhs._values[i];
694 
695  return ret_matrix;
696 }
697 
698 template <typename T>
701 {
702  this->checkShapeEquality(rhs);
703 
704  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
705 
706  for (unsigned int i = 0; i < _n_entries; i++)
707  ret_matrix._values[i] = _values[i] - rhs._values[i];
708 
709  return ret_matrix;
710 }
711 
712 template <typename T>
715 {
716  this->checkShapeEquality(rhs);
717 
718  for (unsigned int i = 0; i < _n_entries; i++)
719  _values[i] += rhs._values[i];
720 
721  return *this;
722 }
723 
724 template <typename T>
726 ColumnMajorMatrixTempl<T>::operator+=(const TypeTensor<T> & rhs)
727 {
728  if ((_n_rows != Moose::dim) || (_n_cols != Moose::dim))
729  mooseError("Cannot perform matrix addition and assignment! The shapes of the two operands are "
730  "not compatible!");
731 
732  for (const auto j : make_range(Moose::dim))
733  for (const auto i : make_range(Moose::dim))
734  (*this)(i, j) += rhs(i, j);
735 
736  return *this;
737 }
738 
739 template <typename T>
742 {
743  this->checkShapeEquality(rhs);
744 
745  for (unsigned int i = 0; i < _n_entries; i++)
746  _values[i] -= rhs._values[i];
747 
748  return *this;
749 }
750 
751 template <typename T>
754 {
755  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
756 
757  for (unsigned int i = 0; i < _n_entries; i++)
758  ret_matrix._values[i] = _values[i] + scalar;
759 
760  return ret_matrix;
761 }
762 
763 template <typename T>
766 {
767  for (unsigned int i = 0; i < _n_entries; i++)
768  _values[i] *= scalar;
769  return *this;
770 }
771 
772 template <typename T>
775 {
776  for (unsigned int i = 0; i < _n_entries; i++)
777  _values[i] /= scalar;
778  return *this;
779 }
780 
781 template <typename T>
784 {
785  for (unsigned int i = 0; i < _n_entries; i++)
786  _values[i] += scalar;
787  return *this;
788 }
789 
790 template <typename T>
791 inline bool
793 {
794  if (_n_entries != rhs._n_entries || _n_rows != rhs._n_rows || _n_cols != rhs._n_cols)
795  return false;
796  return std::equal(_values.begin(), _values.end(), rhs._values.begin());
797 }
798 
799 template <typename T>
800 inline bool
802 {
803  return !(*this == rhs);
804 }
805 
807 
808 namespace MetaPhysicL
809 {
810 template <typename T>
811 struct RawType<ColumnMajorMatrixTempl<T>>
812 {
814 
816  {
817  value_type ret(in.m(), in.n());
818  for (MooseIndex(in.m()) i = 0; i < in.m(); ++i)
819  for (MooseIndex(in.n()) j = 0; j < in.n(); ++j)
820  ret(i, j) = raw_value(in(i, j));
821 
822  return ret;
823  }
824 };
825 }
void inverse(ColumnMajorMatrixTempl< T > &invA) const
Returns inverse of a general matrix.
T tr() const
The trace of the CMM.
ColumnMajorMatrixTempl< T > & operator/=(T scalar)
Scalar Division plus assignment.
friend void dataStore(std::ostream &, ColumnMajorMatrixTempl< T2 > &, void *)
void print_scientific(std::ostream &os)
Prints to file.
void identity()
Turn the matrix into an identity matrix.
This class defines a Tensor that can change its shape.
ColumnMajorMatrixTempl< T > & operator-=(const ColumnMajorMatrixTempl< T > &rhs)
Matrix Matrix Subtraction plus assignment.
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:284
void reshape(const unsigned int rows, const unsigned int cols)
Change the shape of the tensor.
ColumnMajorMatrixTempl< T > abs()
Returns a matrix that is the absolute value of the matrix this was called on.
T norm()
The Euclidean norm of the matrix.
auto raw_value(const Eigen::Map< T > &in)
Definition: ADReal.h:73
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:148
ColumnMajorMatrixTempl< T > kronecker(const ColumnMajorMatrixTempl< T > &rhs) const
Kronecker Product.
ColumnMajorMatrixTempl< T > deviatoric()
Returns a matrix that is the deviatoric of the matrix this was called on.
We need to instantiate the following CompareTypes to tell the compiler that ADReal is a subtype of Ch...
std::basic_ostream< charT, traits > * os
Definition: InfixIterator.h:33
const Number zero
void eigen(ColumnMajorMatrixTempl< T > &eval, ColumnMajorMatrixTempl< T > &evec) const
Returns eigen system solve for a symmetric real matrix.
void checkSquareness() const
Check if matrix is square.
void setDiag(T value)
Set the value of each of the diagonals to the passed in value.
T * rawData()
Returns a reference to the raw data pointer.
void addDiag(T value)
Add to each of the diagonals the passsed in value.
ColumnMajorMatrixTempl< T > operator*(T scalar) const
Scalar multiplication of the ColumnMajorMatrixTempl.
ColumnMajorMatrixTempl< T > & operator*=(T scalar)
Scalar Multiplication plus assignment.
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
void fill(TypeTensor< T > &tensor)
Fills the passed in tensor with the values from this tensor.
unsigned int m() const
Returns the number of columns.
ColumnMajorMatrixTempl< T > operator+(const ColumnMajorMatrixTempl< T > &rhs) const
Matrix Matrix Addition.
ColumnMajorMatrixTempl< T > operator-(const ColumnMajorMatrixTempl< T > &rhs) const
Matrix Matrix Subtraction.
ColumnMajorMatrixTempl< typename RawType< T >::value_type > value_type
T doubleContraction(const ColumnMajorMatrixTempl< T > &rhs) const
Double contraction of two matrices ie A : B = Sum(A_ab * B_ba)
bool operator!=(const ColumnMajorMatrixTempl< T > &rhs) const
ColumnMajorMatrixTempl< Real > ColumnMajorMatrix
void zero()
Zero the matrix.
ColumnMajorMatrixTempl< T > transpose() const
Returns a matrix that is the transpose of the matrix this was called on.
void eigenNonsym(ColumnMajorMatrixTempl< T > &eval_real, ColumnMajorMatrixTempl< T > &eval_img, ColumnMajorMatrixTempl< T > &evec_right, ColumnMajorMatrixTempl< T > &eve_left) const
Returns eigen system solve for a non-symmetric real matrix.
void print()
Print the tensor.
void checkShapeEquality(const ColumnMajorMatrixTempl< T > &rhs) const
Check if matrices are of same shape.
T & operator()(const unsigned int i, const unsigned int j=0)
Get the i,j entry j defaults to zero so you can use it as a column vector.
static value_type value(const ColumnMajorMatrixTempl< T > &in)
ColumnMajorMatrixTempl(const unsigned int rows=Moose::dim, const unsigned int cols=Moose::dim)
Constructor that sets an initial number of entries and shape.
std::vector< T > _values
friend void dataLoad(std::istream &, ColumnMajorMatrixTempl< T2 > &, void *)
IntRange< T > make_range(T beg, T end)
ColumnMajorMatrixTempl< T > & operator=(const TypeTensor< T > &rhs)
Sets the values in this tensor to the values on the RHS.
unsigned int n() const
Returns the number of rows.
void exp(ColumnMajorMatrixTempl< T > &z) const
Returns matrix that is the exponential of the matrix this was called on.
ColumnMajorMatrixTempl< T > & operator+=(const ColumnMajorMatrixTempl< T > &rhs)
Matrix Matrix Addition plus assignment.
unsigned int numEntries() const
The total number of entries in the Tensor.
bool operator==(const ColumnMajorMatrixTempl< T > &rhs) const
Equality operators.