solver_selector.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2019 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_solver_selector_h
#define dealii_solver_selector_h


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

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

#include <deal.II/lac/precondition.h>
#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/solver_richardson.h>
#include <deal.II/lac/vector.h>

DEAL_II_NAMESPACE_OPEN



template <typename VectorType = Vector<double>>
class SolverSelector : public Subscriptor
{
public:
  using vector_type = VectorType;

  SolverSelector() = default;

  SolverSelector(const std::string &name, SolverControl &control);

  virtual ~SolverSelector() override = default;

  template <class Matrix, class Preconditioner>
  void
  solve(const Matrix &        A,
        VectorType &          x,
        const VectorType &    b,
        const Preconditioner &precond) const;

  void
  select(const std::string &name);

  void
  set_control(SolverControl &ctrl);

  void
  set_data(const typename SolverRichardson<VectorType>::AdditionalData &data);

  void
  set_data(const typename SolverCG<VectorType>::AdditionalData &data);

  void
  set_data(const typename SolverMinRes<VectorType>::AdditionalData &data);

  void
  set_data(const typename SolverBicgstab<VectorType>::AdditionalData &data);

  void
  set_data(const typename SolverGMRES<VectorType>::AdditionalData &data);

  void
  set_data(const typename SolverFGMRES<VectorType>::AdditionalData &data);

  static std::string
  get_solver_names();

  DeclException1(ExcSolverDoesNotExist,
                 std::string,
                 << "Solver " << arg1 << " does not exist. Use one of "
                 << std::endl
                 << get_solver_names());



protected:
  SmartPointer<SolverControl, SolverSelector<VectorType>> control;

  std::string solver_name;

private:
  typename SolverRichardson<VectorType>::AdditionalData richardson_data;

  typename SolverCG<VectorType>::AdditionalData cg_data;

  typename SolverMinRes<VectorType>::AdditionalData minres_data;

  typename SolverBicgstab<VectorType>::AdditionalData bicgstab_data;

  typename SolverGMRES<VectorType>::AdditionalData gmres_data;

  typename SolverFGMRES<VectorType>::AdditionalData fgmres_data;
};

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


template <typename VectorType>
SolverSelector<VectorType>::SolverSelector(const std::string &name,
                                           SolverControl &    solver_control)
  : solver_name(name)
  , control(&solver_control)
{}



template <typename VectorType>
void
SolverSelector<VectorType>::select(const std::string &name)
{
  solver_name = name;
}



template <typename VectorType>
template <class Matrix, class Preconditioner>
void
SolverSelector<VectorType>::solve(const Matrix &        A,
                                  VectorType &          x,
                                  const VectorType &    b,
                                  const Preconditioner &precond) const
{
  if (solver_name == "richardson")
    {
      SolverRichardson<VectorType> solver(*control, richardson_data);
      solver.solve(A, x, b, precond);
    }
  else if (solver_name == "cg")
    {
      SolverCG<VectorType> solver(*control, cg_data);
      solver.solve(A, x, b, precond);
    }
  else if (solver_name == "minres")
    {
      SolverMinRes<VectorType> solver(*control, minres_data);
      solver.solve(A, x, b, precond);
    }
  else if (solver_name == "bicgstab")
    {
      SolverBicgstab<VectorType> solver(*control, bicgstab_data);
      solver.solve(A, x, b, precond);
    }
  else if (solver_name == "gmres")
    {
      SolverGMRES<VectorType> solver(*control, gmres_data);
      solver.solve(A, x, b, precond);
    }
  else if (solver_name == "fgmres")
    {
      SolverFGMRES<VectorType> solver(*control, fgmres_data);
      solver.solve(A, x, b, precond);
    }
  else
    Assert(false, ExcSolverDoesNotExist(solver_name));
}



template <typename VectorType>
void
SolverSelector<VectorType>::set_control(SolverControl &ctrl)
{
  control = &ctrl;
}



template <typename VectorType>
std::string
SolverSelector<VectorType>::get_solver_names()
{
  return "richardson|cg|bicgstab|gmres|fgmres|minres";
}



template <typename VectorType>
void
SolverSelector<VectorType>::set_data(
  const typename SolverGMRES<VectorType>::AdditionalData &data)
{
  gmres_data = data;
}



template <typename VectorType>
void
SolverSelector<VectorType>::set_data(
  const typename SolverFGMRES<VectorType>::AdditionalData &data)
{
  fgmres_data = data;
}



template <typename VectorType>
void
SolverSelector<VectorType>::set_data(
  const typename SolverRichardson<VectorType>::AdditionalData &data)
{
  richardson_data = data;
}



template <typename VectorType>
void
SolverSelector<VectorType>::set_data(
  const typename SolverCG<VectorType>::AdditionalData &data)
{
  cg_data = data;
}



template <typename VectorType>
void
SolverSelector<VectorType>::set_data(
  const typename SolverMinRes<VectorType>::AdditionalData &data)
{
  minres_data = data;
}



template <typename VectorType>
void
SolverSelector<VectorType>::set_data(
  const typename SolverBicgstab<VectorType>::AdditionalData &data)
{
  bicgstab_data = data;
}

DEAL_II_NAMESPACE_CLOSE

#endif