matrix_lib.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2002 - 2015 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__matrix_lib_h
#define dealii__matrix_lib_h

#include <deal.II/base/subscriptor.h>
#include <deal.II/lac/vector_memory.h>
#include <deal.II/lac/pointer_matrix.h>
#include <deal.II/lac/solver_richardson.h>

DEAL_II_NAMESPACE_OPEN

template<typename number> class Vector;
template<typename number> class BlockVector;
template<typename number> class SparseMatrix;

template<typename VectorType>
class ProductMatrix : public PointerMatrixBase<VectorType>
{
public:
  ProductMatrix();

  ProductMatrix(VectorMemory<VectorType> &mem);

  template <typename MatrixType1, typename MatrixType2>
  ProductMatrix (const MatrixType1        &m1,
                 const MatrixType2        &m2,
                 VectorMemory<VectorType> &mem);

  ~ProductMatrix();

  template <typename MatrixType1, typename MatrixType2>
  void reinit (const MatrixType1 &m1, const MatrixType2 &m2);

  template <typename MatrixType1, typename MatrixType2>
  void initialize (const MatrixType1 &m1,
                   const MatrixType2 &m2,
                   VectorMemory<VectorType> &mem);

  // Doc in PointerMatrixBase
  void clear();

  virtual void vmult (VectorType       &w,
                      const VectorType &v) const;

  virtual void Tvmult (VectorType       &w,
                       const VectorType &v) const;

  virtual void vmult_add (VectorType       &w,
                          const VectorType &v) const;

  virtual void Tvmult_add (VectorType       &w,
                           const VectorType &v) const;

private:
  PointerMatrixBase<VectorType> *m1;

  PointerMatrixBase<VectorType> *m2;

  SmartPointer<VectorMemory<VectorType>,ProductMatrix<VectorType> > mem;
};


template<typename VectorType>
class ScaledMatrix : public Subscriptor
{
public:
  ScaledMatrix ();
  template <typename MatrixType>
  ScaledMatrix (const MatrixType &M, const double factor);

  ~ScaledMatrix ();
  template <typename MatrixType>
  void initialize (const MatrixType &M, const double factor);

  void clear ();

  void vmult (VectorType &w, const VectorType &v) const;

  void Tvmult (VectorType &w, const VectorType &v) const;

private:
  PointerMatrixBase<VectorType> *m;
  double factor;
};



template<typename number, typename vector_number>
class ProductSparseMatrix : public PointerMatrixBase<Vector<vector_number> >
{
public:
  typedef SparseMatrix<number> MatrixType;

  typedef Vector<vector_number> VectorType;

  ProductSparseMatrix (const MatrixType         &m1,
                       const MatrixType         &m2,
                       VectorMemory<VectorType> &mem);

  ProductSparseMatrix();

  void initialize (const MatrixType         &m1,
                   const MatrixType         &m2,
                   VectorMemory<VectorType> &mem);

  // Doc in PointerMatrixBase
  void clear();

  virtual void vmult (VectorType       &w,
                      const VectorType &v) const;

  virtual void Tvmult (VectorType       &w,
                       const VectorType &v) const;

  virtual void vmult_add (VectorType       &w,
                          const VectorType &v) const;

  virtual void Tvmult_add (VectorType       &w,
                           const VectorType &v) const;

private:
  SmartPointer<const MatrixType,ProductSparseMatrix<number,vector_number> > m1;

  SmartPointer<const MatrixType,ProductSparseMatrix<number,vector_number> > m2;

  SmartPointer<VectorMemory<VectorType>,ProductSparseMatrix<number,vector_number>  > mem;
};


class MeanValueFilter : public Subscriptor
{
public:
  typedef types::global_dof_index size_type;

  MeanValueFilter(const size_type component = numbers::invalid_size_type);

  template <typename number>
  void filter (Vector<number> &v) const;

  template <typename number>
  void filter (BlockVector<number> &v) const;

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

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

  template <typename number>
  void vmult (BlockVector<number>       &dst,
              const BlockVector<number> &src) const;

  template <typename number>
  void vmult_add (BlockVector<number>       &dst,
                  const BlockVector<number> &src) const;


  template <typename VectorType>
  void Tvmult(VectorType &, const VectorType &) const;

  template <typename VectorType>
  void Tvmult_add(VectorType &, const VectorType &) const;

private:
  const size_type component;
};



template<typename VectorType>
class InverseMatrixRichardson : public Subscriptor
{
public:
  InverseMatrixRichardson (SolverControl            &control,
                           VectorMemory<VectorType> &mem);
  ~InverseMatrixRichardson();

  template <typename MatrixType, typename PreconditionerType>
  void initialize (const MatrixType &,
                   const PreconditionerType &);

  SolverControl &control() const;
  void vmult (VectorType &, const VectorType &) const;

  void vmult_add (VectorType &, const VectorType &) const;

  void Tvmult (VectorType &, const VectorType &) const;

  void Tvmult_add (VectorType &, const VectorType &) const;

private:
  VectorMemory<VectorType> &mem;

  mutable SolverRichardson<VectorType> solver;

  PointerMatrixBase<VectorType> *matrix;

  PointerMatrixBase<VectorType> *precondition;
};




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


template<typename VectorType>
inline
ScaledMatrix<VectorType>::ScaledMatrix()
  :
  m(0)
{}



template<typename VectorType>
template<typename MatrixType>
inline
ScaledMatrix<VectorType>::ScaledMatrix(const MatrixType &mat, const double factor)
  :
  m(new_pointer_matrix_base(mat, VectorType())),
  factor(factor)
{}



template<typename VectorType>
template<typename MatrixType>
inline
void
ScaledMatrix<VectorType>::initialize(const MatrixType &mat, const double f)
{
  if (m) delete m;
  m = new_pointer_matrix_base(mat, VectorType());
  factor = f;
}



