iterative_inverse.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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__iterative_inverse_h
#define dealii__iterative_inverse_h

#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/lac/solver_selector.h>
#include <deal.II/lac/vector_memory.h>
#include <deal.II/lac/pointer_matrix.h>

DEAL_II_NAMESPACE_OPEN


template <typename VectorType>
class IterativeInverse : public Subscriptor
{
public:
  template <typename MatrixType, typename PreconditionerType>
  void initialize (const MatrixType &, const PreconditionerType &);

  void clear();

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

  template <class OtherVectorType>
  void vmult (OtherVectorType &dst, const OtherVectorType &src) const;

  SolverSelector<VectorType> solver;

private:
  std_cxx11::shared_ptr<PointerMatrixBase<VectorType> > matrix;

  std_cxx11::shared_ptr<PointerMatrixBase<VectorType> > preconditioner;
};


template <typename VectorType>
template <typename MatrixType, typename PreconditionerType>
inline
void
IterativeInverse<VectorType>::initialize(const MatrixType &m, const PreconditionerType &p)
{
  VectorType v;
  matrix = std_cxx11::shared_ptr<PointerMatrixBase<VectorType> > (new_pointer_matrix_base(m, v));
  preconditioner = std_cxx11::shared_ptr<PointerMatrixBase<VectorType> > (new_pointer_matrix_base(p, v));
}


template <typename VectorType>
inline
void
IterativeInverse<VectorType>::clear()
{
  matrix = 0;
  preconditioner = 0;
}


template <typename VectorType>
inline void
IterativeInverse<VectorType>::vmult (VectorType &dst, const VectorType &src) const
{
  Assert(matrix.get() != 0, ExcNotInitialized());
  Assert(preconditioner.get() != 0, ExcNotInitialized());
  dst = 0.;
  solver.solve(*matrix, dst, src, *preconditioner);
}


template <typename VectorType>
template <class OtherVectorType>
inline void
IterativeInverse<VectorType>::vmult (OtherVectorType &dst, const OtherVectorType &src) const
{
  Assert(matrix.get() != 0, ExcNotInitialized());
  Assert(preconditioner.get() != 0, ExcNotInitialized());
  GrowingVectorMemory<VectorType> mem;
  typename VectorMemory<VectorType>::Pointer sol(mem);
  typename VectorMemory<VectorType>::Pointer rhs(mem);
  sol->reinit(dst);
  *rhs = src;
  solver.solve(*matrix, *sol, *rhs, *preconditioner);
  dst = *sol;
}



DEAL_II_NAMESPACE_CLOSE

#endif