dense_submatrix.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
#ifndef LIBMESH_DENSE_SUBMATRIX_H
#define LIBMESH_DENSE_SUBMATRIX_H
 
// Local Includes
#include "libmesh/libmesh_common.h"
#include "libmesh/dense_matrix.h"
#include "libmesh/dense_subvector.h"
#include "libmesh/int_range.h"
 
// C++ includes
 
namespace libMesh
{
 
template<typename T>
class DenseSubMatrix : public DenseMatrixBase<T>
{
public:
 
  DenseSubMatrix(DenseMatrix<T> & new_parent,
                 const unsigned int ioff=0,
                 const unsigned int joff=0,
                 const unsigned int m=0,
                 const unsigned int n=0);
 
  DenseSubMatrix (DenseSubMatrix &&) = default;
  DenseSubMatrix (const DenseSubMatrix &) = default;
  DenseSubMatrix & operator= (const DenseSubMatrix &) = default;
  DenseSubMatrix & operator= (DenseSubMatrix &&) = default;
  virtual ~DenseSubMatrix() = default;
 
  DenseMatrix<T> & parent () { return _parent_matrix; }
 
  virtual void zero() override final;
 
  T operator() (const unsigned int i,
                const unsigned int j) const;
 
  T & operator() (const unsigned int i,
                  const unsigned int j);
 
  virtual T el(const unsigned int i,
               const unsigned int j) const override final
  { return (*this)(i,j); }
 
  virtual T & el(const unsigned int i,
                 const unsigned int j) override final
  { return (*this)(i,j); }
 
  virtual void left_multiply (const DenseMatrixBase<T> & M2) override final;
 
  virtual void right_multiply (const DenseMatrixBase<T> & M3) override final;
 
  void reposition(const unsigned int ioff,
                  const unsigned int joff,
                  const unsigned int new_m,
                  const unsigned int new_n);
 
  unsigned int i_off() const { return _i_off; }
 
  unsigned int j_off() const { return _j_off; }
 
  void condense(const unsigned int i,
                const unsigned int j,
                const T val,
                DenseSubVector<T> & rhs)
  {
    this->parent().condense(this->i_off()+i,
                            this->j_off()+j,
                            val, rhs.parent());
  }
 
private:
 
  DenseMatrix<T> & _parent_matrix;
 
  unsigned int _i_off;
 
  unsigned int _j_off;
};
 
 
// --------------------------------------------------
// Constructor
template<typename T>
inline
DenseSubMatrix<T>::DenseSubMatrix(DenseMatrix<T> & new_parent,
                                  const unsigned int ioff,
                                  const unsigned int joff,
                                  const unsigned int new_m,
                                  const unsigned int new_n) :
  DenseMatrixBase<T>(new_m,new_n),
  _parent_matrix(new_parent)
{
  this->reposition (ioff, joff, new_m, new_n);
}
 
 
template<typename T>
inline
void DenseSubMatrix<T>::reposition(const unsigned int ioff,
                                   const unsigned int joff,
                                   const unsigned int new_m,
                                   const unsigned int new_n)
{
  _i_off = ioff;
  _j_off = joff;
  this->_m = new_m;
  this->_n = new_n;
 
  // Make sure we still fit in the parent matrix.
  libmesh_assert_less_equal ((this->i_off() + this->m()), _parent_matrix.m());
  libmesh_assert_less_equal ((this->j_off() + this->n()), _parent_matrix.n());
}
 
 
 
template<typename T>
inline
void DenseSubMatrix<T>::zero()
{
  for (auto i : make_range(this->m()))
    for (auto j : make_range(this->n()))
      _parent_matrix(i + this->i_off(),
                     j + this->j_off()) = 0.;
}
 
 
 
template<typename T>
inline
T DenseSubMatrix<T>::operator () (const unsigned int i,
                                  const unsigned int j) const
{
  libmesh_assert_less (i, this->m());
  libmesh_assert_less (j, this->n());
  libmesh_assert_less (i + this->i_off(), _parent_matrix.m());
  libmesh_assert_less (j + this->j_off(), _parent_matrix.n());
 
  return _parent_matrix (i + this->i_off(),
                         j + this->j_off());
}
 
 
template<typename T>
inline
T & DenseSubMatrix<T>::operator () (const unsigned int i,
                                    const unsigned int j)
{
  libmesh_assert_less (i, this->m());
  libmesh_assert_less (j, this->n());
  libmesh_assert_less (i + this->i_off(), _parent_matrix.m());
  libmesh_assert_less (j + this->j_off(), _parent_matrix.n());
 
  return _parent_matrix (i + this->i_off(),
                         j + this->j_off());
}
 
 
// A submatrix is finite iff every component is
template <typename T>
bool isfinite (const DenseSubMatrix<T> & var)
{
  using std::isfinite;
  using libMesh::isfinite; // for T==complex
  const auto m = var.m(), n = var.n();
  for (auto i : make_range(m))
    for (auto j : make_range(n))
      if (!isfinite(var(i,j)))
        return false;
  return true;
}
 
 
// A submatrix is infinite iff some component is infinite but no
// component is NaN.
//
// This is arguably inconsistent with our std::complex overload (and
// the C99 Annex G recommendations for _Complex, and C++ std::complex
// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
// this is probably safer for users.
template <typename T>
bool isinf (const DenseSubMatrix<T> & var)
{
  using std::isinf;
  using libMesh::isinf; // for T==complex
  using std::isnan;
  using libMesh::isnan;
  const auto m = var.m(), n = var.n();
  bool has_inf = false;
  for (auto i : make_range(m))
    for (auto j : make_range(n))
      {
        // NaN anywhere makes us NaN, not inf
        if (isnan(var(i,j)))
          return false;
        has_inf = has_inf || isinf(var(i,j));
      }
  return has_inf;
}
 
 
 
// A submatrix is NaN iff some component is NaN
//
// This is arguably inconsistent with our std::complex overload (and
// the C99 Annex G recommendations for _Complex, and C++ std::complex
// arithmetic), which treats mixed (inf,NaN) pairs as infinite, but
// this is probably safer for users.
template <typename T>
bool isnan (const DenseSubMatrix<T> & var)
{
  using std::isnan;
  using libMesh::isnan; // for T==complex
  const auto m = var.m(), n = var.n();
  for (auto i : make_range(m))
    for (auto j : make_range(n))
      if (isnan(var(i,j)))
        return true;
  return false;
}
 
 
} // namespace libMesh
 
 
#endif // LIBMESH_DENSE_SUBMATRIX_H