Commit 4ffe4a74 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Hypre-solver updated and new preconditioner added

parent d41f7355
...@@ -24,6 +24,9 @@ ...@@ -24,6 +24,9 @@
#ifdef HAVE_ZOLTAN #ifdef HAVE_ZOLTAN
#include <zoltan_cpp.h> #include <zoltan_cpp.h>
#endif #endif
#ifdef MTL_HAS_HYPRE
#include "mpi.h"
#endif
#include "boost/program_options.hpp" #include "boost/program_options.hpp"
namespace AMDiS { namespace AMDiS {
...@@ -36,23 +39,24 @@ namespace AMDiS { ...@@ -36,23 +39,24 @@ namespace AMDiS {
void init(int argc, char **argv, std::string initFileName) void init(int argc, char **argv, std::string initFileName)
{ {
#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_SEQ_PETSC)
#ifdef HAVE_PARALLEL_MTL4 #ifdef HAVE_PARALLEL_MTL4
mtl_environment = new mtl::par::environment(argc, argv); mtl_environment = new mtl::par::environment(argc, argv);
#else #elif defined HAVE_PETSC || defined HAVE_SEQ_PETSC || defined HAVE_PARALLEL_PETSC
PetscInitialize(&argc, &argv, NULL, NULL); PetscInitialize(&argc, &argv, NULL, NULL);
#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) #elif defined MTL_HAS_HYPRE
MPI_Init(&argc, &argv);
#endif
#if defined(HAVE_PARALLEL_DOMAIN_AMDIS)
Parallel::mpi::startRand(); Parallel::mpi::startRand();
#else #else
srand(time(0)); srand(time(0));
#endif #endif
#endif
#ifdef HAVE_ZOLTAN #ifdef HAVE_ZOLTAN
float zoltanVersion = 0.0; float zoltanVersion = 0.0;
Zoltan_Initialize(argc, argv, &zoltanVersion); Zoltan_Initialize(argc, argv, &zoltanVersion);
#endif #endif
#endif
Parameters::clearData(); Parameters::clearData();
...@@ -139,14 +143,15 @@ namespace AMDiS { ...@@ -139,14 +143,15 @@ namespace AMDiS {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Parallel::MeshDistributor::globalMeshDistributor->exitParallelization(); Parallel::MeshDistributor::globalMeshDistributor->exitParallelization();
delete Parallel::MeshDistributor::globalMeshDistributor; delete Parallel::MeshDistributor::globalMeshDistributor;
#ifdef HAVE_PARALLEL_MTL4
if (mtl_environment)
delete mtl_environment;
#endif
#endif #endif
#if (defined HAVE_SEQ_PETSC) || (defined HAVE_PARALLEL_PETSC) #ifdef HAVE_PARALLEL_MTL4
if (mtl_environment)
delete mtl_environment;
#elif defined HAVE_PETSC || defined HAVE_SEQ_PETSC || defined HAVE_PARALLEL_PETSC
PetscFinalize(); PetscFinalize();
#elif defined MTL_HAS_HYPRE
MPI_Finalize();
#endif #endif
} }
......
...@@ -137,6 +137,10 @@ namespace AMDiS { ...@@ -137,6 +137,10 @@ namespace AMDiS {
creator = new DiagonalPreconditioner::Creator; creator = new DiagonalPreconditioner::Creator;
addCreator("diag", creator); addCreator("diag", creator);
creator = new MassLumpingPreconditioner::Creator;
addCreator("lumping", creator);
addCreator("masslumping", creator);
creator = new ILUPreconditioner::Creator; creator = new ILUPreconditioner::Creator;
addCreator("ilu", creator); addCreator("ilu", creator);
......
...@@ -575,7 +575,7 @@ namespace AMDiS { ...@@ -575,7 +575,7 @@ namespace AMDiS {
void print() const; void print() const;
/// ///
int calcMemoryUsage() const; size_t calcMemoryUsage() const;
/// Computes the coefficients of the interpolant of the function fct and /// Computes the coefficients of the interpolant of the function fct and
/// stores these in the DOFVector /// stores these in the DOFVector
......
...@@ -402,9 +402,9 @@ namespace AMDiS { ...@@ -402,9 +402,9 @@ namespace AMDiS {
template<typename T> template<typename T>
int DOFVector<T>::calcMemoryUsage() const size_t DOFVector<T>::calcMemoryUsage() const
{ {
int result = 0; size_t result = 0;
result += sizeof(DOFVector<T>); result += sizeof(DOFVector<T>);
result += vec.size() * sizeof(T); result += vec.size() * sizeof(T);
......
...@@ -43,55 +43,142 @@ namespace AMDiS { ...@@ -43,55 +43,142 @@ namespace AMDiS {
typedef MTLTypes::MTLVector VectorType; typedef MTLTypes::MTLVector VectorType;
typedef MTL4Runner< MatrixType, VectorType > super; typedef MTL4Runner< MatrixType, VectorType > super;
/** Interface to the HYPRE BoomerAMG solver [...]
* Parameters provided by AMDiS:
*
* [solver]->cycle mode:
* 1...V-cycle
* 2...W-cycle
*
* [solver]->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
*
* [solver]->info:
* 0...no printout (default)
* 1...print setup information
* 2...print solve information
* 3...print both setup and solve information
*
* [solver]->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)
* */
Hypre_Runner(LinearSolver* oemPtr) Hypre_Runner(LinearSolver* oemPtr)
: oem(*oemPtr), : oem(*oemPtr),
solver(nullptr) useTransposed(false),
solverCreated(false)
{ {
int cycleMode = 1, interpolation = 0; int cycleMode = -1, interpolation = -1, relaxation = -1;
Parameters::get(oem.getName() + "->cycle mode", cycleMode); Parameters::get(oem.getName() + "->cycle mode", cycleMode);
Parameters::get(oem.getName() + "->interpolation type", interpolation); Parameters::get(oem.getName() + "->interpolation type", interpolation);
Parameters::get(oem.getName() + "->relax type", relaxation);
config.maxIter(oem.getMaxIterations()); config.maxIter(oem.getMaxIterations());
config.interpolation(interpolation); config.interpolationType(interpolation);
config.relaxType(relaxation);
config.cycle(cycleMode); config.cycle(cycleMode);
config.tolerance(oem.getTolerance()); config.tolerance(oem.getRelative());
config.printLevel(oem.getInfo()); config.printLevel(oem.getInfo());
} }
~Hypre_Runner() ~Hypre_Runner()
{ {
if (solver) { exit();
delete solver;
solver = nullptr;
}
} }
/// Implementation of \ref MTL4Runner::init() /// Implementation of \ref MTL4Runner::init()
void init(const typename super::BlockMatrix& A, const MatrixType& fullMatrix) override void init(const typename super::BlockMatrix& A, const MatrixType& mtlMatrix) override
{ {
solver = new mtl::HypreSolverWrapper(fullMatrix); setTransposed(typename MatrixType::orientation());
// TODO: copy matrix directly from DOFMatrix to HYPRE matrix (?)
hypreMatrix.init(mtlMatrix);
HYPRE_IJMatrixGetObject(hypreMatrix, (void**) &matrix);
HYPRE_BoomerAMGCreate(&solver);
mtl::dense_vector<double> swap(1, 0.0);
mtl::HypreParVector x(swap);
HYPRE_BoomerAMGSetup(solver, matrix, x, x);
solverCreated = true;
} }
/// Implementation of \ref MTL4Runner::solve() /// Implementation of \ref MTL4Runner::solve()
int solve(const MatrixType& A , VectorType& x, const VectorType& b) override int solve(const MatrixType& A , VectorType& mtlX, const VectorType& mtlB) override
{ {
int error = (*solver)(x, b, config); mtl::HypreParVector x(mtlX);
mtl::HypreParVector b(mtlB);
config(solver);
int error = 0;
if(useTransposed)
error = HYPRE_BoomerAMGSolveT(solver, matrix, b, x);
else
error = HYPRE_BoomerAMGSolve(solver, matrix, b, x);
mtl::convert(x.getHypreVector(), mtlX);
int num_iter = 0;
HYPRE_BoomerAMGGetNumIterations(solver, &num_iter);
oem.setIterations(num_iter);
double rel_resid = 0.0;
HYPRE_BoomerAMGGetFinalRelativeResidualNorm(solver, &rel_resid);
oem.setRelativeResidual(rel_resid);
oem.setErrorCode(error);
return error; return error;
} }
/// Implementation of \ref OEMRunner::exit() /// Implementation of \ref OEMRunner::exit()
void exit() void exit()
{ {
if (solver) { if (solverCreated)
delete solver; HYPRE_BoomerAMGDestroy(solver);
solver = nullptr; solverCreated = false;
} }
private:
inline void setTransposed(mtl::row_major)
{
useTransposed = false;
}
inline void setTransposed(mtl::col_major)
{
useTransposed = true;
} }
protected: protected:
LinearSolver& oem; LinearSolver& oem;
mtl::HypreSolverWrapper *solver;
private:
HYPRE_Solver solver;
HYPRE_ParCSRMatrix matrix;
mtl::HypreMatrix hypreMatrix;
mtl::AMGConfigurator config; mtl::AMGConfigurator config;
bool useTransposed;
bool solverCreated;
}; };
......
...@@ -29,6 +29,8 @@ ...@@ -29,6 +29,8 @@
#include "DOFMatrix.h" #include "DOFMatrix.h"
#include "CreatorInterface.h" #include "CreatorInterface.h"
#include "itl/masslumping.hpp"
#include <boost/numeric/itl/itl.hpp> #include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/itl/pc/ilu_0.hpp> #include <boost/numeric/itl/pc/ilu_0.hpp>
#include <boost/numeric/itl/pc/ic_0.hpp> #include <boost/numeric/itl/pc/ic_0.hpp>
...@@ -153,6 +155,15 @@ namespace AMDiS { ...@@ -153,6 +155,15 @@ namespace AMDiS {
* Diagonal preconditioner \f$ M^{-1} \f$ for the system \f$ Ax=b \f$ is defined as: \f$ M=diag(A) \f$. * Diagonal preconditioner \f$ M^{-1} \f$ for the system \f$ Ax=b \f$ is defined as: \f$ M=diag(A) \f$.
*/ */
typedef ITL_Preconditioner<itl::pc::diagonal<MTLTypes::MTLMatrix>, MTLTypes::MTLMatrix, MTLTypes::MTLVector > DiagonalPreconditioner; typedef ITL_Preconditioner<itl::pc::diagonal<MTLTypes::MTLMatrix>, MTLTypes::MTLMatrix, MTLTypes::MTLVector > DiagonalPreconditioner;
/**
* \ingroup Solver
* \class AMDiS::DiagonalPreconditioner
* \brief ITL_Preconditioner implementation of diagonal (jacobi) preconditioner \implements ITL_Preconditioner
*
* Diagonal preconditioner \f$ M^{-1} \f$ for the system \f$ Ax=b \f$ is defined as: \f$ M_ii=sum_j(A_ij) \f$.
*/
typedef ITL_Preconditioner<itl::pc::masslumping<MTLTypes::MTLMatrix>, MTLTypes::MTLMatrix, MTLTypes::MTLVector > MassLumpingPreconditioner;
/** /**
......
...@@ -92,6 +92,7 @@ namespace AMDiS { ...@@ -92,6 +92,7 @@ namespace AMDiS {
max_iter(1000), max_iter(1000),
info(0), info(0),
residual(-1.0), residual(-1.0),
rel_residual(-1.0),
print_cycle(100), print_cycle(100),
iterations(-1), iterations(-1),
error(-1), error(-1),
...@@ -204,21 +205,21 @@ namespace AMDiS { ...@@ -204,21 +205,21 @@ namespace AMDiS {
virtual OEMPreconditioner* getLeftPrecon() virtual OEMPreconditioner* getLeftPrecon()
{ {
FUNCNAME("LinearSolver::getLeftPrecon()"); FUNCNAME("LinearSolver::getLeftPrecon()");
ERROR_EXIT("no left preconditioner provided!\n"); ERROR("no left preconditioner provided!\n");
return nullptr; return nullptr;
} }
virtual OEMPreconditioner* getRightPrecon() virtual OEMPreconditioner* getRightPrecon()
{ {
FUNCNAME("LinearSolver::getRightPrecon()"); FUNCNAME("LinearSolver::getRightPrecon()");
ERROR_EXIT("no right preconditioner provided!\n"); ERROR("no right preconditioner provided!\n");
return nullptr; return nullptr;
} }
virtual OEMRunner* getRunner() virtual OEMRunner* getRunner()
{ {
FUNCNAME("LinearSolver::getRunner()"); FUNCNAME("LinearSolver::getRunner()");
ERROR_EXIT("no runner provided!\n"); ERROR("no runner provided!\n");
return nullptr; return nullptr;
} }
...@@ -245,6 +246,11 @@ namespace AMDiS { ...@@ -245,6 +246,11 @@ namespace AMDiS {
relative = rel; relative = rel;
} }
inline void setRelativeResidual(double r)
{
rel_residual = r;
}
/// Sets \ref max_iter /// Sets \ref max_iter
inline void setMaxIterations(int i) inline void setMaxIterations(int i)
{ {
...@@ -294,8 +300,11 @@ namespace AMDiS { ...@@ -294,8 +300,11 @@ namespace AMDiS {
/// info level during solving the system. Set in LinearSolver's constructor. /// info level during solving the system. Set in LinearSolver's constructor.
int info; int info;
/// current residual. /// current absolute residual norm.
double residual; double residual;
/// current relative residual norm.
double rel_residual;
/// Print cycle, after how many iterations the residuum norm is logged. /// Print cycle, after how many iterations the residuum norm is logged.
int print_cycle; int print_cycle;
......
...@@ -52,27 +52,24 @@ int gmres_householder_full(const Matrix &A, Vector &x, const Vector &b, ...@@ -52,27 +52,24 @@ int gmres_householder_full(const Matrix &A, Vector &x, const Vector &b,
if (size(b) == 0) throw mtl::logic_error("empty rhs vector"); if (size(b) == 0) throw mtl::logic_error("empty rhs vector");
const Scalar zero= math::zero(Scalar()), dbl_tol= 1.e-16; // kappa in [1.25, 100] const Scalar zero= math::zero(Scalar()), dbl_tol= 1.e-16;
Scalar rho, bnrm2, temp, alpha, beta; Scalar rho, bnrm2, temp, beta;
Size k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations()))); Size k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
Vector w(b - A *x), r(solve(L,w)); Vector w(b - A *x), r(solve(L,w));
mtl::matrix::multi_vector<Vector> V(Vector(resource(x), zero), kmax+1); mtl::matrix::multi_vector<Vector> V(Vector(resource(x), zero), kmax+1);
mtl::vector::dense_vector<Scalar> sn(kmax, zero), cs(kmax, zero), s(kmax+1, zero), y(kmax, zero); // replicated in distributed solvers mtl::vector::dense_vector<Scalar> sn(kmax, zero), cs(kmax, zero), s(kmax+1, zero), y(kmax, zero); // replicated in distributed solvers
mtl::matrix::dense2D<Scalar> H(kmax, kmax); mtl::matrix::dense2D<Scalar> H(kmax, kmax);
bnrm2 = two_norm(b); bnrm2 = two_norm(b);
if (bnrm2 < dbl_tol) if (bnrm2 < dbl_tol)
bnrm2 = 1.0; bnrm2 = 1.0;
// TODO: if rhs=0 => solution=0 ??
temp = two_norm(r); // norm of preconditioned residual temp = two_norm(r); // norm of preconditioned residual
rho = temp * reciprocal(bnrm2); rho = temp * reciprocal(bnrm2);
if (iter.finished(rho)) // initial guess is good enough solution if (iter.finished(rho)) // initial guess is good enough solution
return iter; return iter;
// u = r + sing(r(0))*||r||*e0 // u = r + sign(r(0))*||r||*e0
beta = signum(r[0])*temp; beta = signum(r[0])*temp;
w = r; w = r;
w[0] += beta; w[0] += beta;
...@@ -114,13 +111,13 @@ int gmres_householder_full(const Matrix &A, Vector &x, const Vector &b, ...@@ -114,13 +111,13 @@ int gmres_householder_full(const Matrix &A, Vector &x, const Vector &b,
irange range_to_end(k+1,imax); irange range_to_end(k+1,imax);
set_to_zero(V.vector(k+1)); set_to_zero(V.vector(k+1));
V.vector(k+1)[range_to_end] = w[range_to_end]; V.vector(k+1)[range_to_end] = w[range_to_end];
alpha = two_norm(V.vector(k+1)); beta = two_norm(V.vector(k+1));
if (alpha != 0.0) { if (beta != 0.0) {
alpha *= signum(w[k+1]); beta *= signum(w[k+1]);
V.vector(k+1)[k+1] += alpha; V.vector(k+1)[k+1] += beta;
V.vector(k+1) *= reciprocal(two_norm(V.vector(k+1))); V.vector(k+1) *= reciprocal(two_norm(V.vector(k+1)));
w[k+1] = -alpha; w[k+1] = -beta;
} }
for (Size i= 0; i < k; i++) { for (Size i= 0; i < k; i++) {
...@@ -133,7 +130,7 @@ int gmres_householder_full(const Matrix &A, Vector &x, const Vector &b, ...@@ -133,7 +130,7 @@ int gmres_householder_full(const Matrix &A, Vector &x, const Vector &b,
s[k+1] = -sn[k]*s[k]; s[k+1] = -sn[k]*s[k];
s[k] = conj(cs[k])*s[k]; s[k] = conj(cs[k])*s[k];
w[k] = cs[k]*w[k] + sn[k]*w[k+1]; w[k] = cs[k]*w[k] + sn[k]*w[k+1];
w[k+1] = 0.0; w[k+1] = 0.0;
irange range(num_rows(H)); irange range(num_rows(H));
......
...@@ -35,8 +35,16 @@ namespace mtl { ...@@ -35,8 +35,16 @@ namespace mtl {
class HypreMatrix class HypreMatrix
{ {
public: public:
HypreMatrix() : initialized(false) {}
template< typename Matrix > template< typename Matrix >
HypreMatrix(const Matrix& matrix) HypreMatrix(const Matrix& matrix)
{
init(matrix);
}
template< typename Matrix >
void init(const Matrix& matrix)
{ {
MPI_Comm comm = MPI_COMM_WORLD; MPI_Comm comm = MPI_COMM_WORLD;
int ilower(0); int ilower(0);
...@@ -69,17 +77,21 @@ namespace mtl { ...@@ -69,17 +77,21 @@ namespace mtl {
delete [] rows; delete [] rows;
delete [] cols; delete [] cols;
initialized = true;
} }
~HypreMatrix() ~HypreMatrix()
{ {
HYPRE_IJMatrixDestroy(ij_matrix); if(initialized)
HYPRE_IJMatrixDestroy(ij_matrix);
initialized = false;
} }
inline operator HYPRE_IJMatrix() { return ij_matrix; } inline operator HYPRE_IJMatrix() { return ij_matrix; }
private: private:
HYPRE_IJMatrix ij_matrix; HYPRE_IJMatrix ij_matrix;
bool initialized;
}; };
class HypreVector class HypreVector
...@@ -145,87 +157,16 @@ namespace mtl { ...@@ -145,87 +157,16 @@ namespace mtl {
HYPRE_ParVector vec; HYPRE_ParVector vec;
HypreVector vec2; HypreVector vec2;
}; };
class HypreSolverWrapper
{
public:
template< typename Matrix >
HypreSolverWrapper(const Matrix& mat)
: matrix(mat)
{
setTransposed(typename Matrix::orientation());
HYPRE_IJMatrixGetObject(matrix, (void**) &A);
HYPRE_BoomerAMGCreate(&solver);
mtl::dense_vector< double > swap(1, 0.0);
HypreParVector x(swap);
HYPRE_BoomerAMGSetup(solver, A, x, x);
}
inline void solve(HypreParVector& x, HypreParVector& b) const
{
if(useTransposed)
HYPRE_BoomerAMGSolveT(solver, A, b, x);