pointer_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2002 - 2018 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_pointer_matrix_h
#define dealii_pointer_matrix_h

#include <deal.II/base/smartpointer.h>
#include <deal.II/base/subscriptor.h>

#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_memory.h>

DEAL_II_NAMESPACE_OPEN

template <typename VectorType>
class VectorMemory;

class IdentityMatrix;
template <typename number>
class FullMatrix;
template <typename number>
class LAPACKFullMatrix;
template <typename number>
class SparseMatrix;
template <typename number>
class BlockSparseMatrix;
template <typename number>
class SparseMatrixEZ;
template <typename number>
class BlockSparseMatrixEZ;
template <typename number>
class TridiagonalMatrix;
template <typename number, typename BlockVectorType>
class BlockMatrixArray;

template <typename VectorType>
class DEAL_II_DEPRECATED PointerMatrixBase : public Subscriptor
{
public:
  using value_type = typename VectorType::value_type;

  virtual ~PointerMatrixBase() override = default;

  virtual void
  clear() = 0;

  virtual void
  vmult(VectorType &dst, const VectorType &src) const = 0;

  virtual void
  Tvmult(VectorType &dst, const VectorType &src) const = 0;

  virtual void
  vmult_add(VectorType &dst, const VectorType &src) const = 0;

  virtual void
  Tvmult_add(VectorType &dst, const VectorType &src) const = 0;
};


template <typename MatrixType, typename VectorType>
class DEAL_II_DEPRECATED PointerMatrix : public PointerMatrixBase<VectorType>
{
public:
  PointerMatrix(const MatrixType *M = nullptr);

  PointerMatrix(const char *name);

  PointerMatrix(const MatrixType *M, const char *name);

  // Use doc from base class
  virtual void
  clear() override;

  bool
  empty() const;

  const PointerMatrix &
  operator=(const MatrixType *M);

  virtual void
  vmult(VectorType &dst, const VectorType &src) const override;

  virtual void
  Tvmult(VectorType &dst, const VectorType &src) const override;

  virtual void
  vmult_add(VectorType &dst, const VectorType &src) const override;

  virtual void
  Tvmult_add(VectorType &dst, const VectorType &src) const override;

private:
  SmartPointer<const MatrixType, PointerMatrix<MatrixType, VectorType>> m;
};


template <typename MatrixType, typename VectorType>
class DEAL_II_DEPRECATED PointerMatrixAux : public PointerMatrixBase<VectorType>
{
public:
  PointerMatrixAux(VectorMemory<VectorType> *mem = 0, const MatrixType *M = 0);

  PointerMatrixAux(VectorMemory<VectorType> *mem, const char *name);

  PointerMatrixAux(VectorMemory<VectorType> *mem,
                   const MatrixType *        M,
                   const char *              name);

  // Use doc from base class
  virtual void
  clear() override;

  bool
  empty() const;

  void
  set_memory(VectorMemory<VectorType> *mem);

  const PointerMatrixAux &
  operator=(const MatrixType *M);

  virtual void
  vmult(VectorType &dst, const VectorType &src) const override;

  virtual void
  Tvmult(VectorType &dst, const VectorType &src) const override;

  virtual void
  vmult_add(VectorType &dst, const VectorType &src) const override;

  virtual void
  Tvmult_add(VectorType &dst, const VectorType &src) const override;

private:
  mutable GrowingVectorMemory<VectorType> my_memory;

  mutable SmartPointer<VectorMemory<VectorType>,
                       PointerMatrixAux<MatrixType, VectorType>>
    mem;

  SmartPointer<const MatrixType, PointerMatrixAux<MatrixType, VectorType>> m;
};



template <typename number>
class DEAL_II_DEPRECATED PointerMatrixVector
  : public PointerMatrixBase<Vector<number>>
{
public:
  PointerMatrixVector(const Vector<number> *M = 0);

  PointerMatrixVector(const char *name);

  PointerMatrixVector(const Vector<number> *M, const char *name);

  // Use doc from base class
  virtual void
  clear();

  bool
  empty() const;

  const PointerMatrixVector &
  operator=(const Vector<number> *M);

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

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

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

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

private:
  SmartPointer<const Vector<number>, PointerMatrixVector<number>> m;
};



