ColumnMajorMatrix.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#pragma once
 
#include "DenseMatrix.h"
#include "MooseError.h"
#include "ADReal.h"
#include "MooseTypes.h"
 
#include "libmesh/type_tensor.h"
#include "libmesh/dense_vector.h"
 
#include "metaphysicl/raw_type.h"
 
// C++ includes
#include <iomanip>
 
template <typename T>
class ColumnMajorMatrixTempl
{
public:
  explicit ColumnMajorMatrixTempl(const unsigned int rows = Moose::dim,
                                  const unsigned int cols = Moose::dim);
 
  ColumnMajorMatrixTempl(const ColumnMajorMatrixTempl<T> & rhs);
 
  explicit ColumnMajorMatrixTempl(const TypeTensor<T> & tensor);
 
  explicit ColumnMajorMatrixTempl(const DenseMatrix<T> & rhs);
 
  explicit ColumnMajorMatrixTempl(const DenseVector<T> & rhs);
 
  ColumnMajorMatrixTempl(const TypeVector<T> & col1,
                         const TypeVector<T> & col2,
                         const TypeVector<T> & col3);
 
  unsigned int numEntries() const;
 
  void reshape(const unsigned int rows, const unsigned int cols);
 
  T & operator()(const unsigned int i, const unsigned int j = 0);
 
  T operator()(const unsigned int i, const unsigned int j = 0) const;
 
  void print();
 
  void print_scientific(std::ostream & os);
 
  void fill(TypeTensor<T> & tensor);
 
  void fill(DenseMatrix<T> & rhs);
 
  void fill(DenseVector<T> & rhs);
 
  ColumnMajorMatrixTempl<T> transpose() const;
 
  ColumnMajorMatrixTempl<T> deviatoric();
 
  ColumnMajorMatrixTempl<T> abs();
 
  void setDiag(T value);
 
  void addDiag(T value);
 
  T tr() const;
 
  void zero();
 
  void identity();
 
  T doubleContraction(const ColumnMajorMatrixTempl<T> & rhs) const;
 
  T norm();
 
  unsigned int n() const;
 
  unsigned int m() const;
 
  void eigen(ColumnMajorMatrixTempl<T> & eval, ColumnMajorMatrixTempl<T> & evec) const;
 
  void eigenNonsym(ColumnMajorMatrixTempl<T> & eval_real,
                   ColumnMajorMatrixTempl<T> & eval_img,
                   ColumnMajorMatrixTempl<T> & evec_right,
                   ColumnMajorMatrixTempl<T> & eve_left) const;
 
  void exp(ColumnMajorMatrixTempl<T> & z) const;
 
  void inverse(ColumnMajorMatrixTempl<T> & invA) const;
 
  T * rawData();
  const T * rawData() const;
 
  ColumnMajorMatrixTempl<T> kronecker(const ColumnMajorMatrixTempl<T> & rhs) const;
 
  ColumnMajorMatrixTempl<T> & operator=(const TypeTensor<T> & rhs);
 
  ColumnMajorMatrixTempl<T> & operator=(const DenseMatrix<T> & rhs);
 
  ColumnMajorMatrixTempl<T> & operator=(const DenseVector<T> & rhs);
 
  template <typename T2>
  ColumnMajorMatrixTempl<T> & operator=(const ColumnMajorMatrixTempl<T2> & rhs);
 
  ColumnMajorMatrixTempl<T> & operator=(const ColumnMajorMatrixTempl<T> & rhs) = default;
 
  ColumnMajorMatrixTempl<T> operator*(T scalar) const;
 
  ColumnMajorMatrixTempl<T> operator*(const TypeVector<T> & rhs) const;
 
  //   /**
  //    * Matrix Vector Multiplication of the TypeTensor Product.  Note that the
  //    * Tensor type is treated as a single dimension Vector for this operation
  //    */
  //   ColumnMajorMatrixTempl operator*(const TypeTensor<T> & rhs) const;
 
  ColumnMajorMatrixTempl<T> operator*(const ColumnMajorMatrixTempl<T> & rhs) const;
 
  ColumnMajorMatrixTempl<T> operator+(const ColumnMajorMatrixTempl<T> & rhs) const;
 
  ColumnMajorMatrixTempl<T> operator-(const ColumnMajorMatrixTempl<T> & rhs) const;
 
  ColumnMajorMatrixTempl<T> & operator+=(const ColumnMajorMatrixTempl<T> & rhs);
 
  ColumnMajorMatrixTempl<T> & operator+=(const TypeTensor<T> & rhs);
 
  ColumnMajorMatrixTempl<T> & operator-=(const ColumnMajorMatrixTempl<T> & rhs);
 
  ColumnMajorMatrixTempl<T> operator+(T scalar) const;
 
