Skip to content

GitLab

Commit cfa2b1d9 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'feature/cleanup_mtl' into 'master' 

cleanup of headers and itl solvers in mtl backend

See merge request !124
parents 21dd5d0c 28b40d8f
Hide whitespace changes
Inline Side-by-side
cmake/modules/AmdisMacros.cmake
View file @ cfa2b1d9
......@@ -24,10 +24,14 @@ endif (NOT BACKEND)
if (BACKEND STREQUAL "MTL" OR BACKEND STREQUAL "MTL4")
find_package(MTL REQUIRED)
set(HAVE_MTL TRUE)
if (MTL_FOUND)
message(STATUS " Found MTL, version: ${MTL_VERSION}")
dune_register_package_flags(LIBRARIES MTL::MTL COMPILE_DEFINITIONS "ENABLE_MTL=1")
endif (MTL_FOUND)
message(STATUS " Found MTL, version: ${MTL_VERSION}")
dune_register_package_flags(LIBRARIES MTL::MTL COMPILE_DEFINITIONS "ENABLE_MTL=1")
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")
find_package(Eigen3 REQUIRED 3.3.5)
set(HAVE_EIGEN TRUE)
......
cmake/modules/CMakeLists.txt
View file @ cfa2b1d9
......@@ -2,6 +2,7 @@ install(FILES
AddAmdisExecutable.cmake
AmdisMacros.cmake
AmdisCXXFeatures.cmake
FindHYPRE.cmake
FindMTL.cmake
FindPETSc.cmake
PkgConfigLinkLibraries.cmake
......
cmake/modules/FindHYPRE.cmake 0 → 100644
View file @ cfa2b1d9
# 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()
cmake/modules/FindMTL.cmake
View file @ cfa2b1d9
......@@ -18,16 +18,15 @@
mark_as_advanced(MTL_FOUND MTL_COMPILE_DEFINITIONS MTL_INCLUDE_DIR)
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
/opt/software/mtl
/opt/sources/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)
......@@ -40,6 +39,10 @@ find_package_handle_standard_args(MTL
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)
add_library(MTL::MTL INTERFACE IMPORTED)
......@@ -48,8 +51,12 @@ if(MTL_FOUND AND NOT TARGET MTL::MTL)
list(APPEND MTL_COMPILE_DEFINITIONS "MTL_HAS_UMFPACK")
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
INTERFACE_INCLUDE_DIRECTORIES "${MTL_INCLUDE_DIR}"
INTERFACE_COMPILE_DEFINITIONS "${MTL_COMPILE_DEFINITIONS}"
)
INTERFACE_COMPILE_DEFINITIONS "${MTL_COMPILE_DEFINITIONS}")
endif()
config.h.cmake
View file @ cfa2b1d9
......@@ -43,11 +43,14 @@
/* Define to true if the MTL library is available */
#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 */
#cmakedefine HAVE_EIGEN 1
/* 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 */
#cmakedefine AMDIS_HAS_CXX_FOLD_EXPRESSIONS 1
......
src/amdis/linearalgebra/mtl/CMakeLists.txt
View file @ cfa2b1d9
install(FILES
Constraints.hpp
HyprePrecon.hpp
ITL_Preconditioner.hpp
ITL_Solver.hpp
KrylovRunner.hpp
MatrixBackend.hpp
Preconditioner.hpp
SolverPrecon.hpp
Traits.hpp
UmfpackRunner.hpp
VectorBackend.hpp
......
src/amdis/linearalgebra/mtl/HyprePrecon.hpp 0 → 100644
View file @ cfa2b1d9
#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
src/amdis/linearalgebra/mtl/ITL_Preconditioner.hpp
View file @ cfa2b1d9
......@@ -6,10 +6,15 @@
#include <boost/numeric/itl/pc/ic_0.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/SolverPrecon.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
{
......@@ -73,6 +78,7 @@ namespace AMDiS
template <class Matrix>
using ICPreconditioner = itl::pc::ic_0<Matrix>;
template <class Traits>
class DefaultCreators< PreconditionerInterface<Traits> >
{
......@@ -107,20 +113,28 @@ namespace AMDiS
Map::addCreator("identity", 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);
}
};
template <class Traits>
itl::pc::solver<PreconditionerInterface<Traits>, typename Traits::Vec, false>
solve(PreconditionerInterface<Traits> const& P, typename Traits::Vec const& vin)
template <class Traits, class Vector>
itl::pc::solver<PreconditionerInterface<Traits>, Vector, false>
solve(PreconditionerInterface<Traits> const& P, Vector const& vin)
{
return {P, vin};
}
template <class Traits>
itl::pc::solver<PreconditionerInterface<Traits>, typename Traits::Vec, true>
adjoint_solve(PreconditionerInterface<Traits> const& P, typename Traits::Vec const& vin)
template <class Traits, class Vector>
itl::pc::solver<PreconditionerInterface<Traits>, Vector, true>
adjoint_solve(PreconditionerInterface<Traits> const& P, Vector const& vin)
{
return {P, vin};
}
......
src/amdis/linearalgebra/mtl/ITL_Solver.hpp
View file @ cfa2b1d9
......@@ -36,10 +36,10 @@ namespace AMDiS
{
public:
explicit cg_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
template <class LinOp, class X, class B, class P, class I>
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
{
public:
explicit cgs_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
template <class LinOp, class X, class B, class P, class I>
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
{
public:
explicit bicg_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
template <class LinOp, class X, class B, class P, class I>
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
{
public:
explicit bicgstab_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
template <class LinOp, class X, class B, class P, class I>
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
{
public:
explicit bicgstab2_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const&, I& iter)
template <class LinOp, class X, class B, class P, class I>
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
{
public:
explicit qmr_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
template <class LinOp, class X, class B, class P, class I>
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
* Quasi-Minimal Residual method \implements ITL_Solver
*
* Solves a linear system by the Transposed-Free Quasi-Minimal Residual method
* (TFQMR). Does not use preconditioning currently.
* (TFQMR).
*/
class tfqmr_solver_type
{
public:
explicit tfqmr_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
template <class LinOp, class X, class B, class P, class I>
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
{
Parameters::get(name + "->ell", ell);
}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
return itl::bicgstab_ell(A, x, b, l, r, iter, ell);
itl::pc::identity<LinOp> id(A);
return itl::bicgstab_ell(A, x, b, precon, id, iter, ell);
}
};
......@@ -205,17 +208,18 @@ namespace AMDiS
Parameters::get(name + "->restart", restart);
Parameters::get(name + "->orthogonalization", ortho);
}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
itl::pc::identity<LinOp> id(A);
switch (ortho)
{
default:
case GRAM_SCHMIDT:
return itl::gmres2(A, x, b, l, r, iter, restart);
return itl::gmres2(A, x, b, id, precon, iter, restart);
break;
case HOUSEHOLDER:
return itl::gmres_householder(A, x, b, l, iter, restart);
return itl::gmres_householder(A, x, b, precon, iter, restart);
break;
}
}
......@@ -240,6 +244,8 @@ namespace AMDiS
* Peter Sonneveld and Martin B. van Gijzen, IDR(s): a family of simple and fast
* algorithms for solving large nonsymmetric linear systems.
* SIAM J. Sci. Comput. Vol. 31, No. 2, pp. 1035-1062 (2008). (copyright SIAM)
*
* NOTE: preconditioners are ignored.
*/
class idr_s_type
{
......@@ -249,10 +255,11 @@ namespace AMDiS
{
Parameters::get(name + "->s", s);
}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, P const&, I& iter)
{
return itl::idr_s(A, x, b, l, r, iter, s);
itl::pc::identity<LinOp> id(A);
return itl::idr_s(A, x, b, id, id, iter, s);
}
};
......@@ -269,10 +276,10 @@ namespace AMDiS
{
public:
explicit minres_solver_type(std::string const& /*name*/) {}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
return itl::minres(A, x, b, l, r, iter);
return itl::minres(A, x, b, precon, iter);
}
};
......@@ -295,10 +302,11 @@ namespace AMDiS
{
Parameters::get(name + "->restart", restart);
}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const& b, L const& l, R const& r, I& iter)
template <class LinOp, class X, class B, class P, class I>
int operator()(LinOp const& A, X& x, B const& b, P const& precon, I& iter)
{
return itl::gcr(A, x, b, l, r, iter, restart);
itl::pc::identity<LinOp> id(A);
return itl::gcr(A, x, b, id, precon, iter, restart);
}
};
......@@ -322,17 +330,18 @@ namespace AMDiS
Parameters::get(name + "->restart", restart);
Parameters::get(name + "->orthogonalization", orthogonalization);
}
template <class LinOp, class X, class B, class L, class R, class I>
int operator()(LinOp const& A, X& x, B const