mg_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2003 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, 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.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_mg_matrix_h
#define dealii_mg_matrix_h

#include <deal.II/base/mg_level_object.h>
#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/multigrid/mg_base.h>

#include <memory>

DEAL_II_NAMESPACE_OPEN


namespace mg
{
  template <typename VectorType = Vector<double> >
  class Matrix
    : public MGMatrixBase<VectorType>
  {
  public:
    Matrix() = default;

    template <typename MatrixType>
    Matrix(const MGLevelObject<MatrixType> &M);

    template <typename MatrixType>
    void
    initialize(const MGLevelObject<MatrixType> &M);

    void reset();

    const LinearOperator<VectorType> &operator[] (unsigned int level) const;

    virtual void vmult (const unsigned int level, VectorType &dst, const VectorType &src) const;
    virtual void vmult_add (const unsigned int level, VectorType &dst, const VectorType &src) const;
    virtual void Tvmult (const unsigned int level, VectorType &dst, const VectorType &src) const;
    virtual void Tvmult_add (const unsigned int level, VectorType &dst, const VectorType &src) const;
    virtual unsigned int get_minlevel() const;
    virtual unsigned int get_maxlevel() const;

    std::size_t memory_consumption () const;
  private:
    MGLevelObject<LinearOperator<VectorType> > matrices;
  };

}


template <typename MatrixType, typename number>
class MGMatrixSelect : public MGMatrixBase<Vector<number> >
{
public:
  MGMatrixSelect (const unsigned int         row    = 0,
                  const unsigned int         col    = 0,
                  MGLevelObject<MatrixType> *matrix = 0);

  void set_matrix (MGLevelObject<MatrixType> *M);

  void select_block (const unsigned int row,
                     const unsigned int col);

  virtual void vmult (const unsigned int level,
                      Vector<number> &dst,
                      const Vector<number> &src) const;

  virtual void vmult_add (const unsigned int level,
                          Vector<number> &dst,
                          const Vector<number> &src) const;

  virtual void Tvmult (const unsigned int level,
                       Vector<number> &dst,
                       const Vector<number> &src) const;

  virtual void Tvmult_add (const unsigned int level,
                           Vector<number> &dst,
                           const Vector<number> &src) const;

private:
  SmartPointer<MGLevelObject<MatrixType>,MGMatrixSelect<MatrixType,number> > matrix;
  unsigned int row;
  unsigned int col;

};

/*----------------------------------------------------------------------*/

namespace mg
{
  template <typename VectorType>
  template <typename MatrixType>
  inline
  void
  Matrix<VectorType>::initialize (const MGLevelObject<MatrixType> &p)
  {
    matrices.resize(p.min_level(), p.max_level());
    for (unsigned int level = p.min_level(); level <= p.max_level(); ++level)
      {
        // Workaround: Unfortunately, not every "p[level]" object has a
        // rich enough interface to populate reinit_(domain|range)_vector.
        // Thus, apply an empty LinearOperator exemplar.
        matrices[level] =
          linear_operator<VectorType>(LinearOperator<VectorType>(), p[level]);
      }
  }



  template <typename VectorType>
  inline
  void
  Matrix<VectorType>::reset ()
  {
    matrices.resize(0,0);
  }



  template <typename VectorType>
  template <typename MatrixType>
  inline
  Matrix<VectorType>::Matrix (const MGLevelObject<MatrixType> &p)
  {
    initialize(p);
  }



  template <typename VectorType>
  inline
  const LinearOperator<VectorType> &
  Matrix<VectorType>::operator[] (unsigned int level) const
  {
    return matrices[level];
  }



  template <typename VectorType>
  void
  Matrix<VectorType>::vmult (const unsigned int level,
                             VectorType         &dst,
                             const VectorType   &src) const
  {
    matrices[level].vmult(dst, src);
  }



  template <typename VectorType>
  void
  Matrix<VectorType>::vmult_add (const unsigned int level,
                                 VectorType         &dst,
                                 const VectorType   &src) const
  {
    matrices[level].vmult_add(dst, src);
  }



  template <typename VectorType>
  void
  Matrix<VectorType>::Tvmult (const unsigned int level,
                              VectorType         &dst,
                              const VectorType   &src) const
  {
    matrices[level].Tvmult(dst, src);
  }



  template <typename VectorType>
  void
  Matrix<VectorType>::Tvmult_add (const unsigned int level,
                                  VectorType         &dst,
                                  const VectorType   &src) const
  {
    matrices[level].Tvmult_add(dst, src);
  }



  template <typename VectorType>
  unsigned int
  Matrix<VectorType>::get_minlevel() const
  {
    return matrices.min_level();
  }



  template <typename VectorType>
  unsigned int
  Matrix<VectorType>::get_maxlevel() const
  {
    return matrices.max_level();
  }



  template <typename VectorType>
  inline
  std::size_t
  Matrix<VectorType>::memory_consumption () const
  {
    return sizeof(*this) + matrices->memory_consumption();
  }
}


/*----------------------------------------------------------------------*/

template <typename MatrixType, typename number>
MGMatrixSelect<MatrixType, number>::
MGMatrixSelect (const unsigned int         row,
                const unsigned int         col,
                MGLevelObject<MatrixType> *p)
  :
  matrix (p, typeid(*this).name()),
  row(row),
  col(col)
{}



template <typename MatrixType, typename number>
void
MGMatrixSelect<MatrixType, number>::set_matrix (MGLevelObject<MatrixType> *p)
{
  matrix = p;
}



template <typename MatrixType, typename number>
void
MGMatrixSelect<MatrixType, number>::
select_block (const unsigned int brow,
              const unsigned int bcol)
{
  row = brow;
  col = bcol;
}



template <typename MatrixType, typename number>
void
MGMatrixSelect<MatrixType, number>::
vmult  (const unsigned int    level,
        Vector<number>       &dst,
        const Vector<number> &src) const
{
  Assert(matrix != 0, ExcNotInitialized());

  const MGLevelObject<MatrixType> &m = *matrix;
  m[level].block(row, col).vmult(dst, src);
}



template <typename MatrixType, typename number>
void
MGMatrixSelect<MatrixType, number>::
vmult_add  (const unsigned int    level,
            Vector<number>       &dst,
            const Vector<number> &src) const
{
  Assert(matrix != 0, ExcNotInitialized());

  const MGLevelObject<MatrixType> &m = *matrix;
  m[level].block(row, col).vmult_add(dst, src);
}



template <typename MatrixType, typename number>
void
MGMatrixSelect<MatrixType, number>::
Tvmult  (const unsigned int    level,
         Vector<number>       &dst,
         const Vector<number> &src) const
{
  Assert(matrix != 0, ExcNotInitialized());

  const MGLevelObject<MatrixType> &m = *matrix;
  m[level].block(row, col).Tvmult(dst, src);
}



template <typename MatrixType, typename number>
void
MGMatrixSelect<MatrixType, number>::
Tvmult_add  (const unsigned int    level,
             Vector<number>       &dst,
             const Vector<number> &src) const
{
  Assert(matrix != 0, ExcNotInitialized());

  const MGLevelObject<MatrixType> &m = *matrix;
  m[level].block(row, col).Tvmult_add(dst, src);
}

DEAL_II_NAMESPACE_CLOSE

#endif