https://mooseframework.inl.gov
ColumnMajorMatrix.h
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://mooseframework.inl.gov
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 "MooseError.h"
14 #include "ADReal.h"
15 #include "MooseTypes.h"
16 
17 #include "libmesh/type_tensor.h"
18 #include "libmesh/dense_vector.h"
19 
20 #include "metaphysicl/raw_type.h"
21 
22 // C++ includes
23 #include <iomanip>
24 
30 template <typename T>
32 {
33 public:
38  explicit ColumnMajorMatrixTempl(const unsigned int rows = Moose::dim,
39  const unsigned int cols = Moose::dim);
40 
45 
49  explicit ColumnMajorMatrixTempl(const TypeTensor<T> & tensor);
50 
51  explicit ColumnMajorMatrixTempl(const DenseMatrix<T> & rhs);
52 
53  explicit ColumnMajorMatrixTempl(const DenseVector<T> & rhs);
54 
58  ColumnMajorMatrixTempl(const TypeVector<T> & col1,
59  const TypeVector<T> & col2,
60  const TypeVector<T> & col3);
61 
66  unsigned int numEntries() const;
67 
72  void reshape(const unsigned int rows, const unsigned int cols);
73 
78  T & operator()(const unsigned int i, const unsigned int j = 0);
79 
86  T operator()(const unsigned int i, const unsigned int j = 0) const;
87 
91  void print();
92 
96  void print_scientific(std::ostream & os);
97 
101  void fill(TypeTensor<T> & tensor);
102 
106  void fill(DenseMatrix<T> & rhs);
107 
111  void fill(DenseVector<T> & rhs);
112 
118 
124 
130 
134  void setDiag(T value);
135 
139  void addDiag(T value);
140 
144  T tr() const;
145 
149  void zero();
150 
154  void identity();
155 
159  T doubleContraction(const ColumnMajorMatrixTempl<T> & rhs) const;
160 
164  T norm();
165 
169  unsigned int n() const;
170 
174  unsigned int m() const;
175 
179  void eigen(ColumnMajorMatrixTempl<T> & eval, ColumnMajorMatrixTempl<T> & evec) const;
180 
184  void eigenNonsym(ColumnMajorMatrixTempl<T> & eval_real,
185  ColumnMajorMatrixTempl<T> & eval_img,
186  ColumnMajorMatrixTempl<T> & evec_right,
187  ColumnMajorMatrixTempl<T> & eve_left) const;
188 
192  void exp(ColumnMajorMatrixTempl<T> & z) const;
193 
197  void inverse(ColumnMajorMatrixTempl<T> & invA) const;
198 
202  T * rawData();
203  const T * rawData() const;
204 
209 
214  ColumnMajorMatrixTempl<T> & operator=(const TypeTensor<T> & rhs);
215 
220  ColumnMajorMatrixTempl<T> & operator=(const DenseMatrix<T> & rhs);
221 
226  ColumnMajorMatrixTempl<T> & operator=(const DenseVector<T> & rhs);
227 
232  template <typename T2>
234 
239 
243  ColumnMajorMatrixTempl<T> operator*(T scalar) const;
244 
248  ColumnMajorMatrixTempl<T> operator*(const TypeVector<T> & rhs) const;
249 
250  // /**
251  // * Matrix Vector Multiplication of the TypeTensor Product. Note that the
252  // * Tensor type is treated as a single dimension Vector for this operation
253  // */
254  // ColumnMajorMatrixTempl operator*(const TypeTensor<T> & rhs) const;
255 
260 
265 
270 
278 
282  ColumnMajorMatrixTempl<T> & operator+=(const TypeTensor<T> & rhs);
283 
291 
295  ColumnMajorMatrixTempl<T> operator+(T scalar) const;
296 
301 
306 
311 
315  void checkSquareness() const;
316 
320  void checkShapeEquality(const ColumnMajorMatrixTempl<T> & rhs) const;
321 
325  bool operator==(const ColumnMajorMatrixTempl<T> & rhs) const;
326  bool operator!=(const ColumnMajorMatrixTempl<T> & rhs) const;
327 
328 protected:
329  unsigned int _n_rows, _n_cols, _n_entries;
330  std::vector<T> _values;
331 
332  template <typename T2>
333  friend void dataStore(std::ostream &, ColumnMajorMatrixTempl<T2> &, void *);
334  template <typename T2>
335  friend void dataLoad(std::istream &, ColumnMajorMatrixTempl<T2> &, void *);
336 };
337 
338 template <typename T>
339 inline unsigned int
341 {
342  return _n_entries;
343 }
344 
345 template <typename T>
346 inline void
347 ColumnMajorMatrixTempl<T>::reshape(unsigned int rows, unsigned int cols)
348 {
349  if (cols * rows == _n_entries)
350  {
351  _n_rows = rows;
352  _n_cols = cols;
353  }
354  else
355  {
356  _n_rows = rows;
357  _n_cols = cols;
358  _n_entries = _n_rows * _n_cols;
359  _values.resize(_n_entries);
360  }
361 }
362 
363 template <typename T>
364 inline T &
365 ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j)
366 {
367  if ((i * j) >= _n_entries)
368  mooseError("Reference outside of ColumnMajorMatrix bounds!");
369 
370  // Row major indexing!
371  return _values[(j * _n_rows) + i];
372 }
373 
374 template <typename T>
375 inline T
376 ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j) const
377 {
378  if ((i * j) >= _n_entries)
379  mooseError("Reference outside of ColumnMajorMatrix bounds!");
380 
381  // Row major indexing!
382  return _values[(j * _n_rows) + i];
383 }
384 
385 template <typename T>
386 inline void
388 {
389  ColumnMajorMatrixTempl<T> & s = (*this);
390 
391  for (unsigned int i = 0; i < _n_rows; i++)
392  {
393  for (unsigned int j = 0; j < _n_cols; j++)
394  Moose::out << std::setw(15) << s(i, j) << " ";
395 
396  Moose::out << std::endl;
397  }
398 }
399 
400 template <typename T>
401 inline void
403 {
404  ColumnMajorMatrixTempl<T> & s = (*this);
405 
406  for (unsigned int i = 0; i < _n_rows; i++)
407  {
408  for (unsigned int j = 0; j < _n_cols; j++)
409  os << std::setw(15) << std::scientific << std::setprecision(8) << s(i, j) << " ";
410 
411  os << std::endl;
412  }
413 }
414 
415 template <typename T>
416 inline void
417 ColumnMajorMatrixTempl<T>::fill(TypeTensor<T> & tensor)
418 {
419  if (Moose::dim * Moose::dim != _n_entries)
420  mooseError(
421  "Cannot fill tensor! The ColumnMajorMatrix doesn't have the same number of entries!");
422 
423  for (const auto j : libMesh::make_range(Moose::dim))
424  for (const auto i : libMesh::make_range(Moose::dim))
425  tensor(i, j) = _values[j * Moose::dim + i];
426 }
427 
428 template <typename T>
429 inline void
430 ColumnMajorMatrixTempl<T>::fill(DenseMatrix<T> & rhs)
431 {
432  if (rhs.n() * rhs.m() != _n_entries)
433  mooseError(
434  "Cannot fill dense matrix! The ColumnMajorMatrix doesn't have the same number of entries!");
435 
436  for (unsigned int j = 0, index = 0; j < rhs.m(); ++j)
437  for (unsigned int i = 0; i < rhs.n(); ++i, ++index)
438  rhs(i, j) = _values[index];
439 }
440 
441 template <typename T>
442 inline void
443 ColumnMajorMatrixTempl<T>::fill(DenseVector<T> & rhs)
444 {
445  if (_n_rows != rhs.size() || _n_cols != 1)
446  mooseError("ColumnMajorMatrix and DenseVector must be the same shape for a fill!");
447 
448  for (unsigned int i = 0; i < _n_rows; ++i)
449  rhs(i) = (*this)(i);
450 }
451 
452 template <typename T>
455 {
456  const ColumnMajorMatrixTempl<T> & s = (*this);
457 
458  ColumnMajorMatrixTempl<T> ret_matrix(_n_cols, _n_rows);
459 
460  for (unsigned int i = 0; i < _n_rows; i++)
461  for (unsigned int j = 0; j < _n_cols; j++)
462  ret_matrix(j, i) = s(i, j);
463 
464  return ret_matrix;
465 }
466 
467 template <typename T>
470 {
471  ColumnMajorMatrixTempl<T> & s = (*this);
472 
473  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols), I(_n_rows, _n_cols);
474 
475  I.identity();
476 
477  for (unsigned int i = 0; i < _n_rows; i++)
478  for (unsigned int j = 0; j < _n_cols; j++)
479  ret_matrix(i, j) = s(i, j) - I(i, j) * (s.tr() / 3.0);
480 
481  return ret_matrix;
482 }
483 
484 template <typename T>
485 inline void
487 {
488  this->checkSquareness();
489 
490  for (unsigned int i = 0; i < _n_rows; i++)
491  (*this)(i, i) = value;
492 }
493 
494 template <typename T>
495 inline void
497 {
498  this->checkSquareness();
499 
500  for (unsigned int i = 0; i < _n_rows; i++)
501  (*this)(i, i) += value;
502 }
503 
504 template <typename T>
505 inline T
507 {
508  this->checkSquareness();
509 
510  T trace = 0;
511 
512  for (unsigned int i = 0; i < _n_rows; i++)
513  trace += (*this)(i, i);
514 
515  return trace;
516 }
517 
518 template <typename T>
519 inline void
521 {
522  for (unsigned int i = 0; i < _n_entries; i++)
523  _values[i] = 0;
524 }
525 
526 template <typename T>
527 inline void
529 {
530  this->checkSquareness();
531 
532  zero();
533 
534  for (unsigned int i = 0; i < _n_rows; i++)
535  (*this)(i, i) = 1;
536 }
537 
538 template <typename T>
539 inline T
541 {
542  this->checkShapeEquality(rhs);
543 
544  T value = 0;
545 
546  for (unsigned int j = 0; j < _n_cols; j++)
547  for (unsigned int i = 0; i < _n_rows; i++)
548  value += (*this)(i, j) * rhs(i, j);
549 
550  return value;
551 }
552 
553 template <typename T>
554 inline unsigned int
556 {
557  return _n_rows;
558 }
559 
560 template <typename T>
561 inline unsigned int
563 {
564  return _n_cols;
565 }
566 
567 template <typename T>
568 inline T *
570 {
571  return &_values[0];
572 }
573 
574 template <typename T>
575 inline const T *
577 {
578  return &_values[0];
579 }
580 
581 template <typename T>
583 ColumnMajorMatrixTempl<T>::operator=(const TypeTensor<T> & rhs)
584 {
585  // Resize the tensor if necessary
586  if ((Moose::dim * Moose::dim) != _n_entries)
587  {
588  _n_entries = Moose::dim * Moose::dim;
589  _values.resize(_n_entries);
590  }
591 
592  // Make sure the shape is correct
593  reshape(Moose::dim, Moose::dim);
594 
595  ColumnMajorMatrixTempl<T> & s = (*this);
596 
597  // Copy the values
598  for (unsigned int j = 0; j < _n_cols; j++)
599  for (unsigned int i = 0; i < _n_cols; i++)
600  s(i, j) = rhs(i, j);
601 
602  return *this;
603 }
604 
605 template <typename T>
606 template <typename T2>
609 {
610  this->reshape(rhs.m(), rhs.n());
611 
612  for (MooseIndex(rhs.m()) i = 0; i < rhs.m(); ++i)
613  for (MooseIndex(rhs.n()) j = 0; j < rhs.n(); ++j)
614  (*this)(i, j) = rhs(i, j);
615 
616  return *this;
617 }
618 
619 template <typename T>
622 {
623  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
624 
625  for (unsigned int i = 0; i < _n_entries; i++)
626  ret_matrix._values[i] = _values[i] * scalar;
627 
628  return ret_matrix;
629 }
630 
631 template <typename T>
633 ColumnMajorMatrixTempl<T>::operator*(const TypeVector<T> & rhs) const
634 {
635  if (_n_cols != Moose::dim)
636  mooseError("Cannot perform matvec operation! The column dimension of "
637  "the ColumnMajorMatrix does not match the TypeVector!");
638 
639  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
640 
641  for (unsigned int i = 0; i < _n_rows; ++i)
642  for (unsigned int j = 0; j < _n_cols; ++j)
643  ret_matrix._values[i] += (*this)(i, j) * rhs(j);
644 
645  return ret_matrix;
646 }
647 
648 // template <typename T>
649 // inline ColumnMajorMatrixTempl<T>
650 // ColumnMajorMatrixTempl<T>::operator*(const TypeTensor<T> & rhs) const
651 // {
652 // mooseAssert(_n_cols == LIBMESH_DIM*LIBMESH_DIM, "Cannot perform matvec operation! The
653 // ColumnMajorMatrixTempl<T> doesn't have the correct shape!");
654 
655 // ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
656 
657 // for (unsigned int i=0; i<_n_rows; ++i)
658 // for (unsigned int j=0; j<_n_cols; ++j)
659 // // Treat the Tensor as a column major column vector
660 // ret_matrix._values[i] += (*this)(i, j) * rhs(j%3, j/3);
661 
662 // return ret_matrix;
663 // }
664 
665 template <typename T>
668 {
669  if (_n_cols != rhs._n_rows)
670  mooseError(
671  "Cannot perform matrix multiply! The shapes of the two operands are not compatible!");
672 
673  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, rhs._n_cols);
674 
675  for (unsigned int i = 0; i < ret_matrix._n_rows; ++i)
676  for (unsigned int j = 0; j < ret_matrix._n_cols; ++j)
677  for (unsigned int k = 0; k < _n_cols; ++k)
678  ret_matrix(i, j) += (*this)(i, k) * rhs(k, j);
679 
680  return ret_matrix;
681 }
682 
683 template <typename T>
686 {
687  this->checkShapeEquality(rhs);
688 
689  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
690 
691  for (unsigned int i = 0; i < _n_entries; i++)
692  ret_matrix._values[i] = _values[i] + rhs._values[i];
693 
694  return ret_matrix;
695 }
696 
697 template <typename T>
700 {
701  this->checkShapeEquality(rhs);
702 
703  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
704 
705  for (unsigned int i = 0; i < _n_entries; i++)
706  ret_matrix._values[i] = _values[i] - rhs._values[i];
707 
708  return ret_matrix;
709 }
710 
711 template <typename T>
714 {
715  this->checkShapeEquality(rhs);
716 
717  for (unsigned int i = 0; i < _n_entries; i++)
718  _values[i] += rhs._values[i];
719 
720  return *this;
721 }
722 
723 template <typename T>
725 ColumnMajorMatrixTempl<T>::operator+=(const TypeTensor<T> & rhs)
726 {
727  if ((_n_rows != Moose::dim) || (_n_cols != Moose::dim))
728  mooseError("Cannot perform matrix addition and assignment! The shapes of the two operands are "
729  "not compatible!");
730 
731  for (const auto j : libMesh::make_range(Moose::dim))
732  for (const auto i : libMesh::make_range(Moose::dim))
733  (*this)(i, j) += rhs(i, j);
734 
735  return *this;
736 }
737 
738 template <typename T>
741 {
742  this->checkShapeEquality(rhs);
743 
744  for (unsigned int i = 0; i < _n_entries; i++)
745  _values[i] -= rhs._values[i];
746 
747  return *this;
748 }
749 
750 template <typename T>
753 {
754  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
755 
756  for (unsigned int i = 0; i < _n_entries; i++)
757  ret_matrix._values[i] = _values[i] + scalar;
758 
759  return ret_matrix;
760 }
761 
762 template <typename T>
765 {
766  for (unsigned int i = 0; i < _n_entries; i++)
767  _values[i] *= scalar;
768  return *this;
769 }
770 
771 template <typename T>
774 {
775  for (unsigned int i = 0; i < _n_entries; i++)
776  _values[i] /= scalar;
777  return *this;
778 }
779 
780 template <typename T>
783 {
784  for (unsigned int i = 0; i < _n_entries; i++)
785  _values[i] += scalar;
786  return *this;
787 }
788 
789 template <typename T>
790 inline bool
792 {
793  if (_n_entries != rhs._n_entries || _n_rows != rhs._n_rows || _n_cols != rhs._n_cols)
794  return false;
795  return std::equal(_values.begin(), _values.end(), rhs._values.begin());
796 }
797 
798 template <typename T>
799 inline bool
801 {
802  return !(*this == rhs);
803 }
804 
806 
807 namespace MetaPhysicL
808 {
809 template <typename T>
810 struct RawType<ColumnMajorMatrixTempl<T>>
811 {
813 
815  {
816  value_type ret(in.m(), in.n());
817  for (MooseIndex(in.m()) i = 0; i < in.m(); ++i)
818  for (MooseIndex(in.n()) j = 0; j < in.n(); ++j)
819  ret(i, j) = raw_value(in(i, j));
820 
821  return ret;
822  }
823 };
824 }
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:302
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: EigenADReal.h:73
static constexpr std::size_t dim
This is the dimension of all vector and tensor datastructures used in MOOSE.
Definition: Moose.h:154
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.