  ColumnMajorMatrixTempl<T> & operator*=(T scalar);
 
  ColumnMajorMatrixTempl<T> & operator/=(T scalar);
 
  ColumnMajorMatrixTempl<T> & operator+=(T scalar);
 
  void checkSquareness() const;
 
  void checkShapeEquality(const ColumnMajorMatrixTempl<T> & rhs) const;
 
  bool operator==(const ColumnMajorMatrixTempl<T> & rhs) const;
  bool operator!=(const ColumnMajorMatrixTempl<T> & rhs) const;
 
protected:
  unsigned int _n_rows, _n_cols, _n_entries;
  std::vector<T> _values;
 
  template <typename T2>
  friend void dataStore(std::ostream &, ColumnMajorMatrixTempl<T2> &, void *);
  template <typename T2>
  friend void dataLoad(std::istream &, ColumnMajorMatrixTempl<T2> &, void *);
};
 
template <typename T>
inline unsigned int
ColumnMajorMatrixTempl<T>::numEntries() const
{
  return _n_entries;
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::reshape(unsigned int rows, unsigned int cols)
{
  if (cols * rows == _n_entries)
  {
    _n_rows = rows;
    _n_cols = cols;
  }
  else
  {
    _n_rows = rows;
    _n_cols = cols;
    _n_entries = _n_rows * _n_cols;
    _values.resize(_n_entries);
  }
}
 
template <typename T>
inline T &
ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j)
{
  if ((i * j) >= _n_entries)
    mooseError("Reference outside of ColumnMajorMatrix bounds!");
 
  // Row major indexing!
  return _values[(j * _n_rows) + i];
}
 
