Commit 28b40d8f authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

cleanup of headers and itl solvers in mtl backend

parent 21dd5d0c
...@@ -24,10 +24,14 @@ endif (NOT BACKEND) ...@@ -24,10 +24,14 @@ endif (NOT BACKEND)
if (BACKEND STREQUAL "MTL" OR BACKEND STREQUAL "MTL4") if (BACKEND STREQUAL "MTL" OR BACKEND STREQUAL "MTL4")
find_package(MTL REQUIRED) find_package(MTL REQUIRED)
set(HAVE_MTL TRUE) set(HAVE_MTL TRUE)
if (MTL_FOUND) message(STATUS " Found MTL, version: ${MTL_VERSION}")
message(STATUS " Found MTL, version: ${MTL_VERSION}") dune_register_package_flags(LIBRARIES MTL::MTL COMPILE_DEFINITIONS "ENABLE_MTL=1")
dune_register_package_flags(LIBRARIES MTL::MTL COMPILE_DEFINITIONS "ENABLE_MTL=1")
endif (MTL_FOUND) find_package(HYPRE)
if (HYPRE_FOUND)
set(HAVE_HYPRE TRUE)
dune_register_package_flags(LIBRARIES HYPRE::HYPRE COMPILE_DEFINITIONS "ENABLE_HYPRE=1")
endif(HYPRE_FOUND)
elseif (BACKEND STREQUAL "EIGEN" OR BACKEND STREQUAL "EIGEN3") elseif (BACKEND STREQUAL "EIGEN" OR BACKEND STREQUAL "EIGEN3")
find_package(Eigen3 REQUIRED 3.3.5) find_package(Eigen3 REQUIRED 3.3.5)
set(HAVE_EIGEN TRUE) set(HAVE_EIGEN TRUE)
......
...@@ -2,6 +2,7 @@ install(FILES ...@@ -2,6 +2,7 @@ install(FILES
AddAmdisExecutable.cmake AddAmdisExecutable.cmake
AmdisMacros.cmake AmdisMacros.cmake
AmdisCXXFeatures.cmake AmdisCXXFeatures.cmake
FindHYPRE.cmake
FindMTL.cmake FindMTL.cmake
FindPETSc.cmake FindPETSc.cmake
PkgConfigLinkLibraries.cmake PkgConfigLinkLibraries.cmake
......
# FindHYPRE.cmake
#
# Finds the HYPRE library
#
# This will define the following variables
#
# HYPRE_FOUND
# HYPRE_INCLUDE_DIR
# HYPRE_LIBRARIES
#
# and the following imported targets
#
# HYPRE::HYPRE
#
# Author: Simon Praetorius <simon.praetorius@tu-dresden.de>
mark_as_advanced(HYPRE_FOUND HYPRE_INCLUDE_DIR HYPRE_LIBRARIES)
find_path(HYPRE_INCLUDE_DIR HYPRE.h
PATHS
/opt/software/hypre
${HYPRE_DIR} ${HYPRE_ROOT}
ENV HYPRE_DIR ENV HYPRE_ROOT
ENV EBROOTHYPRE
PATH_SUFFIXES
hypre include include/hypre
NO_DEFAULT_PATH)
find_path(HYPRE_INCLUDE_DIR HYPRE.h
PATH_SUFFIXES hypre include include/hypre)
find_library(HYPRE_LIBRARY HYPRE
PATHS
/opt/software/hypre
${HYPRE_DIR} ${HYPRE_ROOT}
ENV HYPRE_DIR ENV HYPRE_ROOT
ENV EBROOTHYPRE
PATH_SUFFIXES
lib
NO_DEFAULT_PATH)
find_library(HYPRE_LIBRARY HYPRE
PATH_SUFFIXES lib)
if (HYPRE_LIBRARY)
set(HYPRE_LIBRARIES ${HYPRE_LIBRARY})
get_filename_component(HYPRE_LIBRARY_DIR ${HYPRE_LIBRARY} DIRECTORY)
foreach(_lib_name "HYPRE_FEI" "HYPRE_core")
find_library(_lib ${_lib_name} HINTS ${HYPRE_LIBRARY_DIR} NO_DEFAULT_PATH)
if (_lib)
list(APPEND HYPRE_LIBRARIES ${_lib})
endif (_lib)
unset(_lib)
endforeach()
unset(_lib_name)
endif (HYPRE_LIBRARY)
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(HYPRE
REQUIRED_VARS HYPRE_INCLUDE_DIR HYPRE_LIBRARY
)
# text for feature summary
set_package_properties("HYPRE" PROPERTIES
DESCRIPTION "high performance preconditioners")
if(HYPRE_FOUND AND NOT TARGET HYPRE::HYPRE)
add_library(HYPRE::HYPRE INTERFACE IMPORTED)
set_target_properties(HYPRE::HYPRE PROPERTIES
INTERFACE_INCLUDE_DIRECTORIES "${HYPRE_INCLUDE_DIR}"
INTERFACE_LINK_LIBRARIES "${HYPRE_LIBRARIES}"
)
endif()
...@@ -18,16 +18,15 @@ ...@@ -18,16 +18,15 @@
mark_as_advanced(MTL_FOUND MTL_COMPILE_DEFINITIONS MTL_INCLUDE_DIR) mark_as_advanced(MTL_FOUND MTL_COMPILE_DEFINITIONS MTL_INCLUDE_DIR)
find_path(MTL_INCLUDE_DIR boost/numeric/mtl/mtl.hpp find_path(MTL_INCLUDE_DIR boost/numeric/mtl/mtl.hpp
HINTS
${CMAKE_CURRENT_SOURCE_DIR}/../../libs/mtl4
${CMAKE_CURRENT_SOURCE_DIR}/../../include/mtl4
${MTL_DIR}/../../include
$ENV{MTL_DIR}/../../include
ENV EBROOTMTL
PATHS PATHS
/opt/software/mtl /opt/software/mtl
/opt/sources/mtl /opt/sources/mtl
/opt/development/mtl /opt/development/mtl
${MTL_DIR} ${MTL_ROOT}
ENV MTL_DIR ENV MTL_ROOT
ENV EBROOTMTL
PATH_SUFFIXES
include
) )
set(MTL_COMPILE_DEFINITIONS MTL_ASSERT_FOR_THROW=1) set(MTL_COMPILE_DEFINITIONS MTL_ASSERT_FOR_THROW=1)
...@@ -40,6 +39,10 @@ find_package_handle_standard_args(MTL ...@@ -40,6 +39,10 @@ find_package_handle_standard_args(MTL
REQUIRED_VARS MTL_INCLUDE_DIR REQUIRED_VARS MTL_INCLUDE_DIR
) )
# text for feature summary
set_package_properties("MTL" PROPERTIES
DESCRIPTION "Matrix Template Library")
if(MTL_FOUND AND NOT TARGET MTL::MTL) if(MTL_FOUND AND NOT TARGET MTL::MTL)
add_library(MTL::MTL INTERFACE IMPORTED) add_library(MTL::MTL INTERFACE IMPORTED)
...@@ -48,8 +51,12 @@ if(MTL_FOUND AND NOT TARGET MTL::MTL) ...@@ -48,8 +51,12 @@ if(MTL_FOUND AND NOT TARGET MTL::MTL)
list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK") list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
endif (SuiteSparse_FOUND) endif (SuiteSparse_FOUND)
find_package(HYPRE QUIET)
if (HYPRE_FOUND)
list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_HYPRE")
endif (HYPRE_FOUND)
set_target_properties(MTL::MTL PROPERTIES set_target_properties(MTL::MTL PROPERTIES
INTERFACE_INCLUDE_DIRECTORIES "${MTL_INCLUDE_DIR}" INTERFACE_INCLUDE_DIRECTORIES "${MTL_INCLUDE_DIR}"
INTERFACE_COMPILE_DEFINITIONS "${MTL_COMPILE_DEFINITIONS}" INTERFACE_COMPILE_DEFINITIONS "${MTL_COMPILE_DEFINITIONS}")
)
endif() endif()
...@@ -43,11 +43,14 @@ ...@@ -43,11 +43,14 @@
/* Define to true if the MTL library is available */ /* Define to true if the MTL library is available */
#cmakedefine HAVE_MTL 1 #cmakedefine HAVE_MTL 1
/* Define to true if the HYPRE library is available */
#cmakedefine HAVE_HYPRE ENABLE_HYPRE
/* Define to true if the Eigen3 library is available */ /* Define to true if the Eigen3 library is available */
#cmakedefine HAVE_EIGEN 1 #cmakedefine HAVE_EIGEN 1
/* Define to true if the PETSc library is available */ /* Define to true if the PETSc library is available */
#cmakedefine HAVE_PETSC ENABLE_PETSC #cmakedefine HAVE_PETSC 1
/* some detected compiler features may be used in AMDiS */ /* some detected compiler features may be used in AMDiS */
#cmakedefine AMDIS_HAS_CXX_FOLD_EXPRESSIONS 1 #cmakedefine AMDIS_HAS_CXX_FOLD_EXPRESSIONS 1
......
install(FILES install(FILES
Constraints.hpp Constraints.hpp
HyprePrecon.hpp
ITL_Preconditioner.hpp ITL_Preconditioner.hpp
ITL_Solver.hpp ITL_Solver.hpp
KrylovRunner.hpp KrylovRunner.hpp
MatrixBackend.hpp MatrixBackend.hpp
Preconditioner.hpp Preconditioner.hpp
SolverPrecon.hpp
Traits.hpp Traits.hpp
UmfpackRunner.hpp UmfpackRunner.hpp
VectorBackend.hpp VectorBackend.hpp
......
#pragma once
#if HAVE_HYPRE
#include <dune/common/unused.hh>
#include <amdis/linearalgebra/mtl/itl/hypre.hpp>
namespace AMDiS
{
/// \brief Interface to the HYPRE BoomerAMG preconditioner [...]
/**
* Parameters provided by AMDiS:
*
* [PRECON]->cycle mode:
* 1...V-cycle
* 2...W-cycle
*
* [PRECON]->interpolation type:
* 0...classical modified interpolation
* 1...LS interpolation (for use with GSMG)
* 2...classical modified interpolation for hyperbolic PDEs
* 3...direct interpolation (with separation of weights)
* 4...multipass interpolation
* 5...multipass interpolation (with separation of weights)
* 6...extended+i interpolation
* 7...extended+i (if no common C neighbor) interpolation
* 8...standard interpolation
* 9...standard interpolation (with separation of weights)
* 10..classical block interpolation (for use with nodal systems version only)
* 11..classical block interpolation (for use with nodal systems version only)
* with diagonalized diagonal blocks
* 12..FF interpolation
* 13..FF1 interpolation
* 14..extended interpolation
*
* [PRECON]->info:
* 0...no printout (default)
* 1...print setup information
* 2...print solve information
* 3...print both setup and solve information
*
* [PRECON]->relax type:
* 0...Jacobi
* 1...Gauss-Seidel, sequential (very slow!)
* 2...Gauss-Seidel, interior points in parallel, boundary sequential (slow!)
* 3...hybrid Gauss-Seidel or SOR, forward solve
* 4...hybrid Gauss-Seidel or SOR, backward solve
* 5...hybrid chaotic Gauss-Seidel (works only with OpenMP)
* 6...hybrid symmetric Gauss-Seidel or SSOR
* 9...Gaussian elimination (only on coarsest level)
**/
template <class Traits>
class HyprePrecon
: public PreconditionerInterface<Traits>
{
using Self = HyprePrecon;
using Matrix = typename Traits::Mat;
using Vector = typename Traits::Vec;
using Super = PreconditionerInterface<Traits>;
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:
HyprePrecon(std::string const& prefix)
{
// read HYPRE parameters
int maxIter = 1, cycleMode = -1, interpolation = -1, relaxation = -1, info = 1;
Parameters::get(prefix + "->max iteration", maxIter);
Parameters::get(prefix + "->cycle mode", cycleMode);
Parameters::get(prefix + "->interpolation type", interpolation);
Parameters::get(prefix + "->relax type", relaxation);
Parameters::get(prefix + "->info", info);
config_.maxIter(maxIter);
config_.interpolationType(interpolation);
config_.relaxType(relaxation);
config_.cycle(cycleMode);
config_.printLevel(info);
}
/// Implementation of \ref PreconditionerBase::init()
void init(Matrix const& fullMatrix) override
{
hypreMatrix_.init(fullMatrix);
HYPRE_IJMatrixGetObject(hypreMatrix_, (void**)(&matrix_));
HYPRE_BoomerAMGCreate(&solver_);
// apply solver parameters
config_(solver_);
mtl::dense_vector<double> swap(1, 0.0);
mtl::HypreParVector x(swap);
HYPRE_BoomerAMGSetup(solver_, matrix_, x, x);
solverCreated_ = true;
}
/// Implementation of \ref PreconditionerInterface::exit()
void exit() override
{
if (solverCreated_)
HYPRE_BoomerAMGDestroy(solver_);
solverCreated_ = false;
}
/// Implementation of \ref PreconditionerBase::solve()
void solve(Vector const& mtlB, Vector& mtlX) const override
{
assert(solverCreated_);
mtl::HypreParVector x(mtlB); // use b as initial guess
mtl::HypreParVector b(mtlB);
DUNE_UNUSED int error = HYPRE_BoomerAMGSolve(solver_, matrix_, b, x);
assert(error != 0);
// write output back to MTL vector
mtl::convert(x.getHypreVector(), mtlX);
}
/// Implementation of \ref PreconditionerBase::adjoint_solve()
void adjoint_solve(Vector const& mtlB, Vector& mtlX) const override
{
assert(solverCreated_);
mtl::HypreParVector x(mtlB); // use b as initial guess
mtl::HypreParVector b(mtlB);
DUNE_UNUSED int error = HYPRE_BoomerAMGSolveT(solver_, matrix_, b, x);
assert(error != 0);
// write output back to MTL vector
mtl::convert(x.getHypreVector(), mtlX);
}
private:
mtl::HypreMatrix hypreMatrix_;
HYPRE_ParCSRMatrix matrix_;
HYPRE_Solver solver_;
mtl::AMGConfigurator config_;
bool solverCreated_ = false;
};
} // end namespace AMDiS
#endif // HAVE_HYPRE
...@@ -6,10 +6,15 @@ ...@@ -6,10 +6,15 @@
#include <boost/numeric/itl/pc/ic_0.hpp> #include <boost/numeric/itl/pc/ic_0.hpp>
#include <boost/numeric/mtl/vector/assigner.hpp> #include <boost/numeric/mtl/vector/assigner.hpp>
#include <amdis/linearalgebra/mtl/itl/masslumping.hpp> #include <amdis/CreatorMap.hpp>
#include <amdis/linearalgebra/mtl/Preconditioner.hpp> #include <amdis/linearalgebra/mtl/Preconditioner.hpp>
#include <amdis/linearalgebra/mtl/SolverPrecon.hpp>
#include <amdis/linearalgebra/mtl/Traits.hpp> #include <amdis/linearalgebra/mtl/Traits.hpp>
#include <amdis/CreatorMap.hpp> #include <amdis/linearalgebra/mtl/itl/masslumping.hpp>
#if HAVE_HYPRE
#include <amdis/linearalgebra/mtl/HyprePrecon.hpp>
#endif
namespace AMDiS namespace AMDiS
{ {
...@@ -73,6 +78,7 @@ namespace AMDiS ...@@ -73,6 +78,7 @@ namespace AMDiS
template <class Matrix> template <class Matrix>
using ICPreconditioner = itl::pc::ic_0<Matrix>; using ICPreconditioner = itl::pc::ic_0<Matrix>;
template <class Traits> template <class Traits>
class DefaultCreators< PreconditionerInterface<Traits> > class DefaultCreators< PreconditionerInterface<Traits> >
{ {
...@@ -107,20 +113,28 @@ namespace AMDiS ...@@ -107,20 +113,28 @@ namespace AMDiS
Map::addCreator("identity", pc_id); Map::addCreator("identity", pc_id);
Map::addCreator("no", pc_id); Map::addCreator("no", pc_id);
auto pc_solver = new typename SolverPrecon<Traits>::Creator;
Map::addCreator("solver", pc_solver);
#if HAVE_HYPRE
auto pc_hypre = new typename HyprePrecon<Traits>::Creator;
Map::addCreator("hypre", pc_hypre);
#endif
Map::addCreator("default", pc_id); Map::addCreator("default", pc_id);
} }
}; };
template <class Traits> template <class Traits, class Vector>
itl::pc::solver<PreconditionerInterface<Traits>, typename Traits::Vec, false> itl::pc::solver<PreconditionerInterface<Traits>, Vector, false>
solve(PreconditionerInterface<Traits> const& P, typename Traits::Vec const& vin) solve(PreconditionerInterface<Traits> const& P, Vector const& vin)
{ {
return {P, vin}; return {P, vin};
} }
template <class Traits> template <class Traits, class Vector>
itl::pc::solver<PreconditionerInterface<Traits>, typename Traits::Vec, true> itl::pc::solver<PreconditionerInterface<Traits>, Vector, true>
adjoint_solve(PreconditionerInterface<Traits> const& P, typename Traits::Vec const& vin) adjoint_solve(PreconditionerInterface<Traits> const& P, Vector const& vin)
{ {
return {P, vin}; return {P, vin};
} }
......
...@@ -36,10 +36,10 @@ namespace AMDiS ...@@ -36,10 +36,10 @@ namespace AMDiS
{ {
public: public:
explicit cg_solver_type(std::string const& /*name*/) {} explicit cg_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I> template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter) int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{ {
return itl::cg(A, x, b, l, r, iter); return itl::cg(A, x, b, precon, iter);
} }
}; };
...@@ -56,10 +56,10 @@ namespace AMDiS ...@@ -56,10 +56,10 @@ namespace AMDiS
{ {
public: public:
explicit cgs_solver_type(std::string const& /*name*/) {} explicit cgs_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I> template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter) int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{ {
return itl::cgs(A, x, b, l, iter); return itl::cgs(A, x, b, precon, iter);
} }
}; };
...@@ -76,10 +76,10 @@ namespace AMDiS ...@@ -76,10 +76,10 @@ namespace AMDiS
{ {
public: public:
explicit bicg_solver_type(std::string const& /*name*/) {} explicit bicg_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I> template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter) int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{ {
return itl::bicg(A, x, b, l, iter); return itl::bicg(A, x, b, precon, iter);
} }
}; };
...@@ -96,10 +96,10 @@ namespace AMDiS ...@@ -96,10 +96,10 @@ namespace AMDiS
{ {
public: public:
explicit bicgstab_type(std::string const& /*name*/) {} explicit bicgstab_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I> template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter) int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{ {
return itl::bicgstab(A, x, b, l, iter); return itl::bicgstab(A, x, b, precon, iter);
} }
}; };
...@@ -116,10 +116,10 @@ namespace AMDiS ...@@ -116,10 +116,10 @@ namespace AMDiS
{ {
public: public:
explicit bicgstab2_type(std::string const& /*name*/) {} explicit bicgstab2_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I> template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter) int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{ {
return itl::bicgstab_2(A, x, b, l, iter); return itl::bicgstab_2(A, x, b, precon, iter);
} }
}; };
...@@ -135,10 +135,11 @@ namespace AMDiS ...@@ -135,10 +135,11 @@ namespace AMDiS
{ {
public: public:
explicit qmr_solver_type(std::string const& /*name*/) {} explicit qmr_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I> template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter) int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{ {
return itl::qmr(A, x, b, l, r, iter); itl::pc::identity<LinOp> id(A);
return itl::qmr(A, x, b, precon, id, iter);
} }
}; };
...@@ -149,16 +150,17 @@ namespace AMDiS ...@@ -149,16 +150,17 @@ namespace AMDiS
* Quasi-Minimal Residual method \implements ITL_Solver * Quasi-Minimal Residual method \implements ITL_Solver
* *
* Solves a linear system by the Transposed-Free Quasi-Minimal Residual method * Solves a linear system by the Transposed-Free Quasi-Minimal Residual method
* (TFQMR). Does not use preconditioning currently. * (TFQMR).
*/ */
class tfqmr_solver_type class tfqmr_solver_type
{ {
public: public:
explicit tfqmr_solver_type(std::string const& /*name*/) {} explicit tfqmr_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I> template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter) int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{ {
return itl::tfqmr(A, x, b, l, r, iter); itl::pc::identity<LinOp> id(A);
return itl::tfqmr(A, x, b, precon, id, iter);
} }
}; };
...@@ -179,10 +181,11 @@ namespace AMDiS ...@@ -179,10 +181,11 @@ namespace AMDiS
{ {
Parameters::get(name + "->ell", ell); Parameters::get(name + "->ell", ell);