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 // MOOSE includes
13 #include "Moose.h" // using namespace libMesh
14 #include "MooseError.h"
15 #include "DualReal.h"
16 
17 #include "libmesh/type_tensor.h"
18 #include "libmesh/dense_matrix.h"
19 #include "libmesh/dense_vector.h"
20 
21 // C++ includes
22 #include <iomanip>
23 
29 template <typename T>
31 {
32 public:
37  explicit ColumnMajorMatrixTempl(const unsigned int rows = LIBMESH_DIM,
38  const unsigned int cols = LIBMESH_DIM);
39 
44 
48  explicit ColumnMajorMatrixTempl(const TypeTensor<T> & tensor);
49 
50  explicit ColumnMajorMatrixTempl(const DenseMatrix<T> & rhs);
51 
52  explicit ColumnMajorMatrixTempl(const DenseVector<T> & rhs);
53 
57  ColumnMajorMatrixTempl(const TypeVector<T> & col1,
58  const TypeVector<T> & col2,
59  const TypeVector<T> & col3);
60 
65  unsigned int numEntries() const;
66 
71  void reshape(const unsigned int rows, const unsigned int cols);
72 
77  T & operator()(const unsigned int i, const unsigned int j = 0);
78 
85  T operator()(const unsigned int i, const unsigned int j = 0) const;
86 
90  void print();
91 
95  void print_scientific(std::ostream & os);
96 
100  void fill(TypeTensor<T> & tensor);
101 
105  void fill(DenseMatrix<T> & rhs);
106 
110  void fill(DenseVector<T> & rhs);
111 
117 
123 
129 
133  void setDiag(T value);
134 
138  void addDiag(T value);
139 
143  T tr() const;
144 
148  void zero();
149 
153  void identity();
154 
158  T doubleContraction(const ColumnMajorMatrixTempl<T> & rhs) const;
159 
163  T norm();
164 
168  unsigned int n() const;
169 
173  unsigned int m() const;
174 
178  void eigen(ColumnMajorMatrixTempl<T> & eval, ColumnMajorMatrixTempl<T> & evec) const;
179 
183  void eigenNonsym(ColumnMajorMatrixTempl<T> & eval_real,
184  ColumnMajorMatrixTempl<T> & eval_img,
185  ColumnMajorMatrixTempl<T> & evec_right,
186  ColumnMajorMatrixTempl<T> & eve_left) const;
187 
191  void exp(ColumnMajorMatrixTempl<T> & z) const;
192 
196  void inverse(ColumnMajorMatrixTempl<T> & invA) const;
197 
201  T * rawData();
202  const T * rawData() const;
203 
208 
213  ColumnMajorMatrixTempl<T> & operator=(const TypeTensor<T> & rhs);
214 
219  ColumnMajorMatrixTempl<T> & operator=(const DenseMatrix<T> & rhs);
220 
225  ColumnMajorMatrixTempl<T> & operator=(const DenseVector<T> & rhs);
226 
232 
236  ColumnMajorMatrixTempl<T> operator*(T scalar) const;
237 
241  ColumnMajorMatrixTempl<T> operator*(const TypeVector<T> & rhs) const;
242 
243  // /**
244  // * Matrix Vector Multiplication of the TypeTensor Product. Note that the
245  // * Tensor type is treated as a single dimension Vector for this operation
246  // */
247  // ColumnMajorMatrixTempl operator*(const TypeTensor<T> & rhs) const;
248 
253 
258 
263 
271 
275  ColumnMajorMatrixTempl<T> & operator+=(const TypeTensor<T> & rhs);
276 
284 
288  ColumnMajorMatrixTempl<T> operator+(T scalar) const;
289 
294 
299 
304 
308  void checkSquareness() const;
309 
313  void checkShapeEquality(const ColumnMajorMatrixTempl<T> & rhs) const;
314 
318  bool operator==(const ColumnMajorMatrixTempl<T> & rhs) const;
319  bool operator!=(const ColumnMajorMatrixTempl<T> & rhs) const;
320 
321 protected:
322  unsigned int _n_rows, _n_cols, _n_entries;
323  std::vector<T> _values;
324 };
325 
326 template <typename T>
327 inline unsigned int
329 {
330  return _n_entries;
331 }
332 
333 template <typename T>
334 inline void
335 ColumnMajorMatrixTempl<T>::reshape(unsigned int rows, unsigned int cols)
336 {
337  if (cols * rows == _n_entries)
338  {
339  _n_rows = rows;
340  _n_cols = cols;
341  }
342  else
343  {
344  _n_rows = rows;
345  _n_cols = cols;
346  _n_entries = _n_rows * _n_cols;
347  _values.resize(_n_entries);
348  }
349 }
350 
351 template <typename T>
352 inline T &
353 ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j)
354 {
355  if ((i * j) >= _n_entries)
356  mooseError("Reference outside of ColumnMajorMatrix bounds!");
357 
358  // Row major indexing!
359  return _values[(j * _n_rows) + i];
360 }
361 
362 template <typename T>
363 inline T
364 ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j) const
365 {
366  if ((i * j) >= _n_entries)
367  mooseError("Reference outside of ColumnMajorMatrix bounds!");
368 
369  // Row major indexing!
370  return _values[(j * _n_rows) + i];
371 }
372 
373 template <typename T>
374 inline void
376 {
377  ColumnMajorMatrixTempl<T> & s = (*this);
378 
379  for (unsigned int i = 0; i < _n_rows; i++)
380  {
381  for (unsigned int j = 0; j < _n_cols; j++)
382  Moose::out << std::setw(15) << s(i, j) << " ";
383 
384  Moose::out << std::endl;
385  }
386 }
387 
388 template <typename T>
389 inline void
391 {
392  ColumnMajorMatrixTempl<T> & s = (*this);
393 
394  for (unsigned int i = 0; i < _n_rows; i++)
395  {
396  for (unsigned int j = 0; j < _n_cols; j++)
397  os << std::setw(15) << std::scientific << std::setprecision(8) << s(i, j) << " ";
398 
399  os << std::endl;
400  }
401 }
402 
403 template <typename T>
404 inline void
405 ColumnMajorMatrixTempl<T>::fill(TypeTensor<T> & tensor)
406 {
407  if (LIBMESH_DIM * LIBMESH_DIM != _n_entries)
408  mooseError(
409  "Cannot fill tensor! The ColumnMajorMatrix doesn't have the same number of entries!");
410 
411  for (unsigned int j = 0, index = 0; j < LIBMESH_DIM; ++j)
412  for (unsigned int i = 0; i < LIBMESH_DIM; ++i, ++index)
413  tensor(i, j) = _values[index];
414 }
415 
416 template <typename T>
417 inline void
418 ColumnMajorMatrixTempl<T>::fill(DenseMatrix<T> & rhs)
419 {
420  if (rhs.n() * rhs.m() != _n_entries)
421  mooseError(
422  "Cannot fill dense matrix! The ColumnMajorMatrix doesn't have the same number of entries!");
423 
424  for (unsigned int j = 0, index = 0; j < rhs.m(); ++j)
425  for (unsigned int i = 0; i < rhs.n(); ++i, ++index)
426  rhs(i, j) = _values[index];
427 }
428 
429 template <typename T>
430 inline void
431 ColumnMajorMatrixTempl<T>::fill(DenseVector<T> & rhs)
432 {
433  if (_n_rows != rhs.size() || _n_cols != 1)
434  mooseError("ColumnMajorMatrix and DenseVector must be the same shape for a fill!");
435 
436  for (unsigned int i = 0; i < _n_rows; ++i)
437  rhs(i) = (*this)(i);
438 }
439 
440 template <typename T>
443 {
444  const ColumnMajorMatrixTempl<T> & s = (*this);
445 
446  ColumnMajorMatrixTempl<T> ret_matrix(_n_cols, _n_rows);
447 
448  for (unsigned int i = 0; i < _n_rows; i++)
449  for (unsigned int j = 0; j < _n_cols; j++)
450  ret_matrix(j, i) = s(i, j);
451 
452  return ret_matrix;
453 }
454 
455 template <typename T>
458 {
459  ColumnMajorMatrixTempl<T> & s = (*this);
460 
461  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols), I(_n_rows, _n_cols);
462 
463  I.identity();
464 
465  for (unsigned int i = 0; i < _n_rows; i++)
466  for (unsigned int j = 0; j < _n_cols; j++)
467  ret_matrix(i, j) = s(i, j) - I(i, j) * (s.tr() / 3.0);
468 
469  return ret_matrix;
470 }
471 
472 template <typename T>
473 inline void
475 {
476  this->checkSquareness();
477 
478  for (unsigned int i = 0; i < _n_rows; i++)
479  (*this)(i, i) = value;
480 }
481 
482 template <typename T>
483 inline void
485 {
486  this->checkSquareness();
487 
488  for (unsigned int i = 0; i < _n_rows; i++)
489  (*this)(i, i) += value;
490 }
491 
492 template <typename T>
493 inline T
495 {
496  this->checkSquareness();
497 
498  T trace = 0;
499 
500  for (unsigned int i = 0; i < _n_rows; i++)
501  trace += (*this)(i, i);
502 
503  return trace;
504 }
505 
506 template <typename T>
507 inline void
509 {
510  for (unsigned int i = 0; i < _n_entries; i++)
511  _values[i] = 0;
512 }
513 
514 template <typename T>
515 inline void
517 {
518  this->checkSquareness();
519 
520  zero();
521 
522  for (unsigned int i = 0; i < _n_rows; i++)
523  (*this)(i, i) = 1;
524 }
525 
526 template <typename T>
527 inline T
529 {
530  this->checkShapeEquality(rhs);
531 
532  T value = 0;
533 
534  for (unsigned int j = 0; j < _n_cols; j++)
535  for (unsigned int i = 0; i < _n_rows; i++)
536  value += (*this)(i, j) * rhs(i, j);
537 
538  return value;
539 }
540 
541 template <typename T>
542 inline unsigned int
544 {
545  return _n_rows;
546 }
547 
548 template <typename T>
549 inline unsigned int
551 {
552  return _n_cols;
553 }
554 
555 template <typename T>
556 inline T *
558 {
559  return &_values[0];
560 }
561 
562 template <typename T>
563 inline const T *
565 {
566  return &_values[0];
567 }
568 
569 template <typename T>
571 ColumnMajorMatrixTempl<T>::operator=(const TypeTensor<T> & rhs)
572 {
573  // Resize the tensor if necessary
574  if ((LIBMESH_DIM * LIBMESH_DIM) != _n_entries)
575  {
576  _n_entries = LIBMESH_DIM * LIBMESH_DIM;
577  _values.resize(_n_entries);
578  }
579 
580  // Make sure the shape is correct
581  reshape(LIBMESH_DIM, LIBMESH_DIM);
582 
583  ColumnMajorMatrixTempl<T> & s = (*this);
584 
585  // Copy the values
586  for (unsigned int j = 0; j < _n_cols; j++)
587  for (unsigned int i = 0; i < _n_cols; i++)
588  s(i, j) = rhs(i, j);
589 
590  return *this;
591 }
592 
593 template <typename T>
596 {
597  _n_rows = rhs._n_rows;
598  _n_cols = rhs._n_cols;
599  _n_entries = rhs._n_entries;
600 
601  _values.resize(_n_entries);
602 
603  for (unsigned int i = 0; i < _n_entries; i++)
604  _values[i] = rhs._values[i];
605 
606  return *this;
607 }
608 
609 template <typename T>
611 {
612  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
613 
614  for (unsigned int i = 0; i < _n_entries; i++)
615  ret_matrix._values[i] = _values[i] * scalar;
616 
617  return ret_matrix;
618 }
619 
620 template <typename T>
622 operator*(const TypeVector<T> & rhs) const
623 {
624  if (_n_cols != LIBMESH_DIM)
625  mooseError("Cannot perform matvec operation! The column dimension of "
626  "the ColumnMajorMatrix does not match the TypeVector!");
627 
628  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
629 
630  for (unsigned int i = 0; i < _n_rows; ++i)
631  for (unsigned int j = 0; j < _n_cols; ++j)
632  ret_matrix._values[i] += (*this)(i, j) * rhs(j);
633 
634  return ret_matrix;
635 }
636 
637 // template <typename T>
638 // inline ColumnMajorMatrixTempl<T>
639 // ColumnMajorMatrixTempl<T>::operator*(const TypeTensor<T> & rhs) const
640 // {
641 // mooseAssert(_n_cols == LIBMESH_DIM*LIBMESH_DIM, "Cannot perform matvec operation! The
642 // ColumnMajorMatrixTempl<T> doesn't have the correct shape!");
643 
644 // ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
645 
646 // for (unsigned int i=0; i<_n_rows; ++i)
647 // for (unsigned int j=0; j<_n_cols; ++j)
648 // // Treat the Tensor as a column major column vector
649 // ret_matrix._values[i] += (*this)(i, j) * rhs(j%3, j/3);
650 
651 // return ret_matrix;
652 // }
653 
654 template <typename T>
657 {
658  if (_n_cols != rhs._n_rows)
659  mooseError(
660  "Cannot perform matrix multiply! The shapes of the two operands are not compatible!");
661 
662  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, rhs._n_cols);
663 
664  for (unsigned int i = 0; i < ret_matrix._n_rows; ++i)
665  for (unsigned int j = 0; j < ret_matrix._n_cols; ++j)
666  for (unsigned int k = 0; k < _n_cols; ++k)
667  ret_matrix(i, j) += (*this)(i, k) * rhs(k, j);
668 
669  return ret_matrix;
670 }
671 
672 template <typename T>
675 {
676  this->checkShapeEquality(rhs);
677 
678  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
679 
680  for (unsigned int i = 0; i < _n_entries; i++)
681  ret_matrix._values[i] = _values[i] + rhs._values[i];
682 
683  return ret_matrix;
684 }
685 
686 template <typename T>
689 {
690  this->checkShapeEquality(rhs);
691 
692  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
693 
694  for (unsigned int i = 0; i < _n_entries; i++)
695  ret_matrix._values[i] = _values[i] - rhs._values[i];
696 
697  return ret_matrix;
698 }
699 
700 template <typename T>
703 {
704  this->checkShapeEquality(rhs);
705 
706  for (unsigned int i = 0; i < _n_entries; i++)
707  _values[i] += rhs._values[i];
708 
709  return *this;
710 }
711 
712 template <typename T>
714 ColumnMajorMatrixTempl<T>::operator+=(const TypeTensor<T> & rhs)
715 {
716  if ((_n_rows != LIBMESH_DIM) || (_n_cols != LIBMESH_DIM))
717  mooseError("Cannot perform matrix addition and assignment! The shapes of the two operands are "
718  "not compatible!");
719 
720  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
721  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
722  (*this)(i, j) += rhs(i, j);
723 
724  return *this;
725 }
726 
727 template <typename T>
730 {
731  this->checkShapeEquality(rhs);
732 
733  for (unsigned int i = 0; i < _n_entries; i++)
734  _values[i] -= rhs._values[i];
735 
736  return *this;
737 }
738 
739 template <typename T>
742 {
743  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
744 
745  for (unsigned int i = 0; i < _n_entries; i++)
746  ret_matrix._values[i] = _values[i] + scalar;
747 
748  return ret_matrix;
749 }
750 
751 template <typename T>
754 {
755  for (unsigned int i = 0; i < _n_entries; i++)
756  _values[i] *= scalar;
757  return *this;
758 }
759 
760 template <typename T>
763 {
764  for (unsigned int i = 0; i < _n_entries; i++)
765  _values[i] /= scalar;
766  return *this;
767 }
768 
769 template <typename T>
772 {
773  for (unsigned int i = 0; i < _n_entries; i++)
774  _values[i] += scalar;
775  return *this;
776 }
777 
778 template <typename T>
779 inline bool
781 {
782  if (_n_entries != rhs._n_entries || _n_rows != rhs._n_rows || _n_cols != rhs._n_cols)
783  return false;
784  return std::equal(_values.begin(), _values.end(), rhs._values.begin());
785 }
786 
787 template <typename T>
788 inline bool
790 {
791  return !(*this == rhs);
792 }
793 
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.
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:207
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.
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.
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.
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(const unsigned int rows=LIBMESH_DIM, const unsigned int cols=LIBMESH_DIM)
Constructor that sets an initial number of entries and shape.
ColumnMajorMatrixTempl< T > operator-(const ColumnMajorMatrixTempl< T > &rhs) const
Matrix Matrix Subtraction.
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.
std::vector< T > _values
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.