www.mooseframework.org
ColumnMajorMatrix.h
Go to the documentation of this file.
1 /****************************************************************/
2 /* DO NOT MODIFY THIS HEADER */
3 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
4 /* */
5 /* (c) 2010 Battelle Energy Alliance, LLC */
6 /* ALL RIGHTS RESERVED */
7 /* */
8 /* Prepared by Battelle Energy Alliance, LLC */
9 /* Under Contract No. DE-AC07-05ID14517 */
10 /* With the U. S. Department of Energy */
11 /* */
12 /* See COPYRIGHT for full restrictions */
13 /****************************************************************/
14 
15 #ifndef COLUMNMAJORMATRIX_H
16 #define COLUMNMAJORMATRIX_H
17 
18 // MOOSE includes
19 #include "Moose.h" // using namespace libMesh
20 #include "MooseError.h" // mooseAssert
21 
22 #include "libmesh/type_tensor.h"
23 #include "libmesh/dense_matrix.h"
24 #include "libmesh/dense_vector.h"
25 
26 // C++ includes
27 #include <iomanip>
28 
35 {
36 public:
41  explicit ColumnMajorMatrix(const unsigned int rows = LIBMESH_DIM,
42  const unsigned int cols = LIBMESH_DIM);
43 
48 
52  explicit ColumnMajorMatrix(const TypeTensor<Real> & tensor);
53 
54  explicit ColumnMajorMatrix(const DenseMatrix<Real> & rhs);
55 
56  explicit ColumnMajorMatrix(const DenseVector<Real> & rhs);
57 
61  ColumnMajorMatrix(const TypeVector<Real> & col1,
62  const TypeVector<Real> & col2,
63  const TypeVector<Real> & col3);
64 
69  unsigned int numEntries() const;
70 
75  void reshape(const unsigned int rows, const unsigned int cols);
76 
81  Real & operator()(const unsigned int i, const unsigned int j = 0);
82 
89  Real operator()(const unsigned int i, const unsigned int j = 0) const;
90 
94  void print();
95 
99  void print_scientific(std::ostream & os);
100 
104  void fill(TypeTensor<Real> & tensor);
105 
109  void fill(DenseMatrix<Real> & rhs);
110 
114  void fill(DenseVector<Real> & rhs);
115 
121 
127 
133 
137  void setDiag(Real value);
138 
142  void addDiag(Real value);
143 
147  Real tr() const;
148 
152  void zero();
153 
157  void identity();
158 
162  Real doubleContraction(const ColumnMajorMatrix & rhs) const;
163 
167  Real norm();
168 
172  unsigned int n() const;
173 
177  unsigned int m() const;
178 
182  void eigen(ColumnMajorMatrix & eval, ColumnMajorMatrix & evec) const;
183 
187  void eigenNonsym(ColumnMajorMatrix & eval_real,
188  ColumnMajorMatrix & eval_img,
189  ColumnMajorMatrix & evec_right,
190  ColumnMajorMatrix & eve_left) const;
191 
195  void exp(ColumnMajorMatrix & z) const;
196 
200  void inverse(ColumnMajorMatrix & invA) const;
201 
205  Real * rawData();
206  const Real * rawData() const;
207 
211  ColumnMajorMatrix kronecker(const ColumnMajorMatrix & rhs) const;
212 
217  ColumnMajorMatrix & operator=(const TypeTensor<Real> & rhs);
218 
223  ColumnMajorMatrix & operator=(const DenseMatrix<Real> & rhs);
224 
229  ColumnMajorMatrix & operator=(const DenseVector<Real> & rhs);
230 
236 
240  ColumnMajorMatrix operator*(Real scalar) const;
241 
245  ColumnMajorMatrix operator*(const TypeVector<Real> & rhs) const;
246 
247  // /**
248  // * Matrix Vector Multiplication of the TypeTensor Product. Note that the
249  // * Tensor type is treated as a single dimension Vector for this operation
250  // */
251  // ColumnMajorMatrix operator*(const TypeTensor<Real> & rhs) const;
252 
256  ColumnMajorMatrix operator*(const ColumnMajorMatrix & rhs) const;
257 
261  ColumnMajorMatrix operator+(const ColumnMajorMatrix & rhs) const;
262 
266  ColumnMajorMatrix operator-(const ColumnMajorMatrix & rhs) const;
267 
275 
279  ColumnMajorMatrix & operator+=(const TypeTensor<Real> & rhs);
280 
288 
292  ColumnMajorMatrix operator+(Real scalar) const;
293 
297  ColumnMajorMatrix & operator*=(Real scalar);
298 
302  ColumnMajorMatrix & operator/=(Real scalar);
303 
307  ColumnMajorMatrix & operator+=(Real scalar);
308 
312  bool operator==(const ColumnMajorMatrix & rhs) const;
313  bool operator!=(const ColumnMajorMatrix & rhs) const;
314 
315 protected:
316  unsigned int _n_rows, _n_cols, _n_entries;
317  std::vector<Real> _values;
318 };
319 
320 inline unsigned int
322 {
323  return _n_entries;
324 }
325 
326 inline void
327 ColumnMajorMatrix::reshape(unsigned int rows, unsigned int cols)
328 {
329  if (cols * rows == _n_entries)
330  {
331  _n_rows = rows;
332  _n_cols = cols;
333  }
334  else
335  {
336  _n_rows = rows;
337  _n_cols = cols;
339  _values.resize(_n_entries);
340  }
341 }
342 
343 inline Real &
344 ColumnMajorMatrix::operator()(const unsigned int i, const unsigned int j)
345 {
346  mooseAssert((i * j) < _n_entries, "Reference outside of ColumnMajorMatrix bounds!");
347 
348  // Row major indexing!
349  return _values[(j * _n_rows) + i];
350 }
351 
352 inline Real
353 ColumnMajorMatrix::operator()(const unsigned int i, const unsigned int j) const
354 {
355  mooseAssert((i * j) < _n_entries, "Reference outside of ColumnMajorMatrix bounds!");
356 
357  // Row major indexing!
358  return _values[(j * _n_rows) + i];
359 }
360 
361 inline void
363 {
364  ColumnMajorMatrix & s = (*this);
365 
366  for (unsigned int i = 0; i < _n_rows; i++)
367  {
368  for (unsigned int j = 0; j < _n_cols; j++)
369  Moose::out << std::setw(15) << s(i, j) << " ";
370 
371  Moose::out << std::endl;
372  }
373 }
374 
375 inline void
377 {
378  ColumnMajorMatrix & s = (*this);
379 
380  for (unsigned int i = 0; i < _n_rows; i++)
381  {
382  for (unsigned int j = 0; j < _n_cols; j++)
383  os << std::setw(15) << std::scientific << std::setprecision(8) << s(i, j) << " ";
384 
385  os << std::endl;
386  }
387 }
388 
389 inline void
390 ColumnMajorMatrix::fill(TypeTensor<Real> & tensor)
391 {
392  mooseAssert(
393  LIBMESH_DIM * LIBMESH_DIM == _n_entries,
394  "Cannot fill tensor! The ColumnMajorMatrix doesn't have the same number of entries!");
395 
396  for (unsigned int j = 0, index = 0; j < LIBMESH_DIM; ++j)
397  for (unsigned int i = 0; i < LIBMESH_DIM; ++i, ++index)
398  tensor(i, j) = _values[index];
399 }
400 
401 inline void
402 ColumnMajorMatrix::fill(DenseMatrix<Real> & rhs)
403 {
404  mooseAssert(
405  rhs.n() * rhs.m() == _n_entries,
406  "Cannot fill dense matrix! The ColumnMajorMatrix doesn't have the same number of entries!");
407 
408  for (unsigned int j = 0, index = 0; j < rhs.m(); ++j)
409  for (unsigned int i = 0; i < rhs.n(); ++i, ++index)
410  rhs(i, j) = _values[index];
411 }
412 
413 inline void
414 ColumnMajorMatrix::fill(DenseVector<Real> & rhs)
415 {
416  mooseAssert(_n_rows == rhs.size(), "Vectors must be the same shape for a fill!");
417 
418  for (unsigned int i = 0; i < _n_rows; ++i)
419  rhs(i) = (*this)(i);
420 }
421 
422 inline ColumnMajorMatrix
424 {
425  const ColumnMajorMatrix & s = (*this);
426 
427  ColumnMajorMatrix ret_matrix(_n_cols, _n_rows);
428 
429  for (unsigned int i = 0; i < _n_rows; i++)
430  for (unsigned int j = 0; j < _n_cols; j++)
431  ret_matrix(j, i) = s(i, j);
432 
433  return ret_matrix;
434 }
435 
436 inline ColumnMajorMatrix
438 {
439  ColumnMajorMatrix & s = (*this);
440 
441  ColumnMajorMatrix ret_matrix(_n_rows, _n_cols), I(_n_rows, _n_cols);
442 
443  I.identity();
444 
445  for (unsigned int i = 0; i < _n_rows; i++)
446  for (unsigned int j = 0; j < _n_cols; j++)
447  ret_matrix(i, j) = s(i, j) - I(i, j) * (s.tr() / 3.0);
448 
449  return ret_matrix;
450 }
451 
452 inline ColumnMajorMatrix
454 {
455  ColumnMajorMatrix & s = (*this);
456 
457  ColumnMajorMatrix ret_matrix(_n_rows, _n_cols);
458 
459  for (unsigned int j = 0; j < _n_cols; j++)
460  for (unsigned int i = 0; i < _n_rows; i++)
461  ret_matrix(i, j) = std::abs(s(i, j));
462 
463  return ret_matrix;
464 }
465 
466 inline void
468 {
469  mooseAssert(_n_rows == _n_cols, "Cannot set the diagonal of a non-square matrix!");
470 
471  for (unsigned int i = 0; i < _n_rows; i++)
472  (*this)(i, i) = value;
473 }
474 
475 inline void
477 {
478  mooseAssert(_n_rows == _n_cols, "Cannot add to the diagonal of a non-square matrix!");
479 
480  for (unsigned int i = 0; i < _n_rows; i++)
481  (*this)(i, i) += value;
482 }
483 
484 inline Real
486 {
487  mooseAssert(_n_rows == _n_cols, "Cannot find the trace of a non-square matrix!");
488 
489  Real trace = 0;
490 
491  for (unsigned int i = 0; i < _n_rows; i++)
492  trace += (*this)(i, i);
493 
494  return trace;
495 }
496 
497 inline void
499 {
500  for (unsigned int i = 0; i < _n_entries; i++)
501  _values[i] = 0;
502 }
503 
504 inline void
506 {
507  mooseAssert(_n_rows == _n_cols, "Cannot set the diagonal of a non-square matrix!");
508 
509  zero();
510 
511  for (unsigned int i = 0; i < _n_rows; i++)
512  (*this)(i, i) = 1;
513 }
514 
515 inline Real
517 {
518  mooseAssert(_n_rows == rhs._n_cols && _n_cols == rhs._n_rows,
519  "Matrices must be the same shape for a double contraction!");
520 
521  Real value = 0;
522 
523  for (unsigned int j = 0; j < _n_cols; j++)
524  for (unsigned int i = 0; i < _n_rows; i++)
525  value += (*this)(i, j) * rhs(i, j);
526 
527  return value;
528 }
529 
530 inline Real
532 {
533  return std::sqrt(doubleContraction(*this));
534 }
535 
536 inline unsigned int
538 {
539  return _n_rows;
540 }
541 
542 inline unsigned int
544 {
545  return _n_cols;
546 }
547 
548 inline Real *
550 {
551  return &_values[0];
552 }
553 
554 inline const Real *
556 {
557  return &_values[0];
558 }
559 
560 inline ColumnMajorMatrix &
561 ColumnMajorMatrix::operator=(const TypeTensor<Real> & rhs)
562 {
563  // Resize the tensor if necessary
564  if ((LIBMESH_DIM * LIBMESH_DIM) != _n_entries)
565  {
566  _n_entries = LIBMESH_DIM * LIBMESH_DIM;
567  _values.resize(_n_entries);
568  }
569 
570  // Make sure the shape is correct
571  reshape(LIBMESH_DIM, LIBMESH_DIM);
572 
573  ColumnMajorMatrix & s = (*this);
574 
575  // Copy the values
576  for (unsigned int j = 0; j < _n_cols; j++)
577  for (unsigned int i = 0; i < _n_cols; i++)
578  s(i, j) = rhs(i, j);
579 
580  return *this;
581 }
582 
583 inline ColumnMajorMatrix &
585 {
586  _n_rows = rhs._n_rows;
587  _n_cols = rhs._n_cols;
588  _n_entries = rhs._n_entries;
589 
590  _values.resize(_n_entries);
591 
592  for (unsigned int i = 0; i < _n_entries; i++)
593  _values[i] = rhs._values[i];
594 
595  return *this;
596 }
597 
599 {
600  ColumnMajorMatrix ret_matrix(_n_rows, _n_cols);
601 
602  for (unsigned int i = 0; i < _n_entries; i++)
603  ret_matrix._values[i] = _values[i] * scalar;
604 
605  return ret_matrix;
606 }
607 
608 inline ColumnMajorMatrix ColumnMajorMatrix::operator*(const TypeVector<Real> & rhs) const
609 {
610  mooseAssert(_n_cols == LIBMESH_DIM,
611  "Cannot perform matvec operation! The column dimension of "
612  "the ColumnMajorMatrix does not match the TypeVector!");
613 
614  ColumnMajorMatrix ret_matrix(_n_rows, 1);
615 
616  for (unsigned int i = 0; i < _n_rows; ++i)
617  for (unsigned int j = 0; j < _n_cols; ++j)
618  ret_matrix._values[i] += (*this)(i, j) * rhs(j);
619 
620  return ret_matrix;
621 }
622 
623 // inline ColumnMajorMatrix
624 // ColumnMajorMatrix::operator*(const TypeTensor<Real> & rhs) const
625 // {
626 // mooseAssert(_n_cols == LIBMESH_DIM*LIBMESH_DIM, "Cannot perform matvec operation! The
627 // ColumnMajorMatrix doesn't have the correct shape!");
628 
629 // ColumnMajorMatrix ret_matrix(_n_rows, 1);
630 
631 // for (unsigned int i=0; i<_n_rows; ++i)
632 // for (unsigned int j=0; j<_n_cols; ++j)
633 // // Treat the Tensor as a column major column vector
634 // ret_matrix._values[i] += (*this)(i, j) * rhs(j%3, j/3);
635 
636 // return ret_matrix;
637 // }
638 
640 {
641  mooseAssert(
642  _n_cols == rhs._n_rows,
643  "Cannot perform matrix multiply! The shapes of the two operands are not compatible!");
644 
645  ColumnMajorMatrix ret_matrix(_n_rows, rhs._n_cols);
646 
647  for (unsigned int i = 0; i < ret_matrix._n_rows; ++i)
648  for (unsigned int j = 0; j < ret_matrix._n_cols; ++j)
649  for (unsigned int k = 0; k < _n_cols; ++k)
650  ret_matrix(i, j) += (*this)(i, k) * rhs(k, j);
651 
652  return ret_matrix;
653 }
654 
655 inline ColumnMajorMatrix
657 {
658  mooseAssert(
659  (_n_rows == rhs._n_rows) && (_n_cols == rhs._n_cols),
660  "Cannot perform matrix addition! The shapes of the two operands are not compatible!");
661 
662  ColumnMajorMatrix ret_matrix(_n_rows, _n_cols);
663 
664  for (unsigned int i = 0; i < _n_entries; i++)
665  ret_matrix._values[i] = _values[i] + rhs._values[i];
666 
667  return ret_matrix;
668 }
669 
670 inline ColumnMajorMatrix
672 {
673  mooseAssert(
674  (_n_rows == rhs._n_rows) && (_n_cols == rhs._n_cols),
675  "Cannot perform matrix addition! The shapes of the two operands are not compatible!");
676 
677  ColumnMajorMatrix ret_matrix(_n_rows, _n_cols);
678 
679  for (unsigned int i = 0; i < _n_entries; i++)
680  ret_matrix._values[i] = _values[i] - rhs._values[i];
681 
682  return ret_matrix;
683 }
684 
685 inline ColumnMajorMatrix &
687 {
688  mooseAssert((_n_rows == rhs._n_rows) && (_n_cols == rhs._n_cols),
689  "Cannot perform matrix addition and assignment! The shapes of the two operands are "
690  "not compatible!");
691 
692  for (unsigned int i = 0; i < _n_entries; i++)
693  _values[i] += rhs._values[i];
694 
695  return *this;
696 }
697 
698 inline ColumnMajorMatrix &
699 ColumnMajorMatrix::operator+=(const TypeTensor<Real> & rhs)
700 {
701  mooseAssert((_n_rows == LIBMESH_DIM) && (_n_cols == LIBMESH_DIM),
702  "Cannot perform matrix addition and assignment! The shapes of the two operands are "
703  "not compatible!");
704 
705  for (unsigned int j = 0; j < LIBMESH_DIM; ++j)
706  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
707  (*this)(i, j) += rhs(i, j);
708 
709  return *this;
710 }
711 
712 inline ColumnMajorMatrix &
714 {
715  mooseAssert((_n_rows == rhs._n_rows) && (_n_cols == rhs._n_cols),
716  "Cannot perform matrix subtraction and assignment! The shapes of the two operands "
717  "are not compatible!");
718 
719  for (unsigned int i = 0; i < _n_entries; i++)
720  _values[i] -= rhs._values[i];
721 
722  return *this;
723 }
724 
725 inline ColumnMajorMatrix
727 {
728  ColumnMajorMatrix ret_matrix(_n_rows, _n_cols);
729 
730  for (unsigned int i = 0; i < _n_entries; i++)
731  ret_matrix._values[i] = _values[i] + scalar;
732 
733  return ret_matrix;
734 }
735 
736 inline ColumnMajorMatrix &
738 {
739  for (unsigned int i = 0; i < _n_entries; i++)
740  _values[i] *= scalar;
741  return *this;
742 }
743 
744 inline ColumnMajorMatrix &
746 {
747  for (unsigned int i = 0; i < _n_entries; i++)
748  _values[i] /= scalar;
749  return *this;
750 }
751 
752 inline ColumnMajorMatrix &
754 {
755  for (unsigned int i = 0; i < _n_entries; i++)
756  _values[i] += scalar;
757  return *this;
758 }
759 
760 inline bool
762 {
763  if (_n_entries != rhs._n_entries || _n_rows != rhs._n_rows || _n_cols != rhs._n_cols)
764  return false;
765  return std::equal(_values.begin(), _values.end(), rhs._values.begin());
766 }
767 
768 inline bool
770 {
771  return !(*this == rhs);
772 }
773 
774 #endif // COLUMNMAJORMATRIX_H
ColumnMajorMatrix kronecker(const ColumnMajorMatrix &rhs) const
Kronecker Product.
void zero()
Zero the matrix.
unsigned int n() const
Returns the number of rows.
Real * rawData()
Returns a reference to the raw data pointer.
Real tr() const
The trace of the CMM.
void addDiag(Real value)
Add to each of the diagonals the passsed in value.
bool operator==(const ColumnMajorMatrix &rhs) const
Equality operators.
void setDiag(Real value)
Set the value of each of the diagonals to the passed in value.
unsigned int _n_entries
This class defines a Tensor that can change its shape.
ColumnMajorMatrix operator*(Real scalar) const
Scalar multiplication of the ColumnMajorMatrix.
ColumnMajorMatrix & operator+=(const ColumnMajorMatrix &rhs)
Matrix Matrix Addition plus assignment.
ColumnMajorMatrix operator+(const ColumnMajorMatrix &rhs) const
Matrix Matrix Addition.
ColumnMajorMatrix transpose() const
Returns a matrix that is the transpose of the matrix this was called on.
unsigned int numEntries() const
The total number of entries in the Tensor.
ColumnMajorMatrix(const unsigned int rows=LIBMESH_DIM, const unsigned int cols=LIBMESH_DIM)
Constructor that sets an initial number of entries and shape.
void identity()
Turn the matrix into an identity matrix.
std::vector< Real > _values
ColumnMajorMatrix & operator=(const TypeTensor< Real > &rhs)
Sets the values in this tensor to the values on the RHS.
void fill(TypeTensor< Real > &tensor)
Fills the passed in tensor with the values from this tensor.
ColumnMajorMatrix operator-(const ColumnMajorMatrix &rhs) const
Matrix Matrix Subtraction.
Real & 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.
void reshape(const unsigned int rows, const unsigned int cols)
Change the shape of the tensor.
ColumnMajorMatrix & operator*=(Real scalar)
Scalar Multiplication plus assignment.
ColumnMajorMatrix deviatoric()
Returns a matrix that is the deviatoric of the matrix this was called on.
void eigenNonsym(ColumnMajorMatrix &eval_real, ColumnMajorMatrix &eval_img, ColumnMajorMatrix &evec_right, ColumnMajorMatrix &eve_left) const
Returns eigen system solve for a non-symmetric real matrix.
ColumnMajorMatrix & operator-=(const ColumnMajorMatrix &rhs)
Matrix Matrix Subtraction plus assignment.
bool operator!=(const ColumnMajorMatrix &rhs) const
void exp(ColumnMajorMatrix &z) const
Returns matrix that is the exponential of the matrix this was called on.
Real norm()
The Euclidean norm of the matrix.
void print()
Print the tensor.
ColumnMajorMatrix & operator/=(Real scalar)
Scalar Division plus assignment.
Real doubleContraction(const ColumnMajorMatrix &rhs) const
Double contraction of two matrices ie A : B = Sum(A_ab * B_ba)
void print_scientific(std::ostream &os)
Prints to file.
void inverse(ColumnMajorMatrix &invA) const
Returns inverse of a general matrix.
unsigned int m() const
Returns the number of columns.
void eigen(ColumnMajorMatrix &eval, ColumnMajorMatrix &evec) const
Returns eigen system solve for a symmetric real matrix.
ColumnMajorMatrix abs()
Returns a matrix that is the absolute value of the matrix this was called on.