diff --git a/AMDiS/src/AMDiS.cc b/AMDiS/src/AMDiS.cc
index edce09ecd67d25f05c7083af0c2dfa5fe34ee95c..9ea6d73f083da4d98de20226545a14811d5dbb70 100644
--- a/AMDiS/src/AMDiS.cc
+++ b/AMDiS/src/AMDiS.cc
@@ -24,6 +24,9 @@
#ifdef HAVE_ZOLTAN
#include <zoltan_cpp.h>
#endif
+#ifdef MTL_HAS_HYPRE
+#include "mpi.h"
+#endif
#include "boost/program_options.hpp"
namespace AMDiS {
@@ -36,23 +39,24 @@ namespace AMDiS {
void init(int argc, char **argv, std::string initFileName)
{
-#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) || defined(HAVE_SEQ_PETSC)
#ifdef HAVE_PARALLEL_MTL4
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);
- #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();
- #else
+ #else
srand(time(0));
- #endif
- #endif
-
+ #endif
+
#ifdef HAVE_ZOLTAN
float zoltanVersion = 0.0;
Zoltan_Initialize(argc, argv, &zoltanVersion);
#endif
-#endif
Parameters::clearData();
@@ -139,14 +143,15 @@ namespace AMDiS {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Parallel::MeshDistributor::globalMeshDistributor->exitParallelization();
delete Parallel::MeshDistributor::globalMeshDistributor;
- #ifdef HAVE_PARALLEL_MTL4
- if (mtl_environment)
- delete mtl_environment;
- #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();
+#elif defined MTL_HAS_HYPRE
+ MPI_Finalize();
#endif
}
diff --git a/AMDiS/src/CreatorMap.cc b/AMDiS/src/CreatorMap.cc
index 3ce066e4ab0f4b84acc3ad0a9a9fca4d4eb8d2c9..4b5d50d2f81572b5d096911ebaafbf4be13ae35c 100644
--- a/AMDiS/src/CreatorMap.cc
+++ b/AMDiS/src/CreatorMap.cc
@@ -137,6 +137,10 @@ namespace AMDiS {
creator = new DiagonalPreconditioner::Creator;
addCreator("diag", creator);
+ creator = new MassLumpingPreconditioner::Creator;
+ addCreator("lumping", creator);
+ addCreator("masslumping", creator);
+
creator = new ILUPreconditioner::Creator;
addCreator("ilu", creator);
diff --git a/AMDiS/src/DOFVector.h b/AMDiS/src/DOFVector.h
index a1ff631332836be18d49d13e138826834c16f558..e7ec8e03fae647415f22e1e4d2a817b0fcd63a10 100644
--- a/AMDiS/src/DOFVector.h
+++ b/AMDiS/src/DOFVector.h
@@ -575,7 +575,7 @@ namespace AMDiS {
void print() const;
///
- int calcMemoryUsage() const;
+ size_t calcMemoryUsage() const;
/// Computes the coefficients of the interpolant of the function fct and
/// stores these in the DOFVector
diff --git a/AMDiS/src/DOFVector.hh b/AMDiS/src/DOFVector.hh
index a2c80853a0258d6f75bf9cfce479979ed5af8297..fae3568595f8a5395abb31b167d6446e59ee5848 100644
--- a/AMDiS/src/DOFVector.hh
+++ b/AMDiS/src/DOFVector.hh
@@ -402,9 +402,9 @@ namespace AMDiS {
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 += vec.size() * sizeof(T);
diff --git a/AMDiS/src/solver/HypreSolver.h b/AMDiS/src/solver/HypreSolver.h
index fe1a6e8558cc7b6abb86c1cd4b8c4c3efc94ddac..ffe32a7fe552069d878eb90394fd5337c7811361 100644
--- a/AMDiS/src/solver/HypreSolver.h
+++ b/AMDiS/src/solver/HypreSolver.h
@@ -43,55 +43,142 @@ namespace AMDiS {
typedef MTLTypes::MTLVector VectorType;
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)
: 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() + "->interpolation type", interpolation);
-
+ Parameters::get(oem.getName() + "->relax type", relaxation);
+
config.maxIter(oem.getMaxIterations());
- config.interpolation(interpolation);
+ config.interpolationType(interpolation);
+ config.relaxType(relaxation);
config.cycle(cycleMode);
- config.tolerance(oem.getTolerance());
+ config.tolerance(oem.getRelative());
config.printLevel(oem.getInfo());
}
~Hypre_Runner()
{
- if (solver) {
- delete solver;
- solver = nullptr;
- }
+ exit();
}
/// 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()
- 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;
}
/// Implementation of \ref OEMRunner::exit()
void exit()
{
- if (solver) {
- delete solver;
- solver = nullptr;
- }
+ if (solverCreated)
+ HYPRE_BoomerAMGDestroy(solver);
+ solverCreated = false;
+ }
+
+ private:
+ inline void setTransposed(mtl::row_major)
+ {
+ useTransposed = false;
+ }
+
+ inline void setTransposed(mtl::col_major)
+ {
+ useTransposed = true;
}
protected:
LinearSolver& oem;
- mtl::HypreSolverWrapper *solver;
+
+ private:
+ HYPRE_Solver solver;
+ HYPRE_ParCSRMatrix matrix;
+ mtl::HypreMatrix hypreMatrix;
+
mtl::AMGConfigurator config;
+
+ bool useTransposed;
+ bool solverCreated;
};
diff --git a/AMDiS/src/solver/ITL_Preconditioner.h b/AMDiS/src/solver/ITL_Preconditioner.h
index 86ac020a32190cd88366f15db2dc787121c2b539..ca048c12a116d27bb2989259bf039591e7df487a 100644
--- a/AMDiS/src/solver/ITL_Preconditioner.h
+++ b/AMDiS/src/solver/ITL_Preconditioner.h
@@ -29,6 +29,8 @@
#include "DOFMatrix.h"
#include "CreatorInterface.h"
+#include "itl/masslumping.hpp"
+
#include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/itl/pc/ilu_0.hpp>
#include <boost/numeric/itl/pc/ic_0.hpp>
@@ -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$.
*/
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;
/**
diff --git a/AMDiS/src/solver/LinearSolver.h b/AMDiS/src/solver/LinearSolver.h
index 416a07c9dd40b276283e70fa4130c7de1632891c..a5c02743c1291b0cc3e692cb5ca92ef3e7577c3d 100644
--- a/AMDiS/src/solver/LinearSolver.h
+++ b/AMDiS/src/solver/LinearSolver.h
@@ -92,6 +92,7 @@ namespace AMDiS {
max_iter(1000),
info(0),
residual(-1.0),
+ rel_residual(-1.0),
print_cycle(100),
iterations(-1),
error(-1),
@@ -204,21 +205,21 @@ namespace AMDiS {
virtual OEMPreconditioner* getLeftPrecon()
{
FUNCNAME("LinearSolver::getLeftPrecon()");
- ERROR_EXIT("no left preconditioner provided!\n");
+ ERROR("no left preconditioner provided!\n");
return nullptr;
}
virtual OEMPreconditioner* getRightPrecon()
{
FUNCNAME("LinearSolver::getRightPrecon()");
- ERROR_EXIT("no right preconditioner provided!\n");
+ ERROR("no right preconditioner provided!\n");
return nullptr;
}
virtual OEMRunner* getRunner()
{
FUNCNAME("LinearSolver::getRunner()");
- ERROR_EXIT("no runner provided!\n");
+ ERROR("no runner provided!\n");
return nullptr;
}
@@ -245,6 +246,11 @@ namespace AMDiS {
relative = rel;
}
+ inline void setRelativeResidual(double r)
+ {
+ rel_residual = r;
+ }
+
/// Sets \ref max_iter
inline void setMaxIterations(int i)
{
@@ -294,8 +300,11 @@ namespace AMDiS {
/// info level during solving the system. Set in LinearSolver's constructor.
int info;
- /// current residual.
+ /// current absolute residual norm.
double residual;
+
+ /// current relative residual norm.
+ double rel_residual;
/// Print cycle, after how many iterations the residuum norm is logged.
int print_cycle;
diff --git a/AMDiS/src/solver/itl/gmres_householder.hpp b/AMDiS/src/solver/itl/gmres_householder.hpp
index 749e5f38b4a097c7515d9274e3c804627454eec4..a0498257eca4db9df5f3e61c0378bb2406e777d0 100644
--- a/AMDiS/src/solver/itl/gmres_householder.hpp
+++ b/AMDiS/src/solver/itl/gmres_householder.hpp
@@ -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");
- const Scalar zero= math::zero(Scalar()), dbl_tol= 1.e-16; // kappa in [1.25, 100]
- Scalar rho, bnrm2, temp, alpha, beta;
+ const Scalar zero= math::zero(Scalar()), dbl_tol= 1.e-16;
+ Scalar rho, bnrm2, temp, beta;
Size k, kmax(std::min(size(x), Size(iter.max_iterations() - iter.iterations())));
Vector w(b - A *x), r(solve(L,w));
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::matrix::dense2D<Scalar> H(kmax, kmax);
-
-
bnrm2 = two_norm(b);
if (bnrm2 < dbl_tol)
bnrm2 = 1.0;
- // TODO: if rhs=0 => solution=0 ??
temp = two_norm(r); // norm of preconditioned residual
rho = temp * reciprocal(bnrm2);
if (iter.finished(rho)) // initial guess is good enough solution
return iter;
- // u = r + sing(r(0))*||r||*e0
+ // u = r + sign(r(0))*||r||*e0
beta = signum(r[0])*temp;
w = r;
w[0] += beta;
@@ -114,13 +111,13 @@ int gmres_householder_full(const Matrix &A, Vector &x, const Vector &b,
irange range_to_end(k+1,imax);
set_to_zero(V.vector(k+1));
V.vector(k+1)[range_to_end] = w[range_to_end];
- alpha = two_norm(V.vector(k+1));
- if (alpha != 0.0) {
- alpha *= signum(w[k+1]);
- V.vector(k+1)[k+1] += alpha;
+ beta = two_norm(V.vector(k+1));
+ if (beta != 0.0) {
+ beta *= signum(w[k+1]);
+ V.vector(k+1)[k+1] += beta;
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++) {
@@ -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] = 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;
irange range(num_rows(H));
diff --git a/AMDiS/src/solver/itl/hypre.hpp b/AMDiS/src/solver/itl/hypre.hpp
index 9632ce7f1acafb04dcda0051a4d5cd5dd87c26bf..c377c35bdc858ab752c60e1fd8bc19390e87c81d 100644
--- a/AMDiS/src/solver/itl/hypre.hpp
+++ b/AMDiS/src/solver/itl/hypre.hpp
@@ -35,8 +35,16 @@ namespace mtl {
class HypreMatrix
{
public:
+ HypreMatrix() : initialized(false) {}
+
template< typename Matrix >
HypreMatrix(const Matrix& matrix)
+ {
+ init(matrix);
+ }
+
+ template< typename Matrix >
+ void init(const Matrix& matrix)
{
MPI_Comm comm = MPI_COMM_WORLD;
int ilower(0);
@@ -69,17 +77,21 @@ namespace mtl {
delete [] rows;
delete [] cols;
+ initialized = true;
}
~HypreMatrix()
{
- HYPRE_IJMatrixDestroy(ij_matrix);
+ if(initialized)
+ HYPRE_IJMatrixDestroy(ij_matrix);
+ initialized = false;
}
inline operator HYPRE_IJMatrix() { return ij_matrix; }
private:
HYPRE_IJMatrix ij_matrix;
+ bool initialized;
};
class HypreVector
@@ -145,87 +157,16 @@ namespace mtl {
HYPRE_ParVector vec;
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);
- else
- HYPRE_BoomerAMGSolve(solver, A, b, x);
- }
-
- template< typename Vector >
- int operator()(Vector& mtlX, const Vector& mtlB) const
- {
- HypreParVector x(mtlX);
- HypreParVector b(mtlB);
- solve(x, b);
- convert(x.getHypreVector(), mtlX);
- return 0;
- }
-
- template< typename Vector, typename Configurator >
- int operator()(Vector& mtlX, const Vector& mtlB, const Configurator& config)
- {
- HypreParVector x(mtlX);
- HypreParVector b(mtlB);
- config(solver);
- solve(x, b);
- convert(x.getHypreVector(), mtlX);
- return 0;
- }
-
- ~HypreSolverWrapper()
- {
- HYPRE_BoomerAMGDestroy(solver);
- }
-
-
- private:
- inline void setTransposed(row_major )
- {
- useTransposed = false;
- }
-
- inline void setTransposed(col_major )
- {
- useTransposed = true;
- }
-
- private:
- HYPRE_Solver solver;
- HYPRE_ParCSRMatrix A;
- HypreMatrix matrix;
- bool useTransposed;
-
- };
-
class AMGConfigurator
{
public:
- enum CycleMode { VCYCLE=1, WCYCLE=2 };
enum CycleDirection { DOWN=1, UP=2, COARSEST=3 };
AMGConfigurator()
: _debugFlag(0), _printLevel(0), _logLevel(0),
- _tolerance(1.0e-7), _cycle(VCYCLE), _finestSweep(1),
- _interpolation(0), _maxIter(20)
+ _tolerance(1.0e-7), _cycle(-1), _finestSweep(1),
+ _interpolation(-1), _relaxation(-1), _maxIter(20)
{
for(int i(2); i>= 0; --i)
_sweeps[i] = 1;
@@ -235,8 +176,7 @@ namespace mtl {
AMGConfigurator& printLevel(int pl) { _printLevel = pl; return *this; }
AMGConfigurator& logLevel(int ll) { _logLevel = ll; if (ll > 1) throw std::runtime_error("logging > 1 not supported\n2"); return *this; }
AMGConfigurator& tolerance(double tol) { _tolerance = tol; return *this; }
- AMGConfigurator& cycle(CycleMode c) { _cycle = c; return *this; }
- AMGConfigurator& cycle(int c) { if (c == 1) _cycle = VCYCLE; else _cycle = WCYCLE; return *this; }
+ AMGConfigurator& cycle(int c) { _cycle = c; return *this; }
AMGConfigurator& finestSweep(int nr) { _finestSweep = nr; return *this; }
AMGConfigurator& sweep(CycleDirection dir, int nr)
@@ -245,7 +185,8 @@ namespace mtl {
return *this;
}
- AMGConfigurator& interpolation(int type) {_interpolation = type; return *this; }
+ AMGConfigurator& interpolationType(int type) {_interpolation = type; return *this; }
+ AMGConfigurator& relaxType(int type) {_relaxation = type; return *this; }
AMGConfigurator& maxIter(int iter) { _maxIter = iter; return *this;}
inline void operator()(HYPRE_Solver& solver) const
@@ -254,12 +195,17 @@ namespace mtl {
HYPRE_BoomerAMGSetPrintLevel(solver, _printLevel);
HYPRE_BoomerAMGSetLogging(solver, _logLevel);
HYPRE_BoomerAMGSetTol(solver, _tolerance);
- HYPRE_BoomerAMGSetCycleType(solver, _cycle);
+ if (_cycle >= 0)
+ HYPRE_BoomerAMGSetCycleType(solver, _cycle);
HYPRE_BoomerAMGSetNumSweeps(solver, _finestSweep);
for(int i(3); i >= 1; --i) {
HYPRE_BoomerAMGSetCycleNumSweeps(solver, _sweeps[i-1], i);
}
- HYPRE_BoomerAMGSetInterpType(solver, _interpolation);
+ if (_interpolation >= 0)
+ HYPRE_BoomerAMGSetInterpType(solver, _interpolation);
+ if (_relaxation >= 0)
+ HYPRE_BoomerAMGSetRelaxType(solver, _relaxation);
+
HYPRE_BoomerAMGSetMaxIter(solver, _maxIter);
}
@@ -268,10 +214,11 @@ namespace mtl {
int _printLevel;
int _logLevel;
double _tolerance;
- CycleMode _cycle;
+ int _cycle;
int _finestSweep;
int _sweeps[3];
int _interpolation;
+ int _relaxation;
int _maxIter;
};