Commit 07cf0c03 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Simplify linear algebra backend LinearSolver

parent ac9e1f4c
...@@ -41,7 +41,7 @@ namespace AMDiS ...@@ -41,7 +41,7 @@ namespace AMDiS
} }
/// Creates a object of the type corresponding to key. /// Creates a object of the type corresponding to key.
static CreatorInterface<BaseClass>* getCreator(std::string key, std::string initFileStr) static CreatorInterface<BaseClass>* get(std::string key, std::string initFileStr)
{ {
init(); init();
......
...@@ -3,8 +3,9 @@ ...@@ -3,8 +3,9 @@
#if AMDIS_BACKEND == AMDIS_BACKEND_MTL #if AMDIS_BACKEND == AMDIS_BACKEND_MTL
#include <amdis/linearalgebra/mtl/Constraints.hpp> #include <amdis/linearalgebra/mtl/Constraints.hpp>
#include <amdis/linearalgebra/mtl/ITL_Solver.hpp> #include <amdis/linearalgebra/mtl/LinearSolver.hpp>
#include <amdis/linearalgebra/mtl/ITL_Preconditioner.hpp> #include <amdis/linearalgebra/mtl/Preconditioners.hpp>
#include <amdis/linearalgebra/mtl/Solvers.hpp>
#include <amdis/linearalgebra/mtl/Traits.hpp> #include <amdis/linearalgebra/mtl/Traits.hpp>
#include <amdis/linearalgebra/mtl/MatrixBackend.hpp> #include <amdis/linearalgebra/mtl/MatrixBackend.hpp>
#include <amdis/linearalgebra/mtl/VectorBackend.hpp> #include <amdis/linearalgebra/mtl/VectorBackend.hpp>
...@@ -12,7 +13,8 @@ ...@@ -12,7 +13,8 @@
#elif AMDIS_BACKEND == AMDIS_BACKEND_EIGEN #elif AMDIS_BACKEND == AMDIS_BACKEND_EIGEN
#include <amdis/linearalgebra/eigen/Constraints.hpp> #include <amdis/linearalgebra/eigen/Constraints.hpp>
#include <amdis/linearalgebra/eigen/SolverCreator.hpp> #include <amdis/linearalgebra/eigen/LinearSolver.hpp>
#include <amdis/linearalgebra/eigen/Solvers.hpp>
#include <amdis/linearalgebra/eigen/Traits.hpp> #include <amdis/linearalgebra/eigen/Traits.hpp>
#include <amdis/linearalgebra/eigen/MatrixBackend.hpp> #include <amdis/linearalgebra/eigen/MatrixBackend.hpp>
#include <amdis/linearalgebra/eigen/VectorBackend.hpp> #include <amdis/linearalgebra/eigen/VectorBackend.hpp>
...@@ -20,7 +22,7 @@ ...@@ -20,7 +22,7 @@
#elif AMDIS_BACKEND == AMDIS_BACKEND_PETSC #elif AMDIS_BACKEND == AMDIS_BACKEND_PETSC
#include <amdis/linearalgebra/petsc/Constraints.hpp> #include <amdis/linearalgebra/petsc/Constraints.hpp>
#include <amdis/linearalgebra/petsc/SolverCreator.hpp> #include <amdis/linearalgebra/petsc/LinearSolver.hpp>
#include <amdis/linearalgebra/petsc/Traits.hpp> #include <amdis/linearalgebra/petsc/Traits.hpp>
#include <amdis/linearalgebra/petsc/MatrixBackend.hpp> #include <amdis/linearalgebra/petsc/MatrixBackend.hpp>
#include <amdis/linearalgebra/petsc/VectorBackend.hpp> #include <amdis/linearalgebra/petsc/VectorBackend.hpp>
...@@ -28,9 +30,9 @@ ...@@ -28,9 +30,9 @@
#elif AMDIS_BACKEND == AMDIS_BACKEND_ISTL #elif AMDIS_BACKEND == AMDIS_BACKEND_ISTL
#include <amdis/linearalgebra/istl/Constraints.hpp> #include <amdis/linearalgebra/istl/Constraints.hpp>
#include <amdis/linearalgebra/istl/ISTLSolver.hpp> #include <amdis/linearalgebra/istl/LinearSolver.hpp>
#include <amdis/linearalgebra/istl/PreconCreator.hpp> #include <amdis/linearalgebra/istl/Preconditioners.hpp>
#include <amdis/linearalgebra/istl/SolverCreator.hpp> #include <amdis/linearalgebra/istl/Solvers.hpp>
#include <amdis/linearalgebra/istl/Traits.hpp> #include <amdis/linearalgebra/istl/Traits.hpp>
#include <amdis/linearalgebra/istl/MatrixBackend.hpp> #include <amdis/linearalgebra/istl/MatrixBackend.hpp>
#include <amdis/linearalgebra/istl/VectorBackend.hpp> #include <amdis/linearalgebra/istl/VectorBackend.hpp>
...@@ -39,5 +41,4 @@ ...@@ -39,5 +41,4 @@
#include <amdis/linearalgebra/MatrixFacade.hpp> #include <amdis/linearalgebra/MatrixFacade.hpp>
#include <amdis/linearalgebra/VectorFacade.hpp> #include <amdis/linearalgebra/VectorFacade.hpp>
#include <amdis/linearalgebra/LinearSolver.hpp> #include <amdis/linearalgebra/LinearSolverInterface.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp>
...@@ -79,7 +79,7 @@ namespace AMDiS ...@@ -79,7 +79,7 @@ namespace AMDiS
using PartitionSet = typename LinAlgTraits::PartitionSet; using PartitionSet = typename LinAlgTraits::PartitionSet;
public: public:
using LinearSolver = LinearSolverInterface<Mat,Vec>; using LinearSolver = LinearSolverInterface<Mat,Vec,Vec>;
using SystemMatrix = BiLinearForm<GlobalBasis, GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>; using SystemMatrix = BiLinearForm<GlobalBasis, GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>;
using SystemVector = LinearForm<GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>; using SystemVector = LinearForm<GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>;
using SolutionVector = DOFVector<GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>; using SolutionVector = DOFVector<GlobalBasis, typename Traits::CoefficientType, LinAlgTraits>;
......
...@@ -228,10 +228,8 @@ void ProblemStat<Traits>::createSolver() ...@@ -228,10 +228,8 @@ void ProblemStat<Traits>::createSolver()
std::string solverName = "default"; std::string solverName = "default";
Parameters::get(name_ + "->solver", solverName); Parameters::get(name_ + "->solver", solverName);
auto solverCreator using Solver = AMDiS::LinearSolver<Mat,Vec,Vec>;
= named(CreatorMap<LinearSolver>::getCreator(solverName, name_ + "->solver")); linearSolver_ = std::make_shared<Solver>(solverName, name_ + "->solver");
linearSolver_ = solverCreator->createWithString(name_ + "->solver");
} }
...@@ -348,28 +346,26 @@ template <class Traits> ...@@ -348,28 +346,26 @@ template <class Traits>
void ProblemStat<Traits>:: void ProblemStat<Traits>::
solve(AdaptInfo& /*adaptInfo*/, bool createMatrixData, bool storeMatrixData) solve(AdaptInfo& /*adaptInfo*/, bool createMatrixData, bool storeMatrixData)
{ {
Dune::Timer t; Dune::InverseOperatorResult stat;
SolverInfo solverInfo(name_ + "->solver");
solverInfo.setCreateMatrixData(createMatrixData);
solverInfo.setStoreMatrixData(storeMatrixData);
solution_->resize(); solution_->resize();
linearSolver_->solve(*systemMatrix_, *solution_, *rhs_, solverInfo);
if (solverInfo.info() > 0) { if (createMatrixData)
msg("solution of discrete system needed {} seconds", t.elapsed()); linearSolver_->init(systemMatrix_->impl());
if (solverInfo.absResidual() >= 0.0) { // solve the linear system
if (solverInfo.relResidual() >= 0.0) linearSolver_->apply(solution_->impl(), rhs_->impl(), stat);
msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
solverInfo.absResidual(), solverInfo.relResidual()); if (!storeMatrixData)
else linearSolver_->finish();
msg("Residual norm: ||b-Ax|| = {}", solverInfo.absResidual());
} info(2, "solution of discrete system needed {} seconds", stat.elapsed);
} if (stat.reduction >= 0.0)
info(1, "Residual norm: ||b-Ax||/||b-Ax^0|| = {}", stat.reduction);
test_exit(!solverInfo.doBreak() || !solverInfo.error(), "Could not solver the linear system!"); bool ignoreConverged = false;
Parameters::get(name_ + "->solver->ignore converged", ignoreConverged);
test_exit(stat.converged || ignoreConverged, "Could not solver the linear system!");
} }
......
...@@ -5,12 +5,9 @@ install(FILES ...@@ -5,12 +5,9 @@ install(FILES
Constraints.hpp Constraints.hpp
DOFMapping.hpp DOFMapping.hpp
DOFMapping.inc.hpp DOFMapping.inc.hpp
LinearSolver.hpp
LinearSolverInterface.hpp LinearSolverInterface.hpp
MatrixFacade.hpp MatrixFacade.hpp
ParallelIndexSet.hpp ParallelIndexSet.hpp
RunnerInterface.hpp
SolverInfo.hpp
SparsityPattern.hpp SparsityPattern.hpp
SymmetryStructure.hpp SymmetryStructure.hpp
Traits.hpp Traits.hpp
......
#pragma once
#include <string>
#include <dune/common/timer.hh>
// AMDiS includes
#include <amdis/CreatorInterface.hpp>
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/LinearSolverInterface.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp>
namespace AMDiS
{
/** \ingroup Solver
*
* \brief
* Wrapper class for various solvers. These solvers
* are parametrized by MatrixType and VectorType. The different
* algorithms, like krylov subspace methods or other external
* solvers where the backend provides an interface, can be assigned
* by different Runner objects.
**/
template <class Mat, class Vec, class Runner>
class LinearSolver
: public LinearSolverInterface<Mat,Vec>
{
using Self = LinearSolver;
using Super = LinearSolverInterface<Mat,Vec>;
public:
/// A creator to be used instead of the constructor.
struct Creator : CreatorInterfaceName<Super>
{
std::unique_ptr<Super> createWithString(std::string prefix) override
{
return std::make_unique<Self>(prefix);
}
};
public:
/// Constructor
template <class... Args>
explicit LinearSolver(std::string prefix, Args&&... args)
: runner_(prefix, FWD(args)...)
{}
Runner* runner() override
{
return &runner_;
}
protected:
/// Implements \ref LinearSolverInterface::solveImpl()
void solveImpl(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo) override
{
Dune::Timer t;
if (solverInfo.doCreateMatrixData()) {
// init matrix or wrap block-matrix or ...
runner_.init(A.matrix());
}
if (solverInfo.info() > 0)
msg("prepare solver needed {} seconds", t.elapsed());
int error = runner_.solve(A.matrix(), x.vector(), b.vector(), solverInfo);
solverInfo.setError(error);
if (!solverInfo.doStoreMatrixData())
runner_.exit();
}
private:
/// redirect the implementation to a runner implementing a \ref RunnerInterface
Runner runner_;
};
} // end namespace AMDiS
#pragma once #pragma once
#include <amdis/Output.hpp> #include <dune/istl/solver.hh>
#include <amdis/linearalgebra/RunnerInterface.hpp>
/**
* \defgroup Solver Solver module
* @{ <img src="solver.png"> @}
*
* \brief
* Contains all classes needed for solving linear and non linear equation
* systems.
*/
namespace AMDiS namespace AMDiS
{ {
class SolverInfo; template <class M, class X, class Y = X>
template <class T, template <class> class MatrixImpl>
class MatrixFacade;
template <class T, template <class> class VectorImpl>
class VectorFacade;
/// Abstract base class for linear solvers
template <class Mat, class Vec>
class LinearSolverInterface class LinearSolverInterface
{ {
public: public:
/// Destructor. /// \brief Prepare the solve (and preconditioner), e.g. make a factorization
virtual ~LinearSolverInterface() = default; /// of the matrix, or extract its diagonal in a jacobian precon.
virtual void init(M const& A) = 0;
/// Public method to call in order to solve a linear system Ax = b.
/**
* The method redirects to the specific linear solver and prints statistics
* and error estimations at the end.
*
* The parameters correspond to
* \p A A [block-]matrix that represents the system-matrix.
* \p x A [block-]vector for the unknown components.
* \p b A [block-]vector for the right-hand side of the linear system.
**/
template <class TA, class TX, class TY, template <class> class MI, template <class> class VI>
void solve(MatrixFacade<TA,MI> const& A, VectorFacade<TX,VI>& x, VectorFacade<TY,VI> const& b, SolverInfo& solverInfo)
{
solveImpl(A.impl(), x.impl(), b.impl(), solverInfo);
}
virtual RunnerInterface<Mat,Vec>* runner() /// \brief Cleanup the solver, e.g. free the previously created factorization.
{ virtual void finish() = 0;
error_exit("Must be implemented by derived class.");
return nullptr;
}
protected: /// \brief Apply the inverse operator to the rhs vector b
/// main methods that all solvers must implement virtual void apply(X& x, Y const& b, Dune::InverseOperatorResult& res) = 0;
virtual void solveImpl(Mat const& A, Vec& x, Vec const& b, SolverInfo& solverInfo) = 0;
}; };
} // end namespace AMDiS } // end namespace AMDiS
#pragma once
#include <amdis/Output.hpp>
namespace AMDiS
{
class SolverInfo;
/// Interface for Runner / Worker types used in solver classes
template <class Mat, class Vec>
class RunnerInterface
{
using M = typename Mat::BaseMatrix;
using X = typename Vec::BaseVector;
using Y = typename Vec::BaseVector;
public:
/// virtual destructor
virtual ~RunnerInterface() = default;
/// Is called at the beginning of a solution procedure
virtual void init(M const& A) = 0;
/// Is called at the end of a solution procedure
virtual void exit() = 0;
/// Solve the system A*x = b
virtual int solve(M const& A, X& x, Y const& b, SolverInfo& solverInfo) = 0;
/// Solve the system A*x = b
virtual int adjointSolve(M const& A, X& x, Y const& b, SolverInfo& solverInfo)
{
error_exit("Must be implemented by derived class.");
return 0;
}
};
} // end namespace AMDiS
#pragma once
#include <string>
#include <amdis/Initfile.hpp>
namespace AMDiS
{
/// Class that stores information about the solution process, like tolerances
/// and iteration counts and is passed to the solver and filled there.
class SolverInfo
{
public:
/// The constructor reads needed parameters and sets solver \p prefix.
/**
* Reads parameters for a solver with name 'NAME':
* NAME->info \ref info_
* NAME->break if tolerance not reached \ref breakTolNotReached_
**/
explicit SolverInfo(std::string const& prefix)
: prefix_(prefix)
{
Parameters::get(prefix + "->info", info_);
Parameters::get(prefix + "->break if tolerance not reached", breakTolNotReached_);
}
/** \name getting methods
* \{
*/
/// Returns error code in last run of an iterative solver
int error() const
{
return error_;
}
/// Returns info
int info() const
{
return info_;
}
/// Returns \ref absResidual_
double absResidual() const
{
return absResidual_;
}
/// Returns \ref relResidual_
double relResidual() const
{
return relResidual_;
}
/// Returns \ref createMatrixData
bool doCreateMatrixData() const
{
return createMatrixData_;
}
/// Returns \ref storeMatrixData
bool doStoreMatrixData() const
{
return storeMatrixData_;
}
bool doBreak() const
{
return breakTolNotReached_;
}
/** \} */
/** \name setting methods
* \{
*/
/// Sets \ref absResidual_
void setAbsResidual(double r)
{
absResidual_ = r;
}
/// Sets \ref relResidual_
void setRelResidual(double r)
{
relResidual_ = r;
}
/// Sets \ref info_
void setInfo(int i)
{
info_ = i;
}
/// Sets \ref error_
void setError(int e)
{
error_ = e;
}
/// Sets \ref createMatrixData_
void setCreateMatrixData(bool b)
{
createMatrixData_ = b;
}
/// Sets \ref storeMatrixData_
void setStoreMatrixData(bool b)
{
storeMatrixData_ = b;
}
/** \} */
private:
/// The initfile prefix to read parameters
std::string prefix_;
/// Throw an error if tolerance could not be reached
bool breakTolNotReached_ = false;
/// The solver verbosity level
int info_ = 0;
/// The absolute residual, default (-1): not set
double absResidual_ = -1.0;
/// The relative residual, relative to the two_norm(b), default (-1): not set
double relResidual_ = -1.0;
/// The error-code, default (-1): not set
int error_ = -1;
/// If true, the matrix will be initialized and the
/// corresponding runner of the system receives the
/// matrix in the init() method.
bool createMatrixData_ = true;
/// If false, the exit() method of the runner will be called.
bool storeMatrixData_ = false;
};
} // end namespace AMDiS
...@@ -2,11 +2,12 @@ install(FILES ...@@ -2,11 +2,12 @@ install(FILES
Constraints.hpp Constraints.hpp
DirectRunner.hpp DirectRunner.hpp
IterativeRunner.hpp IterativeRunner.hpp
LinearSolver.hpp
MatrixBackend.hpp MatrixBackend.hpp
MatrixSize.hpp MatrixSize.hpp
PreconConfig.hpp PreconConfig.hpp
SolverConfig.hpp SolverConfig.hpp
SolverCreator.hpp Solvers.hpp
Traits.hpp Traits.hpp
VectorBackend.hpp VectorBackend.hpp
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linearalgebra/eigen) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linearalgebra/eigen)
#pragma once #pragma once
#include <algorithm> #include <memory>
#include <string> #include <string>
#include <amdis/linearalgebra/RunnerInterface.hpp> #include <amdis/CreatorInterface.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp> #include <amdis/Initfile.hpp>
#include <amdis/linearalgebra/LinearSolverInterface.hpp>
#include <amdis/linearalgebra/eigen/SolverConfig.hpp> #include <amdis/linearalgebra/eigen/SolverConfig.hpp>
namespace AMDiS namespace AMDiS
...@@ -12,19 +13,25 @@ namespace AMDiS ...@@ -12,19 +13,25 @@ namespace AMDiS
/** /**
* \ingroup Solver * \ingroup Solver
* \class AMDiS::DirectRunner * \class AMDiS::DirectRunner
* \brief \implements RunnerInterface for the (external) direct solvers * \brief \implements Dune::InverseOperator for the (external) direct solvers
*/ */
template <class Mat, class Vec, template <class> class Solver> template <class M, class X, class Y, template <class> class Solver>
class DirectRunner class DirectRunner
: public RunnerInterface<Mat,Vec> : public LinearSolverInterface<M,X,Y>
{ {
using M = typename Mat::BaseMatrix; using Self = DirectRunner;
using X = typename Vec::BaseVector;
using Y = typename Vec::BaseVector;
protected:
using EigenSolver = Solver<M>; using EigenSolver = Solver<M>;
public:
struct Creator final : CreatorInterfaceName<LinearSolverInterface<M,X,Y>>
{
std::unique_ptr<LinearSolverInterface<M,X,Y>>
createWithString(std::string prefix) final
{
return std::make_unique<Self>(prefix);
}
};
public: public:
/// Constructor. /// Constructor.
DirectRunner(std::string const& prefix) DirectRunner(std::string const& prefix)
...@@ -34,7 +41,7 @@ namespace AMDiS ...@@ -34,7 +41,7 @@ namespace AMDiS
Parameters::get(prefix + "->reuse pattern", reusePattern_); Parameters::get(prefix + "->reuse pattern", reusePattern_);
} }
/// Implements \ref RunnerInterface::init() /// initialize the matrix and maybe compute factorization