Commit 5b2dce34 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

startet to adopt istl solvers to new interface

parent cba97849
......@@ -55,6 +55,10 @@ namespace AMDiS
virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
SolverInfo& solverInfo) = 0;
};
template <class Matrix, class Vector>
struct SolverCreator;
} // end namespace AMDiS
#pragma once
#include <dune/amdis/linear_algebra/istl/ISTL_Solver.hpp>
#include <dune/amdis/linear_algebra/istl/Preconditioner.hpp>
namespace AMDiS
{
template <class Matrix, class VectorX, class VectorB>
class ISTLRunner
: public RunnerInterface<Matrix, VectorX, VectorB>
{
using LinOperator = Dune::MatrixAdapter<Matrix, VectorX, VectorB>;
using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
using PreconRunner = ISTLPreconditioner<Matrix, VectorX, VectorB>;
using PreconBase = PreconditionerInterface<Matrix, VectorX, VectorB>;
public:
ISTLRunner(std::string solverName, std::string prefix)
: solverName(solverName)
, prefix(prefix)
, solverCreator(prefix)
{
// get creator string for preconditioner
std::string preconName = "no";
Parameters::get(prefix + "->left precon", preconName);
if (preconName.empty() || preconName == "no")
Parameters::get(prefix + "->right precon", preconName);
if (preconName.empty() || preconName == "no")
Parameters::get(prefix + "->precon->name", preconName);
precon = std::make_shared<PreconRunner>(preconName, prefix + "->precon");
}
/// Implementes \ref RunnerInterface::init()
virtual void init(Matrix const& A) override
{
precon->init(A);
Matrix& _A = const_cast<Matrix&>(A); // Unschoen!!!
linOperator = std::make_shared<LinOperator>(_A);
solver = solverCreator.create(solverName, *linOperator, *precon);
}
/// Implementes \ref RunnerInterface::exit()
virtual void exit() override
{
solver.reset();
linOperator.reset();
precon->exit();
precon.reset();
}
/// Implementes \ref RunnerInterface::solve()
virtual int solve(Matrix const& A, VectorX& x, VectorB const& b,
SolverInfo& solverInfo) override
{
// storing some statistics
Dune::InverseOperatorResult statistics;
// solve the linear system
VectorB _b = b;
solver->apply(x, _b, statistics);
solverInfo.setRelResidual(statistics.reduction);
}
/// Implementes \ref RunnerInterface::getLeftPrecon()
virtual std::shared_ptr<PreconBase> getLeftPrecon() override
{
return preconAdapter;
}
private:
std::string solverName;
std::string prefix;
ISTLSolverCreator<Matrix, VectorX, VectorB> solverCreator;
std::shared_ptr<PreconRunner> precon;
std::shared_ptr<LinOperator> linOperator;
std::shared_ptr<BaseSolver> solver;
};
} // end namespace AMDiS
#pragma once
#include <dune/istl/preconditioners.hh>
namespace AMDiS
{
// creators for preconditioners, depending on matrix-type
// ---------------------------------------------------------------------------
template <class Matrix, class VectorX, class VectorB = VectorX>
class ISTLPreconCreator;
template <class... MatrixRows, class VectorX, class VectorB>
class ISTLPreconCreator<Dune::MultiTypeBlockMatrix<MatrixRows...>, VectorX, VectorB>
{
using BasePreconditioner = Dune::Preconditioner<VectorX, VectorB>;
public:
ISTLPreconCreator(std::string prefix)
: prefix(prefix)
{}
template <class Matrix>
std::unique_ptr<BasePreconditioner> create(std::string /*name*/, Matrix const& A)
{
double w = 1.0;
Parameters::get(prefix + "->relaxation", w);
AMDIS_WARNING("Only Richardson preconditioner available for MultiTypeBlockMatrix!");
return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w);
}
private:
std::string prefix;
};
template <class Block, class Alloc, class VectorX, class VectorB>
class ISTLPreconCreator<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>
{
using BasePreconditioner = Dune::Preconditioner<VectorX, VectorB>;
public:
ISTLPreconCreator(std::string prefix)
: prefix(prefix)
{}
template <class Matrix>
std::unique_ptr<BasePreconditioner> create(std::string name, Matrix const& A)
{
double w = 1.0;
Parameters::get(prefix + "->relaxation", w);
if (name == "diag" || name == "jacobi")
return std::make_unique<Dune::SeqJac<Matrix, VectorX, VectorB>>(A, 1, w);
else if (name == "gs" || name == "gauss_seidel")
return std::make_unique<Dune::SeqGS<Matrix, VectorX, VectorB>>(A, 1, w);
else if (name == "sor")
return std::make_unique<Dune::SeqSOR<Matrix, VectorX, VectorB>>(A, 1, w);
else if (name == "ssor")
return std::make_unique<Dune::SeqSSOR<Matrix, VectorX, VectorB>>(A, 1, w);
else if (name == "ilu" || name == "ilu0")
return std::make_unique<Dune::SeqILU0<Matrix, VectorX, VectorB>>(A, w);
else
return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w);
}
private:
std::string prefix;
};
} // end namespace AMDiS
#pragma once
#include <dune/istl/solvers.hh>
#include <dune/istl/umfpack.hh>
namespace AMDiS
{
// creators for solvers, depending on matrix-type
// ---------------------------------------------------------------------------
template <class Matrix, class VectorX, class VectorB = VectorX>
class ISTLSolverCreator;
template <class... MatrixRows, class VectorX, class VectorB>
class ISTLSolverCreator<Dune::MultiTypeBlockMatrix<MatrixRows...>, VectorX, VectorB>
{
using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
public:
ISTLSolverCreator(std::string prefix)
: prefix(prefix)
, info(2)
{
Parameters::get(prefix + "->info", info);
}
template <class Matrix, class Precon>
std::unique_ptr<BaseSolver> create(std::string name, Matrix& A, Precon& P)
{
size_t maxIter = 500;
Parameters::get(prefix + "->max iteration", maxIter);
double rTol = 1.e-6;
Parameters::get(prefix + "->reduction", rTol);
if (name == "cg")
return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "bicgstab")
return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "minres")
return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "gmres")
{
int restart = 30;
Parameters::get(prefix + "->restart", restart);
return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info);
}
else
AMDIS_ERROR_EXIT("Unknown solver name!");
}
private:
std::string prefix;
int info;
};
template <class Block, class Alloc, class VectorX, class VectorB>
class ISTLSolverCreator<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>
{
using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
using Matrix = Dune::BCRSMatrix<Block,Alloc>;
public:
ISTLSolverCreator(std::string prefix)
: prefix(prefix)
, info(2)
{
Parameters::get(prefix + "->info", info);
}
template <class LinOperator, class Precon>
std::unique_ptr<BaseSolver> create(std::string name, LinOperator& A, Precon& P)
{
size_t maxIter = 500;
Parameters::get(prefix + "->max iteration", maxIter);
double rTol = 1.e-6;
Parameters::get(prefix + "->relative tolerance", rTol);
if (name == "cg")
return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "bicgstab")
return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "minres")
return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "gmres")
{
int restart = 30;
Parameters::get(prefix + "->restart", restart);
return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info);
}
#if HAVE_SUITESPARSE_UMFPACK
else if (name == "umfpack" || name == "direct")
return std::make_unique<Dune::UMFPack<Matrix>>(A, info);
#endif
else
AMDIS_ERROR_EXIT("Unknown solver name!");
}
private:
std::string prefix;
int info;
};
} // end namespace AMDiS
\ No newline at end of file
#pragma once
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/umfpack.hh>
#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/linear_algebra/istl/Preconditioner.hpp>
#include <dune/amdis/linear_algebra/LinearSolverInterface.hpp>
#include <dune/amdis/linear_algebra/istl/ISTLRunner.hpp>
namespace AMDiS
{
// creators for preconditioners, depending on matrix-type
// A general solver class for istl solvers
// ---------------------------------------------------------------------------
template <class Matrix, class VectorX, class VectorB = VectorX>
class PreconCreator;
template <class... MatrixRows, class VectorX, class VectorB>
class PreconCreator<Dune::MultiTypeBlockMatrix<MatrixRows...>, VectorX, VectorB>
template <class Matrix, class VectorX, class VectorB, class Runner>
class LinearSolver
: public LinearSolverInterface<Matrix, VectorX, VectorB>
{
using BasePreconditioner = Dune::Preconditioner<VectorX, VectorB>;
using Super = LinearSolverInterface<Matrix, Vector, Vector>;
using RunnerBase = typename Super::RunnerBase;
public:
PreconCreator(std::string prefix)
: prefix(prefix)
/// Constructor
LinearSolver(std::string name, std::string prefix)
: runner(std::make_shared<Runner>(name, prefix))
{}
template <class Matrix>
std::unique_ptr<BasePreconditioner> create(std::string /*name*/, Matrix const& A)
{
double w = 1.0;
Parameters::get(prefix + "->relaxation", w);
AMDIS_WARNING("Only Richardson preconditioner available for MultiTypeBlockMatrix!");
return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w);
}
private:
std::string prefix;
};
template <class Block, class Alloc, class VectorX, class VectorB>
class PreconCreator<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>
{
using BasePreconditioner = Dune::Preconditioner<VectorX, VectorB>;
public:
PreconCreator(std::string prefix)
: prefix(prefix)
{}
template <class Matrix>
std::unique_ptr<BasePreconditioner> create(std::string name, Matrix const& A)
/// Implements \ref LinearSolverInterface::solveImpl()
virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
SolverInfo& solverInfo) override
{
double w = 1.0;
Parameters::get(prefix + "->relaxation", w);
if (name == "diag" || name == "jacobi")
return std::make_unique<Dune::SeqJac<Matrix, VectorX, VectorB>>(A, 1, w);
else if (name == "gs" || name == "gauss_seidel")
return std::make_unique<Dune::SeqGS<Matrix, VectorX, VectorB>>(A, 1, w);
else if (name == "sor")
return std::make_unique<Dune::SeqSOR<Matrix, VectorX, VectorB>>(A, 1, w);
else if (name == "ssor")
return std::make_unique<Dune::SeqSSOR<Matrix, VectorX, VectorB>>(A, 1, w);
else if (name == "ilu" || name == "ilu0")
return std::make_unique<Dune::SeqILU0<Matrix, VectorX, VectorB>>(A, w);
else
return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w);
}
private:
std::string prefix;
};
// creators for solvers, depending on matrix-type
// ---------------------------------------------------------------------------
template <class Matrix, class VectorX, class VectorB = VectorX>
class SolverCreator;
template <class... MatrixRows, class VectorX, class VectorB>
class SolverCreator<Dune::MultiTypeBlockMatrix<MatrixRows...>, VectorX, VectorB>
{
using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
public:
SolverCreator(std::string prefix)
: prefix(prefix)
, info(2)
{
Parameters::get(prefix + "->info", info);
}
template <class Matrix, class Precon>
std::unique_ptr<BaseSolver> create(std::string name, Matrix& A, Precon& P)
{
size_t maxIter = 500;
Parameters::get(prefix + "->max iteration", maxIter);
Timer t;
if (solverInfo.doCreateMatrixData()) {
// init matrix or wrap block-matrix or ...
runner->init(A);
}
double rTol = 1.e-6;
Parameters::get(prefix + "->reduction", rTol);
// solve the linear system
int error = runner->solve(A, x, b, solverInfo);
solverInfo.setError(error);
if (name == "cg")
return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "bicgstab")
return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "minres")
return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "gmres")
{
int restart = 30;
Parameters::get(prefix + "->restart", restart);
return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info);
}
else
AMDIS_ERROR_EXIT("Unknown solver name!");
if (!solverInfo.doStoreMatrixData())
runner->exit();
}
private:
std::string prefix;
int info;
private:
/// redirect the implementation to a runner. Implements a \ref RunnerInterface
std::shared_ptr<Runner> runner;
};
template <class Block, class Alloc, class VectorX, class VectorB>
class SolverCreator<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>
template <class Matrix, class VectorX, class VectorB>
class SolverCreator<Matrix, VectorX, VectorB>
{
using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
using Matrix = Dune::BCRSMatrix<Block,Alloc>;
using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
public:
SolverCreator(std::string prefix)
: prefix(prefix)
, info(2)
{
Parameters::get(prefix + "->info", info);
}
template <class Runner>
using Solver = LinearSolver<Matrix, VectorX, VectorB, Runner>;
template <class LinOperator, class Precon>
std::unique_ptr<BaseSolver> create(std::string name, LinOperator& A, Precon& P)
{
size_t maxIter = 500;
Parameters::get(prefix + "->max iteration", maxIter);
/// Instantiate a new linear solver
static std::unique_ptr<SolverBase> create(std::string name, std::string prefix)
{
using ISTL = ISTLRunner<Matrix, VectorX, VectorB>;
double rTol = 1.e-6;
Parameters::get(prefix + "->reduction", rTol);
if (name == "cg")
return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "bicgstab")
return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "minres")
return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
else if (name == "gmres")
{
int restart = 30;
Parameters::get(prefix + "->restart", restart);
return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info);
}
#if HAVE_SUITESPARSE_UMFPACK
else if (name == "umfpack" || name == "direct")
return std::make_unique<Dune::UMFPack<Matrix>>(A, info);
#endif
if (name == "cg" || name == "bicgstab" || name == "minres" || name == "gmres")
return std::make_unique<Solver<ISTL>>(name, prefix);
else
AMDIS_ERROR_EXIT("Unknown solver name!");
AMDIS_ERROR_EXIT("Unknown solver name!"); // TODO: add direct solver interface here
}
private:
std::string prefix;
int info;
};
// A general solver class for istl solvers
// ---------------------------------------------------------------------------
template <class Matrix, class VectorX, class VectorB = VectorX>
class LinearSolver
{
using BasePreconditioner = Dune::Preconditioner<VectorX, VectorB>;
using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
public:
LinearSolver(std::string name)
: name(name)
, solverName("cg")
, preconName("no")
{
// get creator string for solver
Parameters::get(name + "->name", solverName);
// get creator string for preconditioner
Parameters::get(name + "->left precon", preconName);
if (preconName.empty() || preconName == "no")
Parameters::get(name + "->right precon", preconName);
if (preconName.empty() || preconName == "no")
Parameters::get(name + "->precon->name", preconName);
}
void solve(AdaptInfo& adaptInfo, Matrix& A, VectorX& x, VectorB& b)
{
// create preconditioner
PreconCreator<Matrix, VectorX, VectorB> preconCreator(name + "->precon");
std::unique_ptr<BasePreconditioner> precon
= preconCreator.create(preconName, A);
// create a solver
Dune::MatrixAdapter<Matrix, VectorX, VectorB> linOperator(A);
PreconAdapter<BasePreconditioner> preconAdapter(*precon);
SolverCreator<Matrix, VectorX, VectorB> solverCreator(name);
std::unique_ptr<BaseSolver> solver
= solverCreator.create(solverName, linOperator, preconAdapter);
// storing some statistics
Dune::InverseOperatorResult statistics;
// solve the linear system
solver->apply(x, b, statistics);
adaptInfo.setSolverIterations(statistics.iterations);
adaptInfo.setSolverResidual(statistics.reduction);
}
private:
std::string name;
std::string solverName;
std::string preconName;
};
} // end namespace AMDiS
......@@ -5,6 +5,7 @@
#include <dune/istl/preconditioner.hh>
#include <dune/istl/solvercategory.hh>
#include <dune/amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
namespace AMDiS {
......@@ -22,27 +23,58 @@ namespace detail
} // end namespace detail
template <class Precon>
class PreconAdapter : public Dune::Preconditioner<typename Precon::domain_type, typename Precon::range_type>
template <class Matrix, class VectorX, class VectorB>
class ISTLPreconditioner
: public Dune::Preconditioner<VectorX, VectorB>
, public PreconditionerInterface<Matrix, VectorX, VectorB>
{
public:
enum {category=Dune::SolverCategory::sequential};
using X = typename Precon::domain_type;
using Y = typename Precon::range_type;
PreconAdapter(Precon& precon)
: precon(precon)
ISTLPreconditioner(std::string preconName, std::string prefix)
: preconName(preconName)
, prefix(prefix)