precondition_selector.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2022 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_precondition_selector_h
#define dealii_precondition_selector_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/subscriptor.h>
 
#include <string>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
template <class number>
class Vector;
template <class number>
class SparseMatrix;
#endif
 
 
template <typename MatrixType = SparseMatrix<double>,
          typename VectorType = ::Vector<double>>
class PreconditionSelector : public Subscriptor
{
public:
  using size_type = typename MatrixType::size_type;
 
  PreconditionSelector(const std::string &                    preconditioning,
                       const typename VectorType::value_type &omega = 1.);
 
  virtual ~PreconditionSelector() override;
 
  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