template <typename VectorType, typename MatrixType>
inline PointerMatrixBase<VectorType> *
new_pointer_matrix_base(MatrixType &matrix,
                        const VectorType &,
                        const char *name = "PointerMatrixAux")
{
  return new PointerMatrixAux<MatrixType, VectorType>(nullptr, &matrix, name);
}

template <typename numberv>
PointerMatrixBase<Vector<numberv>> *
new_pointer_matrix_base(const IdentityMatrix &matrix,
                        const Vector<numberv> &,
                        const char *name = "PointerMatrix")
{
  return new PointerMatrix<IdentityMatrix, Vector<numberv>>(&matrix, name);
}


template <typename numberv, typename numberm>
PointerMatrixBase<Vector<numberv>> *
new_pointer_matrix_base(const FullMatrix<numberm> &matrix,
                        const Vector<numberv> &,
                        const char *name = "PointerMatrix")
{
  return new PointerMatrix<FullMatrix<numberm>, Vector<numberv>>(&matrix, name);
}


template <typename numberv, typename numberm>
PointerMatrixBase<Vector<numberv>> *
new_pointer_matrix_base(const LAPACKFullMatrix<numberm> &matrix,
                        const Vector<numberv> &,
                        const char *name = "PointerMatrix")
{
  return new PointerMatrix<LAPACKFullMatrix<numberm>, Vector<numberv>>(&matrix,
                                                                       name);
}


template <typename numberv, typename numberm>
PointerMatrixBase<Vector<numberv>> *
new_pointer_matrix_base(const SparseMatrix<numberm> &matrix,
                        const Vector<numberv> &,
                        const char *name = "PointerMatrix")
{
  return new PointerMatrix<SparseMatrix<numberm>, Vector<numberv>>(&matrix,
                                                                   name);
}


template <typename VectorType, typename numberm>
PointerMatrixBase<VectorType> *
new_pointer_matrix_base(const BlockSparseMatrix<numberm> &matrix,
                        const VectorType &,
                        const char *name = "PointerMatrix")
{
  return new PointerMatrix<BlockSparseMatrix<numberm>, VectorType>(&matrix,
                                                                   name);
}


template <typename numberv, typename numberm>
PointerMatrixBase<Vector<numberv>> *
new_pointer_matrix_base(const SparseMatrixEZ<numberm> &matrix,
                        const Vector<numberv> &,
                        const char *name = "PointerMatrix")
{
  return new PointerMatrix<SparseMatrixEZ<numberm>, Vector<numberv>>(&matrix,
                                                                     name);
}


template <typename VectorType, typename numberm>
PointerMatrixBase<VectorType> *
new_pointer_matrix_base(const BlockSparseMatrixEZ<numberm> &matrix,
                        const VectorType &,
                        const char *name = "PointerMatrix")
{
  return new PointerMatrix<BlockSparseMatrixEZ<numberm>, VectorType>(&matrix,
                                                                     name);
}


template <typename numberv, typename numberm, typename BLOCK_VectorType>
PointerMatrixBase<BLOCK_VectorType> *
new_pointer_matrix_base(
  const BlockMatrixArray<numberm, BLOCK_VectorType> &matrix,
  const BLOCK_VectorType &,
  const char *name = "PointerMatrix")
{
  return new PointerMatrix<BlockMatrixArray<numberm, BLOCK_VectorType>,
                           BlockVector<numberv>>(&matrix, name);
}


template <typename numberv, typename numberm>
PointerMatrixBase<Vector<numberv>> *
new_pointer_matrix_base(const TridiagonalMatrix<numberm> &matrix,
                        const Vector<numberv> &,
                        const char *name = "PointerMatrix")
{
  return new PointerMatrix<TridiagonalMatrix<numberm>, Vector<numberv>>(&matrix,
                                                                        name);
}



//---------------------------------------------------------------------------


template <typename MatrixType, typename VectorType>
PointerMatrix<MatrixType, VectorType>::PointerMatrix(const MatrixType *M)
  : m(M, typeid(*this).name())
{}


