From 59ec2774bc4a7f8f1b46cecf945e631c08e56e7f Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Mon, 6 Jun 2011 10:41:09 +0000
Subject: [PATCH] Word-around for a bug in solving Rosenbrock system using
UMFPACK.
---
AMDiS/src/DOFMatrix.cc | 3 +-
AMDiS/src/Debug.cc | 1 +
AMDiS/src/ITL_OEMSolver.h | 22 ++++----
AMDiS/src/ITL_OEMSolver.hh | 70 +++++++++++++-------------
AMDiS/src/MTL4Solver.h | 68 ++++++++++++-------------
AMDiS/src/OEMSolver.h | 49 ++----------------
AMDiS/src/ProblemStat.cc | 2 +
AMDiS/src/UmfPackSolver.h | 19 ++++---
AMDiS/src/time/RosenbrockStationary.cc | 3 +-
9 files changed, 97 insertions(+), 140 deletions(-)
diff --git a/AMDiS/src/DOFMatrix.cc b/AMDiS/src/DOFMatrix.cc
index 6b19f136..c1cbe8ab 100644
--- a/AMDiS/src/DOFMatrix.cc
+++ b/AMDiS/src/DOFMatrix.cc
@@ -241,7 +241,8 @@ namespace AMDiS {
} else {
for (int j = 0; j < nCol; j++) {
DegreeOfFreedom col = colIndices[j];
- ins[row][col] += elMat[i][j];
+ if (fabs(elMat[i][j]) > 1e-10)
+ ins[row][col] += elMat[i][j];
}
}
}
diff --git a/AMDiS/src/Debug.cc b/AMDiS/src/Debug.cc
index ca4f4ff1..305ae867 100644
--- a/AMDiS/src/Debug.cc
+++ b/AMDiS/src/Debug.cc
@@ -437,6 +437,7 @@ namespace AMDiS {
std::ofstream out;
out.open(filename.c_str());
+ out.precision(20);
for (cursor_type cursor = begin<major>(mat.getBaseMatrix()),
cend = end<major>(mat.getBaseMatrix()); cursor != cend; ++cursor)
diff --git a/AMDiS/src/ITL_OEMSolver.h b/AMDiS/src/ITL_OEMSolver.h
index 6889a75e..853afdbf 100644
--- a/AMDiS/src/ITL_OEMSolver.h
+++ b/AMDiS/src/ITL_OEMSolver.h
@@ -45,11 +45,11 @@ namespace AMDiS {
template <typename ITLSolver>
class ITL_OEMSolver : public MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver > > {
- protected:
- public:
+ 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 > >(name) {}
+ MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, ITL_OEMSolver_runner< ITLSolver > >(name) {}
~ITL_OEMSolver() {}
@@ -77,19 +77,19 @@ namespace AMDiS {
*/
template <typename ITLSolver>
class ITL_OEMSolver_para
- : public MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > > // ITL_OEMSolver<ITLSolver>
+ : public MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > > // ITL_OEMSolver<ITLSolver>
{
public:
- /// The constructor reads needed parameters and sets solvers \ref name.
- ITL_OEMSolver_para(std::string name)
- : MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > >(name)
- {
+ /// The constructor reads needed parameters and sets solvers \ref name.
+ ITL_OEMSolver_para(std::string name)
+ : MTL4Solver< MTLMatrix, MTLVector, ITL_OEMSolver_para_runner< ITLSolver > >(name)
+ {
- }
+ }
- ~ITL_OEMSolver_para() {}
+ ~ITL_OEMSolver_para() {}
- /// Set parameter of iterative solver
+ /// Set parameter of iterative solver
class Creator : public OEMSolverCreator
{
diff --git a/AMDiS/src/ITL_OEMSolver.hh b/AMDiS/src/ITL_OEMSolver.hh
index 248584d3..a515a62e 100644
--- a/AMDiS/src/ITL_OEMSolver.hh
+++ b/AMDiS/src/ITL_OEMSolver.hh
@@ -9,51 +9,51 @@ namespace AMDiS {
/// Pointer to the right preconditioner
ITL_BasePreconditioner< Vector >* r;
PreconPair():
- l(NULL), r(NULL)
- {}
+ l(NULL), r(NULL)
+ {}
};
struct BaseITL_runner {
OEMSolver& oem;
BaseITL_runner(OEMSolver* oemPtr):
- oem(*oemPtr)
+ oem(*oemPtr)
{
TEST_EXIT(oemPtr != NULL)("need a real oemsolver\n");
}
template< typename Vector , typename Matrix >
- void setPrecon(PreconPair< Vector >& pair, const Matrix& matrix)
- {
- /// Creator for the left preconditioner
- CreatorInterface< ITL_BasePreconditioner< Vector > >* leftCreator(NULL);
-
- /// Creator for the right preconditioner
- CreatorInterface< ITL_BasePreconditioner< Vector > >* rightCreator(NULL);
-
- std::string preconType("no");
- Parameters::get(oem.getName() + "->left precon", preconType);
-
- leftCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
- TEST_EXIT(leftCreator != NULL)
- ("there is no creator for the given left preconditioner");
-
- preconType = "no";
- Parameters::get(oem.getName() + "->right precon", preconType);
- rightCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
- TEST_EXIT(rightCreator != NULL)
- ("there is no creator for the given right preconditioner");
-
- pair.l = leftCreator->create(matrix);
- pair.r = rightCreator->create(matrix);
- }
+ void setPrecon(PreconPair< Vector >& pair, const Matrix& matrix)
+ {
+ /// Creator for the left preconditioner
+ CreatorInterface< ITL_BasePreconditioner< Vector > >* leftCreator(NULL);
+
+ /// Creator for the right preconditioner
+ CreatorInterface< ITL_BasePreconditioner< Vector > >* rightCreator(NULL);
+
+ std::string preconType("no");
+ Parameters::get(oem.getName() + "->left precon", preconType);
+
+ leftCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
+ TEST_EXIT(leftCreator != NULL)
+ ("there is no creator for the given left preconditioner");
+
+ preconType = "no";
+ Parameters::get(oem.getName() + "->right precon", preconType);
+ rightCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
+ TEST_EXIT(rightCreator != NULL)
+ ("there is no creator for the given right preconditioner");
+
+ pair.l = leftCreator->create(matrix);
+ pair.r = rightCreator->create(matrix);
+ }
};
template< typename ITLSolver >
struct ITL_OEMSolver_para_runner : BaseITL_runner {
ITL_OEMSolver_para_runner(OEMSolver* oem_):
- BaseITL_runner(oem_),
- ell(1)
+ BaseITL_runner(oem_),
+ ell(1)
{
ell=oem.getMaxIterations();
Parameters::get(oem.getName() + "->ell", ell);
@@ -63,11 +63,11 @@ namespace AMDiS {
int solve(const Matrix& A, Vector& x, Vector& b)
{
itl::cyclic_iteration<typename Matrix::value_type> iter(b, oem.getMaxIterations(), oem.getRelative(),
- oem.getTolerance(), oem.getPrint_cycle());
+ oem.getTolerance(), oem.getPrint_cycle());
static PreconPair< Vector > preconPair;
if (preconPair.l == NULL || preconPair.r == NULL || !oem.getMultipleRhs())
- BaseITL_runner::setPrecon(preconPair, A);
+ BaseITL_runner::setPrecon(preconPair, A);
int error = ITLSolver()(A, x, b, *(preconPair.l), *(preconPair.r), iter, ell);
oem.setIterations(iter.iterations());
@@ -75,7 +75,7 @@ namespace AMDiS {
oem.setErrorCode(error);
return error;
}
- private:
+ private:
int ell;
};
@@ -84,18 +84,18 @@ namespace AMDiS {
ITL_OEMSolver_runner(OEMSolver* oem):
- BaseITL_runner(oem)
+ BaseITL_runner(oem)
{}
template< typename Matrix, typename Vector>
int solve(const Matrix& A, Vector& x, Vector& b)
{
itl::cyclic_iteration<typename Matrix::value_type> iter(b, oem.getMaxIterations(), oem.getRelative(),
- oem.getTolerance(), oem.getPrint_cycle());
+ oem.getTolerance(), oem.getPrint_cycle());
static PreconPair< Vector > preconPair;
if (preconPair.l == NULL || preconPair.r == NULL || !oem.getMultipleRhs())
- BaseITL_runner::setPrecon(preconPair, A);
+ BaseITL_runner::setPrecon(preconPair, A);
int error = ITLSolver()(A, x, b, *(preconPair.l), *(preconPair.r), iter);
oem.setErrorCode(error);
diff --git a/AMDiS/src/MTL4Solver.h b/AMDiS/src/MTL4Solver.h
index b6ba1907..1d26715f 100644
--- a/AMDiS/src/MTL4Solver.h
+++ b/AMDiS/src/MTL4Solver.h
@@ -25,53 +25,49 @@
#include <iostream>
namespace AMDiS {
- template< typename MTLMatrix, typename MTLVector, typename Worker >
- class MTL4Solver : public OEMSolver {
-
- protected:
- Worker worker;
-
- template< typename Matrix, typename Vector, typename Mapper >
- int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper)
- {
- MTLMatrix matrix(mapper.nRow(), mapper.nCol());
- set_to_zero(matrix);
- MatMap< const Matrix, Mapper > matMap(A,mapper);
- matrix << matMap;
+ template< typename MTLMatrix, typename MTLVector, typename Worker >
+ class MTL4Solver : public OEMSolver {
+ protected:
+ Worker worker;
- MTLVector mtl_x(mapper.nRow());
- set_to_zero(mtl_x);
- VecMap< Vector, Mapper > xVecMap(x, mapper);
- mtl_x << xVecMap;
+ template< typename Matrix, typename Vector, typename Mapper >
+ int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper)
+ {
+ MTLMatrix matrix(mapper.nRow(), mapper.nCol());
+ set_to_zero(matrix);
+ MatMap< const Matrix, Mapper > matMap(A,mapper);
+ matrix << matMap;
- MTLVector mtl_b(mapper.nRow());
- set_to_zero(mtl_b);
- VecMap< Vector, Mapper> bVecMap(b, mapper);
- mtl_b << bVecMap;
+ MTLVector mtl_x(mapper.nRow());
+ set_to_zero(mtl_x);
+ VecMap< Vector, Mapper > xVecMap(x, mapper);
+ mtl_x << xVecMap;
- error = worker.solve(matrix, mtl_x, mtl_b);
+ MTLVector mtl_b(mapper.nRow());
+ set_to_zero(mtl_b);
+ VecMap< Vector, Mapper> bVecMap(b, mapper);
+ mtl_b << bVecMap;
- mtl_x >> xVecMap;
- return error;
- }
+ error = worker.solve(matrix, mtl_x, mtl_b);
- public:
- MTL4Solver(std::string n):
- OEMSolver(n),
- worker(this)
- {
- FUNCNAME("MTL4Solver()");
+ mtl_x >> xVecMap;
+ return error;
+ }
- }
+ public:
+ MTL4Solver(std::string n):
+ OEMSolver(n),
+ worker(this)
+ {}
- virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
+ virtual 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);
+ }
};
diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h
index 8d30802b..329bad38 100644
--- a/AMDiS/src/OEMSolver.h
+++ b/AMDiS/src/OEMSolver.h
@@ -73,9 +73,7 @@ namespace AMDiS {
print_cycle(100),
iterations(-1),
error(-1),
- multipleRhs(false)/*,
- leftPrecon(NULL),
- rightPrecon(NULL)*/
+ multipleRhs(false)
{}
///
@@ -99,47 +97,13 @@ namespace AMDiS {
SystemVector& x,
SystemVector& b,
VectorialMapper&) = 0;
-#if 0
- {
- int ns = x.getSize(), // Number of systems.
- size = x.getUsedSize(); // Size of all DOFVectors
-
- // Copy vectors
- mtl::dense_vector<value_type> xx(size), bb(size);
- xx << x;
- bb << b
-/* for (int i = 0, counter = 0; i < ns; i++) {
- DOFVector<double>::Iterator it_b(b.getDOFVector(i), USED_DOFS);
- DOFVector<double>::Iterator it_x(x.getDOFVector(i), USED_DOFS);
-
- for (it_b.reset(), it_x.reset(); !it_b.end(); ++it_b, ++it_x) {
- bb[counter] = *it_b;
- xx[counter] = *it_x;
- counter++;
- }
- }
- */
- // Solver on DOFVector for single system
- DOFMatrix::base_matrix_type baseA;
- baseA << A;
- int r = solveSystem(baseA, xx, bb);
-
- for (int i = 0, counter = 0; i < ns; i++) {
- DOFVector<double>::Iterator it(x.getDOFVector(i), USED_DOFS);
- for (it.reset(); !it.end(); ++it)
- *it = xx[counter++];
- }
-
- return r;
- }
-#endif
inline int solveSystem(const SolverMatrix< Matrix< DOFMatrix* > >& A,
SystemVector& x,
SystemVector& b)
{
- VectorialMapper mapper(A);
- return solveSystem(A, x, b, mapper);
+ VectorialMapper mapper(A);
+ return solveSystem(A, x, b, mapper);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
@@ -306,13 +270,6 @@ namespace AMDiS {
* as they can do the factorization of the matrix only once.
*/
bool multipleRhs;
-#if 0
- /// Left preconditioner (ignored by some iterative solvers and by UmfPack)
- PreconCreator *leftPrecon;
-
- /// Right preconditioner (ignored by some iterative solvers and by UmfPack)
- PreconCreator *rightPrecon;
-#endif
};
/**
diff --git a/AMDiS/src/ProblemStat.cc b/AMDiS/src/ProblemStat.cc
index fd6dd0cd..288a6e2b 100644
--- a/AMDiS/src/ProblemStat.cc
+++ b/AMDiS/src/ProblemStat.cc
@@ -634,6 +634,8 @@ namespace AMDiS {
{
FUNCNAME("ProblemStat::buildAfterCoarsen()");
+ VtkWriter::writeFile(rhs->getDOFVector(0), "test.vtu");
+
// buildAfterCoarsen_sebastianMode(adaptInfo, flag);
// return;
diff --git a/AMDiS/src/UmfPackSolver.h b/AMDiS/src/UmfPackSolver.h
index 795153ce..783cdb8b 100644
--- a/AMDiS/src/UmfPackSolver.h
+++ b/AMDiS/src/UmfPackSolver.h
@@ -52,6 +52,8 @@ namespace AMDiS {
store_symbolic(0),
symmetric_strategy(0)
{
+ FUNCNAME("Umfpack_runner::Umfpack_runner()");
+
TEST_EXIT_DBG(oem_ != NULL)("Need real OEMSolver\n");
Parameters::get(oem.getName() + "->store symbolic", store_symbolic);
Parameters::get(oem.getName() + "->symmetric strategy", symmetric_strategy);
@@ -60,15 +62,14 @@ namespace AMDiS {
template< typename Vector>
int solve(const Matrix& A, Vector& x, Vector& b)
{
- FUNCNAME("(UmfPackSolver.h)::solve()");
- if (!solver) {
+ FUNCNAME("Umfpack_runner::solve()");
+ if (!solver) {
if (!symmetric_strategy)
solver = new mtl::matrix::umfpack::solver<matrix_type>(A);
else
solver = new mtl::matrix::umfpack::solver<matrix_type>(A, UMFPACK_STRATEGY_SYMMETRIC);
} else {
if (!oem.getMultipleRhs()) {
-
if (store_symbolic)
solver->update_numeric();
else
@@ -77,16 +78,14 @@ namespace AMDiS {
}
int code = (*solver)(x, b);
- if( oem.getInfo() > 0) {
+ if (oem.getInfo() > 0) {
mtl::dense_vector<typename matrix_type::value_type> r(b);
r -= A * x;
double residual = two_norm(r);
oem.setResidual(residual);
oem.setErrorCode(code);
std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n";
- if (residual > oem.getTolerance()) {
- ERROR_EXIT("Tolerance tol=%e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());
- }
+ TEST_EXIT(residual <= oem.getTolerance())("Tolerance tol = %e could not be reached!\n Set tolerance by '->solver->tolerance:' \n", oem.getTolerance());
}
return code;
@@ -94,11 +93,11 @@ namespace AMDiS {
~Umfpack_runner()
{
- if(!solver)
+ if (!solver)
delete solver;
}
- private:
+ private:
mtl::matrix::umfpack::solver<matrix_type> *solver;
int store_symbolic;
@@ -108,7 +107,7 @@ namespace AMDiS {
using namespace MTLTypes;
class UmfPackSolver : public MTL4Solver<MTLMatrix, MTLVector, Umfpack_runner< MTLMatrix > >
{
- protected:
+ protected:
public:
/// Creator class used in the OEMSolverMap.
diff --git a/AMDiS/src/time/RosenbrockStationary.cc b/AMDiS/src/time/RosenbrockStationary.cc
index 7b029b06..87ec8da9 100644
--- a/AMDiS/src/time/RosenbrockStationary.cc
+++ b/AMDiS/src/time/RosenbrockStationary.cc
@@ -14,6 +14,7 @@
#include "ProblemStat.h"
#include "SystemVector.h"
#include "OEMSolver.h"
+#include "Debug.h"
namespace AMDiS {
@@ -79,7 +80,7 @@ namespace AMDiS {
}
ProblemStatSeq::buildAfterCoarsen(adaptInfo, flag, (i == 0), asmVector);
- solver->setMultipleRhs(i != 0);
+ // solver->setMultipleRhs(i != 0);
ProblemStatSeq::solve(adaptInfo);
*(stageSolutions[i]) = *solution;
--
GitLab