template<typename VectorType>
inline
void
ScaledMatrix<VectorType>::clear()
{
  if (m) delete m;
  m = 0;
}



template<typename VectorType>
inline
ScaledMatrix<VectorType>::~ScaledMatrix()
{
  clear ();
}


template<typename VectorType>
inline
void
ScaledMatrix<VectorType>::vmult (VectorType &w, const VectorType &v) const
{
  m->vmult(w, v);
  w *= factor;
}


template<typename VectorType>
inline
void
ScaledMatrix<VectorType>::Tvmult (VectorType &w, const VectorType &v) const
{
  m->Tvmult(w, v);
  w *= factor;
}


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

template<typename VectorType>
ProductMatrix<VectorType>::ProductMatrix ()
  : m1(0), m2(0), mem(0)
{}


template<typename VectorType>
ProductMatrix<VectorType>::ProductMatrix (VectorMemory<VectorType> &m)
  : m1(0), m2(0), mem(&m)
{}


template<typename VectorType>
template<typename MatrixType1, typename MatrixType2>
ProductMatrix<VectorType>::ProductMatrix (const MatrixType1        &mat1,
                                          const MatrixType2        &mat2,
                                          VectorMemory<VectorType> &m)
  : mem(&m)
{
  m1 = new PointerMatrix<MatrixType1, VectorType>(&mat1, typeid(*this).name());
  m2 = new PointerMatrix<MatrixType2, VectorType>(&mat2, typeid(*this).name());
}


template<typename VectorType>
template<typename MatrixType1, typename MatrixType2>
void
ProductMatrix<VectorType>::reinit (const MatrixType1 &mat1, const MatrixType2 &mat2)
{
  if (m1) delete m1;
  if (m2) delete m2;
  m1 = new PointerMatrix<MatrixType1, VectorType>(&mat1, typeid(*this).name());
  m2 = new PointerMatrix<MatrixType2, VectorType>(&mat2, typeid(*this).name());
}


template<typename VectorType>
template<typename MatrixType1, typename MatrixType2>
void
ProductMatrix<VectorType>::initialize (const MatrixType1        &mat1,
                                       const MatrixType2        &mat2,
                                       VectorMemory<VectorType> &memory)
{
  mem = &memory;
  if (m1) delete m1;
  if (m2) delete m2;
  m1 = new PointerMatrix<MatrixType1, VectorType>(&mat1, typeid(*this).name());
  m2 = new PointerMatrix<MatrixType2, VectorType>(&mat2, typeid(*this).name());
}


template<typename VectorType>
ProductMatrix<VectorType>::~ProductMatrix ()
{
  if (m1) delete m1;
  if (m2) delete m2;
}


template<typename VectorType>
void
ProductMatrix<VectorType>::clear ()
{
  if (m1) delete m1;
  m1 = 0;
  if (m2) delete m2;
  m2 = 0;
}


template<typename VectorType>
void
ProductMatrix<VectorType>::vmult (VectorType &dst, const VectorType &src) const
{
  Assert (mem != 0, ExcNotInitialized());
  Assert (m1 != 0, ExcNotInitialized());
  Assert (m2 != 0, ExcNotInitialized());

  VectorType *v = mem->alloc();
  v->reinit(dst);
  m2->vmult (*v, src);
  m1->vmult (dst, *v);
  mem->free(v);
}


template<typename VectorType>
void
ProductMatrix<VectorType>::vmult_add (VectorType &dst, const VectorType &src) const
{
  Assert (mem != 0, ExcNotInitialized());
  Assert (m1 != 0, ExcNotInitialized());
  Assert (m2 != 0, ExcNotInitialized());

  VectorType *v = mem->alloc();
  v->reinit(dst);
  m2->vmult (*v, src);
  m1->vmult_add (dst, *v);
  mem->free(v);
}


template<typename VectorType>
void
ProductMatrix<VectorType>::Tvmult (VectorType &dst, const VectorType &src) const
{
  Assert (mem != 0, ExcNotInitialized());
  Assert (m1 != 0, ExcNotInitialized());
  Assert (m2 != 0, ExcNotInitialized());

  VectorType *v = mem->alloc();
  v->reinit(dst);
  m1->Tvmult (*v, src);
  m2->Tvmult (dst, *v);
  mem->free(v);
}


template<typename VectorType>
void
ProductMatrix<VectorType>::Tvmult_add (VectorType &dst, const VectorType &src) const
{
  Assert (mem != 0, ExcNotInitialized());
  Assert (m1 != 0, ExcNotInitialized());
  Assert (m2 != 0, ExcNotInitialized());

  VectorType *v = mem->alloc();
  v->reinit(dst);
  m1->Tvmult (*v, src);
  m2->Tvmult_add (dst, *v);
  mem->free(v);
}



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

template <typename VectorType>
inline void
MeanValueFilter::Tvmult(VectorType &, const VectorType &) const
{
  Assert(false, ExcNotImplemented());
}


template <typename VectorType>
inline void
MeanValueFilter::Tvmult_add(VectorType &, const VectorType &) const
{
  Assert(false, ExcNotImplemented());
}

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

template <typename VectorType>
template <typename MatrixType, typename PreconditionerType>
inline void
InverseMatrixRichardson<VectorType>::initialize (const MatrixType &m,
                                                 const PreconditionerType &p)
{
  if (matrix != 0)
    delete matrix;
  matrix = new PointerMatrix<MatrixType, VectorType>(&m);
  if (precondition != 0)
    delete precondition;
  precondition = new PointerMatrix<PreconditionerType, VectorType>(&p);
}


DEAL_II_NAMESPACE_CLOSE

#endif