template <typename MatrixType, typename VectorType>
PointerMatrix<MatrixType, VectorType>::PointerMatrix(const char *name)
  : m(nullptr, name)
{}


template <typename MatrixType, typename VectorType>
PointerMatrix<MatrixType, VectorType>::PointerMatrix(const MatrixType *M,
                                                     const char *      name)
  : m(M, name)
{}


template <typename MatrixType, typename VectorType>
inline void
PointerMatrix<MatrixType, VectorType>::clear()
{
  m = nullptr;
}


template <typename MatrixType, typename VectorType>
inline const PointerMatrix<MatrixType, VectorType> &
PointerMatrix<MatrixType, VectorType>::operator=(const MatrixType *M)
{
  m = M;
  return *this;
}


template <typename MatrixType, typename VectorType>
inline bool
PointerMatrix<MatrixType, VectorType>::empty() const
{
  if (m == nullptr)
    return true;
  return m->empty();
}

template <typename MatrixType, typename VectorType>
inline void
PointerMatrix<MatrixType, VectorType>::vmult(VectorType &      dst,
                                             const VectorType &src) const
{
  Assert(m != nullptr, ExcNotInitialized());
  m->vmult(dst, src);
}


template <typename MatrixType, typename VectorType>
inline void
PointerMatrix<MatrixType, VectorType>::Tvmult(VectorType &      dst,
                                              const VectorType &src) const
{
  Assert(m != nullptr, ExcNotInitialized());
  m->Tvmult(dst, src);
}


template <typename MatrixType, typename VectorType>
inline void
PointerMatrix<MatrixType, VectorType>::vmult_add(VectorType &      dst,
                                                 const VectorType &src) const
{
  Assert(m != nullptr, ExcNotInitialized());
  m->vmult_add(dst, src);
}


template <typename MatrixType, typename VectorType>
inline void
PointerMatrix<MatrixType, VectorType>::Tvmult_add(VectorType &      dst,
                                                  const VectorType &src) const
{
  Assert(m != nullptr, ExcNotInitialized());
  m->Tvmult_add(dst, src);
}



//----------------------------------------------------------------------//


template <typename MatrixType, typename VectorType>
PointerMatrixAux<MatrixType, VectorType>::PointerMatrixAux(
  VectorMemory<VectorType> *mem,
  const MatrixType *        M)
  : mem(mem, typeid(*this).name())
  , m(M, typeid(*this).name())
{
  if (mem == 0)
    mem = &my_memory;
}


template <typename MatrixType, typename VectorType>
PointerMatrixAux<MatrixType, VectorType>::PointerMatrixAux(
  VectorMemory<VectorType> *mem,
  const char *              name)
  : mem(mem, name)
  , m(0, name)
{
  if (mem == 0)
    mem = &my_memory;
}


template <typename MatrixType, typename VectorType>
PointerMatrixAux<MatrixType, VectorType>::PointerMatrixAux(
  VectorMemory<VectorType> *mem,
  const MatrixType *        M,
  const char *              name)
  : mem(mem, name)
  , m(M, name)
{
  if (mem == nullptr)
    mem = &my_memory;
}


template <typename MatrixType, typename VectorType>
inline void
PointerMatrixAux<MatrixType, VectorType>::clear()
{
  m = nullptr;
}


template <typename MatrixType, typename VectorType>
inline const PointerMatrixAux<MatrixType, VectorType> &
PointerMatrixAux<MatrixType, VectorType>::operator=(const MatrixType *M)
{
  m = M;
  return *this;
}


template <typename MatrixType, typename VectorType>
inline void
PointerMatrixAux<MatrixType, VectorType>::set_memory(
  VectorMemory<VectorType> *M)
{
  mem = M;
  if (mem == 0)
    mem = &my_memory;
}


template <typename MatrixType, typename VectorType>
inline bool
PointerMatrixAux<MatrixType, VectorType>::empty() const
{
  if (m == 0)
    return true;
  return m->empty();
}