template <typename T>
inline T
ColumnMajorMatrixTempl<T>::operator()(const unsigned int i, const unsigned int j) const
{
  if ((i * j) >= _n_entries)
    mooseError("Reference outside of ColumnMajorMatrix bounds!");
 
  // Row major indexing!
  return _values[(j * _n_rows) + i];
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::print()
{
  ColumnMajorMatrixTempl<T> & s = (*this);
 
  for (unsigned int i = 0; i < _n_rows; i++)
  {
    for (unsigned int j = 0; j < _n_cols; j++)
      Moose::out << std::setw(15) << s(i, j) << " ";
 
    Moose::out << std::endl;
  }
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::print_scientific(std::ostream & os)
{
  ColumnMajorMatrixTempl<T> & s = (*this);
 
  for (unsigned int i = 0; i < _n_rows; i++)
  {
    for (unsigned int j = 0; j < _n_cols; j++)
      os << std::setw(15) << std::scientific << std::setprecision(8) << s(i, j) << " ";
 
    os << std::endl;
  }
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::fill(TypeTensor<T> & tensor)
{
  if (Moose::dim * Moose::dim != _n_entries)
    mooseError(
        "Cannot fill tensor! The ColumnMajorMatrix doesn't have the same number of entries!");
 
  for (const auto j : libMesh::make_range(Moose::dim))
    for (const auto i : libMesh::make_range(Moose::dim))
      tensor(i, j) = _values[j * Moose::dim + i];
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::fill(DenseMatrix<T> & rhs)
{
  if (rhs.n() * rhs.m() != _n_entries)
    mooseError(
        "Cannot fill dense matrix! The ColumnMajorMatrix doesn't have the same number of entries!");
 
  for (unsigned int j = 0, index = 0; j < rhs.m(); ++j)
    for (unsigned int i = 0; i < rhs.n(); ++i, ++index)
      rhs(i, j) = _values[index];
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::fill(DenseVector<T> & rhs)
{
  if (_n_rows != rhs.size() || _n_cols != 1)
    mooseError("ColumnMajorMatrix and DenseVector must be the same shape for a fill!");
 
  for (unsigned int i = 0; i < _n_rows; ++i)
    rhs(i) = (*this)(i);
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T>
ColumnMajorMatrixTempl<T>::transpose() const
{
  const ColumnMajorMatrixTempl<T> & s = (*this);
 
  ColumnMajorMatrixTempl<T> ret_matrix(_n_cols, _n_rows);
 
  for (unsigned int i = 0; i < _n_rows; i++)
    for (unsigned int j = 0; j < _n_cols; j++)
      ret_matrix(j, i) = s(i, j);
 
  return ret_matrix;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T>
ColumnMajorMatrixTempl<T>::deviatoric()
{
  ColumnMajorMatrixTempl<T> & s = (*this);
 
  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols), I(_n_rows, _n_cols);
 
  I.identity();
 
  for (unsigned int i = 0; i < _n_rows; i++)
    for (unsigned int j = 0; j < _n_cols; j++)
      ret_matrix(i, j) = s(i, j) - I(i, j) * (s.tr() / 3.0);
 
  return ret_matrix;
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::setDiag(T value)
{
  this->checkSquareness();
 
  for (unsigned int i = 0; i < _n_rows; i++)
    (*this)(i, i) = value;
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::addDiag(T value)
{
  this->checkSquareness();
 
  for (unsigned int i = 0; i < _n_rows; i++)
    (*this)(i, i) += value;
}
 
template <typename T>
inline T
ColumnMajorMatrixTempl<T>::tr() const
{
  this->checkSquareness();
 
  T trace = 0;
 
  for (unsigned int i = 0; i < _n_rows; i++)
    trace += (*this)(i, i);
 
  return trace;
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::zero()
{
  for (unsigned int i = 0; i < _n_entries; i++)
    _values[i] = 0;
}
 
template <typename T>
inline void
ColumnMajorMatrixTempl<T>::identity()
{
  this->checkSquareness();
 
  zero();
 
  for (unsigned int i = 0; i < _n_rows; i++)
    (*this)(i, i) = 1;
}
 
template <typename T>
inline T
ColumnMajorMatrixTempl<T>::doubleContraction(const ColumnMajorMatrixTempl<T> & rhs) const
{
  this->checkShapeEquality(rhs);
 
  T value = 0;
 
  for (unsigned int j = 0; j < _n_cols; j++)
    for (unsigned int i = 0; i < _n_rows; i++)
      value += (*this)(i, j) * rhs(i, j);
 
  return value;
}
 
template <typename T>
inline unsigned int
ColumnMajorMatrixTempl<T>::n() const
{
  return _n_rows;
}
 
template <typename T>
inline unsigned int
ColumnMajorMatrixTempl<T>::m() const
{
  return _n_cols;
}
 
template <typename T>
inline T *
ColumnMajorMatrixTempl<T>::rawData()
{
  return &_values[0];
}
 
template <typename T>
inline const T *
ColumnMajorMatrixTempl<T>::rawData() const
{
  return &_values[0];
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T> &
ColumnMajorMatrixTempl<T>::operator=(const TypeTensor<T> & rhs)
{
  // Resize the tensor if necessary
  if ((Moose::dim * Moose::dim) != _n_entries)
  {
    _n_entries = Moose::dim * Moose::dim;
    _values.resize(_n_entries);
  }
 
  // Make sure the shape is correct
  reshape(Moose::dim, Moose::dim);
 
  ColumnMajorMatrixTempl<T> & s = (*this);
 
  // Copy the values
  for (unsigned int j = 0; j < _n_cols; j++)
    for (unsigned int i = 0; i < _n_cols; i++)
      s(i, j) = rhs(i, j);
 
  return *this;
}
 
template <typename T>
template <typename T2>
inline ColumnMajorMatrixTempl<T> &
ColumnMajorMatrixTempl<T>::operator=(const ColumnMajorMatrixTempl<T2> & rhs)
{
  this->reshape(rhs.m(), rhs.n());
 
  for (MooseIndex(rhs.m()) i = 0; i < rhs.m(); ++i)
    for (MooseIndex(rhs.n()) j = 0; j < rhs.n(); ++j)
      (*this)(i, j) = rhs(i, j);
 
  return *this;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T>
ColumnMajorMatrixTempl<T>::operator*(T scalar) const
{
  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
 
  for (unsigned int i = 0; i < _n_entries; i++)
    ret_matrix._values[i] = _values[i] * scalar;
 
  return ret_matrix;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T>
ColumnMajorMatrixTempl<T>::operator*(const TypeVector<T> & rhs) const
{
  if (_n_cols != Moose::dim)
    mooseError("Cannot perform matvec operation! The column dimension of "
               "the ColumnMajorMatrix does not match the TypeVector!");
 
  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
 
  for (unsigned int i = 0; i < _n_rows; ++i)
    for (unsigned int j = 0; j < _n_cols; ++j)
      ret_matrix._values[i] += (*this)(i, j) * rhs(j);
 
  return ret_matrix;
}
 
// template <typename T>
// inline ColumnMajorMatrixTempl<T>
// ColumnMajorMatrixTempl<T>::operator*(const TypeTensor<T> & rhs) const
// {
//   mooseAssert(_n_cols == LIBMESH_DIM*LIBMESH_DIM, "Cannot perform matvec operation!  The
//   ColumnMajorMatrixTempl<T> doesn't have the correct shape!");
 
//   ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, 1);
 
//   for (unsigned int i=0; i<_n_rows; ++i)
//     for (unsigned int j=0; j<_n_cols; ++j)
//       // Treat the Tensor as a column major column vector
//       ret_matrix._values[i] += (*this)(i, j) * rhs(j%3, j/3);
 
//   return ret_matrix;
// }
 
template <typename T>
inline ColumnMajorMatrixTempl<T>
ColumnMajorMatrixTempl<T>::operator*(const ColumnMajorMatrixTempl<T> & rhs) const
{
  if (_n_cols != rhs._n_rows)
    mooseError(
        "Cannot perform matrix multiply! The shapes of the two operands are not compatible!");
 
  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, rhs._n_cols);
 
  for (unsigned int i = 0; i < ret_matrix._n_rows; ++i)
    for (unsigned int j = 0; j < ret_matrix._n_cols; ++j)
      for (unsigned int k = 0; k < _n_cols; ++k)
        ret_matrix(i, j) += (*this)(i, k) * rhs(k, j);
 
  return ret_matrix;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T>
ColumnMajorMatrixTempl<T>::operator+(const ColumnMajorMatrixTempl<T> & rhs) const
{
  this->checkShapeEquality(rhs);
 
  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
 
  for (unsigned int i = 0; i < _n_entries; i++)
    ret_matrix._values[i] = _values[i] + rhs._values[i];
 
  return ret_matrix;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T>
ColumnMajorMatrixTempl<T>::operator-(const ColumnMajorMatrixTempl<T> & rhs) const
{
  this->checkShapeEquality(rhs);
 
  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
 
  for (unsigned int i = 0; i < _n_entries; i++)
    ret_matrix._values[i] = _values[i] - rhs._values[i];
 
  return ret_matrix;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T> &
ColumnMajorMatrixTempl<T>::operator+=(const ColumnMajorMatrixTempl<T> & rhs)
{
  this->checkShapeEquality(rhs);
 
  for (unsigned int i = 0; i < _n_entries; i++)
    _values[i] += rhs._values[i];
 
  return *this;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T> &
ColumnMajorMatrixTempl<T>::operator+=(const TypeTensor<T> & rhs)
{
  if ((_n_rows != Moose::dim) || (_n_cols != Moose::dim))
    mooseError("Cannot perform matrix addition and assignment! The shapes of the two operands are "
               "not compatible!");
 
  for (const auto j : libMesh::make_range(Moose::dim))
    for (const auto i : libMesh::make_range(Moose::dim))
      (*this)(i, j) += rhs(i, j);
 
  return *this;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T> &
ColumnMajorMatrixTempl<T>::operator-=(const ColumnMajorMatrixTempl<T> & rhs)
{
  this->checkShapeEquality(rhs);
 
  for (unsigned int i = 0; i < _n_entries; i++)
    _values[i] -= rhs._values[i];
 
  return *this;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T>
ColumnMajorMatrixTempl<T>::operator+(T scalar) const
{
  ColumnMajorMatrixTempl<T> ret_matrix(_n_rows, _n_cols);
 
  for (unsigned int i = 0; i < _n_entries; i++)
    ret_matrix._values[i] = _values[i] + scalar;
 
  return ret_matrix;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T> &
ColumnMajorMatrixTempl<T>::operator*=(T scalar)
{
  for (unsigned int i = 0; i < _n_entries; i++)
    _values[i] *= scalar;
  return *this;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T> &
ColumnMajorMatrixTempl<T>::operator/=(T scalar)
{
  for (unsigned int i = 0; i < _n_entries; i++)
    _values[i] /= scalar;
  return *this;
}
 
template <typename T>
inline ColumnMajorMatrixTempl<T> &
ColumnMajorMatrixTempl<T>::operator+=(T scalar)
{
  for (unsigned int i = 0; i < _n_entries; i++)
    _values[i] += scalar;
  return *this;
}
 
template <typename T>
inline bool
ColumnMajorMatrixTempl<T>::operator==(const ColumnMajorMatrixTempl<T> & rhs) const
{
  if (_n_entries != rhs._n_entries || _n_rows != rhs._n_rows || _n_cols != rhs._n_cols)
    return false;
  return std::equal(_values.begin(), _values.end(), rhs._values.begin());
}
 
template <typename T>
inline bool
ColumnMajorMatrixTempl<T>::operator!=(const ColumnMajorMatrixTempl<T> & rhs) const
{
  return !(*this == rhs);
}
 
typedef ColumnMajorMatrixTempl<Real> ColumnMajorMatrix;
 
namespace MetaPhysicL
{
template <typename T>
struct RawType<ColumnMajorMatrixTempl<T>>
{
  typedef ColumnMajorMatrixTempl<typename RawType<T>::value_type> value_type;
 
  static value_type value(const ColumnMajorMatrixTempl<T> & in)
  {
    value_type ret(in.m(), in.n());
    for (MooseIndex(in.m()) i = 0; i < in.m(); ++i)
      for (MooseIndex(in.n()) j = 0; j < in.n(); ++j)
        ret(i, j) = raw_value(in(i, j));
 
    return ret;
  }
};
}