diff --git a/AMDiS/src/ITL_OEMSolver.h b/AMDiS/src/ITL_OEMSolver.h
index 8e0d6b389ac2fd2e91a45cb54f01c87c8358e340..ed7945349dfbc5f1047b34fb9966dd011c99afd1 100644
--- a/AMDiS/src/ITL_OEMSolver.h
+++ b/AMDiS/src/ITL_OEMSolver.h
@@ -44,19 +44,19 @@ namespace AMDiS {
template <typename ITLSolver>
class ITL_OEMSolver : public MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector> > {
-
+
protected:
public:
/// The constructor reads needed parameters and sets solvers \ref name.
ITL_OEMSolver(std::string name) :
MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >(name) {}
-
- int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+
+ int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
SystemVector& b,
VectorialMapper& m)
{
- return MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >::solve(A,x,b,m);
+ return MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >::solve(A,x,b,m);
}
@@ -74,7 +74,7 @@ namespace AMDiS {
return new ITL_OEMSolver<ITLSolver>(this->name);
}
};
-
+
};
/**
@@ -91,22 +91,20 @@ namespace AMDiS {
/// The constructor reads needed parameters and sets solvers \ref name.
ITL_OEMSolver_para(std::string name)
: MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >(name)
- {
-
- }
+ {}
int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
SystemVector& b,
VectorialMapper& m)
{
- return MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >::solve(A,x,b,m);
+ return MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver, MTLTypes::MTLMatrix, MTLTypes::MTLVector > >::solve(A,x,b,m);
}
-
+
~ITL_OEMSolver_para() {}
-
+
/// Set parameter of iterative solver
-
+
class Creator : public OEMSolverCreator
{
public:
diff --git a/AMDiS/src/ITL_OEMSolver.hh b/AMDiS/src/ITL_OEMSolver.hh
index 045ccdb7439456f8174bb6709ae0e7ea0674d58d..db6e3688eff2518084ebe70a16ee93e71db56f2a 100644
--- a/AMDiS/src/ITL_OEMSolver.hh
+++ b/AMDiS/src/ITL_OEMSolver.hh
@@ -58,13 +58,13 @@ namespace AMDiS {
BaseITL_runner(oem_),
ell(1)
{
- ell=oem.getMaxIterations();
+ ell = oem.getMaxIterations();
Parameters::get(oem.getName() + "->ell", ell);
}
void init(const Matrix& A)
{
- BaseITL_runner::setPrecon(preconPair, A);
+ BaseITL_runner::setPrecon(preconPair, A);
}
int solve(const Matrix& A, Vector& x, Vector& b)
@@ -84,6 +84,19 @@ namespace AMDiS {
return error;
}
+
+ void exit()
+ {
+ if (preconPair.l != NULL) {
+ delete preconPair.l;
+ preconPair.l = NULL;
+ }
+
+ if (preconPair.r != NULL) {
+ delete preconPair.r;
+ preconPair.r = NULL;
+ }
+ }
private:
PreconPair< Vector > preconPair;
ITLSolver solver;
@@ -91,35 +104,47 @@ namespace AMDiS {
};
template< typename ITLSolver, typename Matrix, typename Vector >
- struct ITL_OEMSolver_runner : BaseITL_runner {
-
-
+ struct ITL_OEMSolver_runner : BaseITL_runner {
ITL_OEMSolver_runner(OEMSolver* oem):
BaseITL_runner(oem)
{}
-
+
void init(const Matrix& A)
{
- BaseITL_runner::setPrecon(preconPair, A);
+ BaseITL_runner::setPrecon(preconPair, A);
}
+
int solve(const Matrix& A , Vector& x, Vector& b)
{
FUNCNAME("ITL_OEMSolver_runner::solve()");
TEST_EXIT(preconPair.l != NULL)("there is no left preconditioner");
TEST_EXIT(preconPair.r != NULL)("there is no right preconditioner");
-
+
Vector r(A * x - b);
itl::cyclic_iteration<typename Matrix::value_type> iter(r, oem.getMaxIterations(), oem.getRelative(),
oem.getTolerance(), oem.getPrint_cycle());
-
+
int error = solver(A, x, b, *(preconPair.l), *(preconPair.r), iter);
oem.setErrorCode(error);
oem.setIterations(iter.iterations());
oem.setResidual(iter.resid());
-
+
return error;
}
- private:
+
+ void exit()
+ {
+ if (preconPair.l != NULL) {
+ delete preconPair.l;
+ preconPair.l = NULL;
+ }
+
+ if (preconPair.r != NULL) {
+ delete preconPair.r;
+ preconPair.r = NULL;
+ }
+ }
+ private:
ITLSolver solver;
PreconPair< Vector > preconPair;
};
diff --git a/AMDiS/src/ITL_Preconditioner.h b/AMDiS/src/ITL_Preconditioner.h
index cde50924135b299530a3eccca8275368070072c5..12a094eafe074b8a60c3c8095b4910cdae73d59b 100644
--- a/AMDiS/src/ITL_Preconditioner.h
+++ b/AMDiS/src/ITL_Preconditioner.h
@@ -34,10 +34,6 @@
namespace AMDiS {
- // ============================================================================
- // ===== class ITL_Preconditioner ==============================================
- // ============================================================================
-
/**
* \ingroup Solver
@@ -48,26 +44,25 @@ namespace AMDiS {
class ITL_BasePreconditioner
{
public:
- virtual ~ITL_BasePreconditioner() {}
-
- // Old style AMDiS
- // virtual Vec member_solve(const Vec& vin) const =0 ;
- // virtual Vec member_adjoint_solve(const Vec& vin) const =0 ;
-
- // New style without copying the result vector
- virtual void solve(const Vec& x, Vec& y) const= 0;
- virtual void adjoint_solve(const Vec& x, Vec& y) const= 0;
+ virtual ~ITL_BasePreconditioner() {}
+
+ // Old style AMDiS
+ // virtual Vec member_solve(const Vec& vin) const =0 ;
+ // virtual Vec member_adjoint_solve(const Vec& vin) const =0 ;
+
+ // New style without copying the result vector
+ virtual void solve(const Vec& x, Vec& y) const= 0;
+ virtual void adjoint_solve(const Vec& x, Vec& y) const= 0;
};
-
-
-
+
+
+
template< typename Vec >
inline //Vec
itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, false>
solve(const ITL_BasePreconditioner< Vec >& P, const Vec& vin)
{
return itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, false>(P, vin);
- //return P.member_solve(vin);
}
template< typename Vec >
@@ -76,7 +71,6 @@ namespace AMDiS {
adjoint_solve(const ITL_BasePreconditioner< Vec >& P, const Vec& vin)
{
return itl::pc::solver<ITL_BasePreconditioner< Vec >, Vec, true>(P, vin);
- //return P.member_adjoint_solve(vin);
}
/**
@@ -89,31 +83,16 @@ namespace AMDiS {
{
public:
ITL_Preconditioner(const MTLMatrix& A) : precon(A) {}
-#if 0
- MTLVector
- member_solve(const MTLVector& vin) const
- {
- return solve(precon, vin);
- }
-
- MTLVector
- member_adjoint_solve(const MTLVector& vin) const
- {
- return adjoint_solve(precon, vin);
- }
-#else
-
void solve(const MTLVector& vin, MTLVector& vout) const
{
precon.solve(vin, vout);
- }
-
+ }
+
void adjoint_solve(const MTLVector& vin, MTLVector& vout) const
{
precon.adjoint_solve(vin, vout);
- }
-#endif
-
+ }
+
/// Creator class used in the OEMSolverMap.
class Creator : public CreatorInterface< ITL_BasePreconditioner< MTLVector > >
{
@@ -127,16 +106,13 @@ namespace AMDiS {
return new ITL_Preconditioner<Preconditioner, MTLVector, MTLMatrix>(A);
}
};
-
+
private:
- Preconditioner precon;
+ Preconditioner precon;
};
+
- // ============================================================================
- // ===== class DiagonalPreconditioner =========================================
- // ============================================================================
-
/**
* \ingroup Solver
*
@@ -168,26 +144,22 @@ namespace AMDiS {
};
- // ============================================================================
- // ===== class ILUPreconditioner ==============================================
- // ============================================================================
-
/**
* \ingroup Solver
*
* \brief
* ILU (Incomplete LU factorization) preconditioner.
*
- * The preconditioner is used from ITL. It corresponds for instance to "Iterative Methods for
- * Sparce Linear Systems", second edition, Yousef Saad. The preconditioner is
- * described in chapter 10.3 (algorithm 10.4).
+ * The preconditioner is used from ITL. It corresponds for instance to
+ * "Iterative Methods for Sparce Linear Systems", second edition, Yousef Saad.
+ * The preconditioner is described in chapter 10.3 (algorithm 10.4).
*/
class ILUPreconditioner
- : public ITL_Preconditioner< itl::pc::ilu_0<MTLTypes::MTLMatrix>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
+ : public ITL_Preconditioner< itl::pc::ilu_0<MTLTypes::MTLMatrix>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
{};
-
+
class ICPreconditioner
- : public ITL_Preconditioner< itl::pc::ic_0<MTLTypes::MTLMatrix>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
+ : public ITL_Preconditioner< itl::pc::ic_0<MTLTypes::MTLMatrix>, MTLTypes::MTLVector, MTLTypes::MTLMatrix >
{};
diff --git a/AMDiS/src/MTL4Solver.h b/AMDiS/src/MTL4Solver.h
index 1f483fc0191bcfd17c26dde98834b53a638cbc49..e9b88ae88635ccaf2bb18cebe74189c73b1065c3 100644
--- a/AMDiS/src/MTL4Solver.h
+++ b/AMDiS/src/MTL4Solver.h
@@ -28,7 +28,7 @@
#include <boost/numeric/mtl/utility/is_distributed.hpp>
#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(HAVE_PARALLEL_MTL4)
-# include <boost/numeric/mtl/par/distribution.hpp>
+#include <boost/numeric/mtl/par/distribution.hpp>
#endif
namespace AMDiS {
@@ -71,7 +71,7 @@ namespace AMDiS {
FUNCNAME("MTL4Solver::solve()");
Timer t;
- if(num_rows(matrix) == 0 || !getMultipleRhs() ) {
+ if (num_rows(matrix) == 0 || !getMultipleRhs()) {
init(mapper, mtl::traits::is_distributed<MTLMatrix>());
set_to_zero(matrix);
@@ -81,13 +81,13 @@ namespace AMDiS {
}
MTLVector mtl_x;
- init(mtl_x,mtl::traits::is_distributed<MTLVector>());
+ init(mtl_x, mtl::traits::is_distributed<MTLVector>());
set_to_zero(mtl_x);
VecMap< Vector, Mapper > xVecMap(x, mapper);
mtl_x << xVecMap;
MTLVector mtl_b;
- init(mtl_b,mtl::traits::is_distributed<MTLVector>());
+ init(mtl_b, mtl::traits::is_distributed<MTLVector>());
set_to_zero(mtl_b);
VecMap< Vector, Mapper> bVecMap(b, mapper);
mtl_b << bVecMap;
@@ -104,6 +104,9 @@ namespace AMDiS {
MSG("MTL4Solver: ||b-Ax||= %e\n", residual);
mtl_x >> xVecMap;
+
+ worker.exit();
+
return error;
}
diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h
index 0888c1ca3d3fd64492fe4806507b3971dbed2182..c9ce0532b5f21d096ea59b32fa56bad39c19e9ef 100644
--- a/AMDiS/src/OEMSolver.h
+++ b/AMDiS/src/OEMSolver.h
@@ -45,7 +45,6 @@
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
#include "parallel/ParallelMapper.h"
#endif
-//#include "ITL_Preconditioner.h"
namespace AMDiS {
@@ -99,7 +98,8 @@ namespace AMDiS {
VectorialMapper&)
{
FUNCNAME("OEMSolver::solveSystem()");
- TEST_EXIT(false)("This linear solver is not suitable for sequentiell problems\n");
+ TEST_EXIT(false)
+ ("This linear solver is not suitable for sequentiell problems\n");
return -1;
}
@@ -118,7 +118,8 @@ namespace AMDiS {
ParallelMapper&)
{
FUNCNAME("OEMSolver::solveSystem()");
- TEST_EXIT(false)("This linear solver is not suitable for domaindecomposition problems\n");
+ TEST_EXIT(false)
+ ("This linear solver is not suitable for domaindecomposition problems\n");
return -1;
}
#endif
@@ -269,11 +270,9 @@ namespace AMDiS {
/// Error code in last solver (not set by UmfPack)
int error;
- /** \brief
- * If true, the solver has to solve multiple right hand sides with the same
- * matrix. Some solvers, e.g. direct once, may take advanteges from this knowledge,
- * as they can do the factorization of the matrix only once.
- */
+ /// If true, the solver has to solve multiple right hand sides with the same
+ /// matrix. Some solvers, e.g. direct once, may take advanteges from this knowledge,
+ /// as they can do the factorization of the matrix only once.
bool multipleRhs;
};
@@ -295,10 +294,8 @@ namespace AMDiS {
}
protected:
- /** \brief
- * Pointer to the problem the solver will belong to. Needed as parameter
- * when constructing an OEMSolver
- */
+ /// Pointer to the problem the solver will belong to. Needed as parameter
+ /// when constructing an OEMSolver
std::string name;
};
}
diff --git a/AMDiS/src/UmfPackSolver.h b/AMDiS/src/UmfPackSolver.h
index 320b489a9d4bbc1f856b6b8f914dc12031fd2f19..7059879aa0030499ff42d08ac22cf45d7f6272f0 100644
--- a/AMDiS/src/UmfPackSolver.h
+++ b/AMDiS/src/UmfPackSolver.h
@@ -94,9 +94,11 @@ namespace AMDiS {
TEST_EXIT(residual <= oem.getTolerance())("Tolerance tol = %e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());
}
return code;
-
}
+ void exit()
+ {}
+
~Umfpack_runner()
{
if (solver != NULL)
@@ -135,17 +137,15 @@ namespace AMDiS {
/// Constructor
UmfPackSolver(std::string name)
: MTL4Solver< MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >(name)
- {
- }
-
+ {}
+
int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
SystemVector& b,
VectorialMapper& m)
{
- return solve(A,x,b,m);
- }
-
+ return solve(A,x,b,m);
+ }
private:
};