template <typename MatrixType, typename VectorType>
inline void
PointerMatrixAux<MatrixType, VectorType>::vmult(VectorType &      dst,
                                                const VectorType &src) const
{
  if (mem == nullptr)
    mem = &my_memory;
  Assert(mem != nullptr, ExcNotInitialized());
  Assert(m != nullptr, ExcNotInitialized());
  m->vmult(dst, src);
}


template <typename MatrixType, typename VectorType>
inline void
PointerMatrixAux<MatrixType, VectorType>::Tvmult(VectorType &      dst,
                                                 const VectorType &src) const
{
  if (mem == nullptr)
    mem = &my_memory;
  Assert(mem != nullptr, ExcNotInitialized());
  Assert(m != nullptr, ExcNotInitialized());
  m->Tvmult(dst, src);
}


template <typename MatrixType, typename VectorType>
inline void
PointerMatrixAux<MatrixType, VectorType>::vmult_add(VectorType &      dst,
                                                    const VectorType &src) const
{
  if (mem == nullptr)
    mem = &my_memory;
  Assert(mem != nullptr, ExcNotInitialized());
  Assert(m != nullptr, ExcNotInitialized());
  VectorType *v = mem->alloc();
  v->reinit(dst);
  m->vmult(*v, src);
  dst += *v;
  mem->free(v);
}


template <typename MatrixType, typename VectorType>
inline void
PointerMatrixAux<MatrixType, VectorType>::Tvmult_add(
  VectorType &      dst,
  const VectorType &src) const
{
  if (mem == nullptr)
    mem = &my_memory;
  Assert(mem != nullptr, ExcNotInitialized());
  Assert(m != nullptr, ExcNotInitialized());
  VectorType *v = mem->alloc();
  v->reinit(dst);
  m->Tvmult(*v, src);
  dst += *v;
  mem->free(v);
}


//----------------------------------------------------------------------//


template <typename number>
PointerMatrixVector<number>::PointerMatrixVector(const Vector<number> *M)
  : m(M, typeid(*this).name())
{}


template <typename number>
PointerMatrixVector<number>::PointerMatrixVector(const char *name)
  : m(0, name)
{}


template <typename number>
PointerMatrixVector<number>::PointerMatrixVector(const Vector<number> *M,
                                                 const char *          name)
  : m(M, name)
{}


template <typename number>
inline void
PointerMatrixVector<number>::clear()
{
  m = nullptr;
}


template <typename number>
inline const PointerMatrixVector<number> &
PointerMatrixVector<number>::operator=(const Vector<number> *M)
{
  m = M;
  return *this;
}


template <typename number>
inline bool
PointerMatrixVector<number>::empty() const
{
  if (m == 0)
    return true;
  return m->empty();
}

template <typename number>
inline void
PointerMatrixVector<number>::vmult(Vector<number> &      dst,
                                   const Vector<number> &src) const
{
  Assert(m != nullptr, ExcNotInitialized());
  Assert(dst.size() == 1, ExcDimensionMismatch(dst.size(), 1));

  dst(0) = *m * src;
}


template <typename number>
inline void
PointerMatrixVector<number>::Tvmult(Vector<number> &      dst,
                                    const Vector<number> &src) const
{
  Assert(m != nullptr, ExcNotInitialized());
  Assert(src.size() == 1, ExcDimensionMismatch(src.size(), 1));

  dst.equ(src(0), *m);
}


template <typename number>
inline void
PointerMatrixVector<number>::vmult_add(Vector<number> &      dst,
                                       const Vector<number> &src) const
{
  Assert(m != nullptr, ExcNotInitialized());
  Assert(dst.size() == 1, ExcDimensionMismatch(dst.size(), 1));

  dst(0) += *m * src;
}


template <typename number>
inline void
PointerMatrixVector<number>::Tvmult_add(Vector<number> &      dst,
                                        const Vector<number> &src) const
{
  Assert(m != nullptr, ExcNotInitialized());
  Assert(src.size() == 1, ExcDimensionMismatch(src.size(), 1));

  dst.add(src(0), *m);
}


DEAL_II_NAMESPACE_CLOSE

#endif