precondition_selector.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_precondition_selector_h
#define dealii_precondition_selector_h


#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
#include <string>

DEAL_II_NAMESPACE_OPEN

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


template <typename MatrixType = SparseMatrix<double>,
          typename VectorType = ::Vector<double> >
class PreconditionSelector : public Subscriptor
{
public:
  typedef typename MatrixType::size_type size_type;

  PreconditionSelector (const std::string                     &preconditioning,
                        const typename VectorType::value_type &omega=1.);

  virtual ~PreconditionSelector();

  void use_matrix(const MatrixType &M);

  size_type m () const;

  size_type n () const;

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

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

  static std::string get_precondition_names();

  DeclException0 (ExcNoMatrixGivenToUse);

protected:

  std::string preconditioning;

private:
  SmartPointer<const MatrixType,PreconditionSelector<MatrixType,VectorType> > A;

  const typename VectorType::value_type omega;
};

/* --------------------- Inline and template functions ------------------- */


template <typename MatrixType, typename VectorType>
PreconditionSelector<MatrixType,VectorType>
::PreconditionSelector(const std::string                     &preconditioning,
                       const typename VectorType::value_type &omega) :
  preconditioning(preconditioning),
  omega(omega)  {}


template <typename MatrixType, typename VectorType>
PreconditionSelector<MatrixType,VectorType>::~PreconditionSelector()
{
  // release the matrix A
  A=nullptr;
}


template <typename MatrixType, typename VectorType>
void PreconditionSelector<MatrixType,VectorType>::use_matrix(const MatrixType &M)
{
  A=&M;
}


template <typename MatrixType, typename VectorType>
inline typename PreconditionSelector<MatrixType,VectorType>::size_type
PreconditionSelector<MatrixType,VectorType>::m () const
{
  Assert(A!=nullptr, ExcNoMatrixGivenToUse());
  return A->m();
}


template <typename MatrixType, typename VectorType>
inline typename PreconditionSelector<MatrixType,VectorType>::size_type
PreconditionSelector<MatrixType,VectorType>::n () const
{
  Assert(A!=nullptr, ExcNoMatrixGivenToUse());
  return A->n();
}



template <typename MatrixType, typename VectorType>
void PreconditionSelector<MatrixType,VectorType>::vmult (VectorType &dst,
                                                         const VectorType &src) const
{
  if (preconditioning=="none")
    {
      dst=src;
    }
  else
    {
      Assert(A!=nullptr, ExcNoMatrixGivenToUse());

      if (preconditioning=="jacobi")
        {
          A->precondition_Jacobi(dst,src,omega);
        }
      else if (preconditioning=="sor")
        {
          A->precondition_SOR(dst,src,omega);
        }
      else if (preconditioning=="ssor")
        {
          A->precondition_SSOR(dst,src,omega);
        }
      else
        Assert(false,ExcNotImplemented());
    }
}


template <typename MatrixType, typename VectorType>
void PreconditionSelector<MatrixType,VectorType>::Tvmult (VectorType &dst,
                                                          const VectorType &src) const
{
  if (preconditioning=="none")
    {
      dst=src;
    }
  else
    {
      Assert(A!=nullptr, ExcNoMatrixGivenToUse());

      if (preconditioning=="jacobi")
        {
          A->precondition_Jacobi(dst,src,omega); // Symmetric operation
        }
      else if (preconditioning=="sor")
        {
          A->precondition_TSOR(dst,src,omega);
        }
      else if (preconditioning=="ssor")
        {
          A->precondition_SSOR(dst,src,omega); // Symmetric operation
        }
      else
        Assert(false,ExcNotImplemented());
    }
}


template <typename MatrixType, typename VectorType>
std::string PreconditionSelector<MatrixType,VectorType>::get_precondition_names()
{
  return "none|jacobi|sor|ssor";
}


DEAL_II_NAMESPACE_CLOSE

#endif