From f0f8d33cc520f6cea0da971ee3320d2e5bfd346a Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Thu, 2 Oct 2008 11:12:44 +0000
Subject: [PATCH] * Several changes to make Rosenbrock method possible to solve
systems with equal matrix but different right hand side vectors
---
AMDiS/src/BiCGSolver.h | 2 +-
AMDiS/src/BiCGSolver.hh | 2 +-
AMDiS/src/BiCGStab.h | 2 +-
AMDiS/src/BiCGStab.hh | 2 +-
AMDiS/src/BiCGStab2.h | 2 +-
AMDiS/src/BiCGStab2.hh | 2 +-
AMDiS/src/BiCGStab_M.h | 2 +-
AMDiS/src/BiCGStab_M.hh | 2 +-
AMDiS/src/CGSolver.h | 2 +-
AMDiS/src/CGSolver.hh | 2 +-
AMDiS/src/Flag.h | 52 +++++++++++++--------------
AMDiS/src/GMResSolver.h | 3 +-
AMDiS/src/GMResSolver.hh | 3 +-
AMDiS/src/GMResSolver2.h | 3 +-
AMDiS/src/GMResSolver2.hh | 4 +--
AMDiS/src/MultiGridWrapper.h | 10 +++---
AMDiS/src/ODirSolver.h | 3 +-
AMDiS/src/ODirSolver.hh | 2 +-
AMDiS/src/OEMSolver.h | 39 ++++++++++----------
AMDiS/src/OResSolver.h | 3 +-
AMDiS/src/OResSolver.hh | 3 +-
AMDiS/src/PardisoSolver.cc | 3 +-
AMDiS/src/PardisoSolver.h | 3 +-
AMDiS/src/ProblemInstat.h | 22 ++++++------
AMDiS/src/ProblemIterationInterface.h | 12 ++++---
AMDiS/src/ProblemScal.cc | 5 +--
AMDiS/src/ProblemScal.h | 12 ++++---
AMDiS/src/ProblemStatBase.h | 7 ++--
AMDiS/src/ProblemVec.cc | 18 ++++++----
AMDiS/src/ProblemVec.h | 9 +++--
AMDiS/src/SolutionDataStorage.h | 3 +-
AMDiS/src/SolutionDataStorage.hh | 10 ++++--
AMDiS/src/StandardProblemIteration.cc | 11 ++++--
AMDiS/src/StandardProblemIteration.h | 4 +--
AMDiS/src/TFQMR.h | 3 +-
AMDiS/src/TFQMR.hh | 5 +--
AMDiS/src/UmfPackSolver.h | 17 +++++++--
AMDiS/src/UmfPackSolver.hh | 32 ++++++++++++-----
AMDiS/src/VecSymSolver.h | 3 +-
AMDiS/src/VecSymSolver.hh | 2 +-
40 files changed, 197 insertions(+), 129 deletions(-)
diff --git a/AMDiS/src/BiCGSolver.h b/AMDiS/src/BiCGSolver.h
index 2094d717..e428ad4a 100644
--- a/AMDiS/src/BiCGSolver.h
+++ b/AMDiS/src/BiCGSolver.h
@@ -77,7 +77,7 @@ namespace AMDiS {
/** \brief
* realisation of OEMSolver::solveSystem
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b, bool reuseMatrix);
/** \brief
* realisation of OEMSolver::init
diff --git a/AMDiS/src/BiCGSolver.hh b/AMDiS/src/BiCGSolver.hh
index b540eb3f..14965297 100644
--- a/AMDiS/src/BiCGSolver.hh
+++ b/AMDiS/src/BiCGSolver.hh
@@ -35,7 +35,7 @@ namespace AMDiS {
template<typename VectorType>
int BiCGSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x, VectorType *b)
+ VectorType *x, VectorType *b, bool reuseMatrix)
{
FUNCNAME("BiCGSolver::solveSystem");
VectorType *t = h;
diff --git a/AMDiS/src/BiCGStab.h b/AMDiS/src/BiCGStab.h
index fc8e6ce6..e965c80d 100644
--- a/AMDiS/src/BiCGStab.h
+++ b/AMDiS/src/BiCGStab.h
@@ -80,7 +80,7 @@ namespace AMDiS
* realisation of OEMSolver::solveSystem
*/
int solveSystem(MatVecMultiplier<VectorType> *mv,
- VectorType *x, VectorType *b);
+ VectorType *x, VectorType *b, bool reuseMatrix);
/** \brief
* realisation of OEMSolver::init
diff --git a/AMDiS/src/BiCGStab.hh b/AMDiS/src/BiCGStab.hh
index 1058576d..7479cfdf 100644
--- a/AMDiS/src/BiCGStab.hh
+++ b/AMDiS/src/BiCGStab.hh
@@ -42,7 +42,7 @@ namespace AMDiS
template<typename VectorType>
int BiCGStab<VectorType>::solveSystem(MatVecMultiplier<VectorType> *mv,
- VectorType *x, VectorType *b)
+ VectorType *x, VectorType *b, bool reuseMatrix)
{
FUNCNAME("BiCGStab::solveSystem");
diff --git a/AMDiS/src/BiCGStab2.h b/AMDiS/src/BiCGStab2.h
index 8767dd96..054a5bfc 100644
--- a/AMDiS/src/BiCGStab2.h
+++ b/AMDiS/src/BiCGStab2.h
@@ -79,7 +79,7 @@ namespace AMDiS
* realisation of OEMSolver::solveSystem
*/
int solveSystem(MatVecMultiplier<VectorType> *mv,
- VectorType *x, VectorType *b);
+ VectorType *x, VectorType *b, bool reuseMatrix);
/** \brief
* realisation of OEMSolver::init
diff --git a/AMDiS/src/BiCGStab2.hh b/AMDiS/src/BiCGStab2.hh
index 139fff73..e0bdedbc 100644
--- a/AMDiS/src/BiCGStab2.hh
+++ b/AMDiS/src/BiCGStab2.hh
@@ -42,7 +42,7 @@ namespace AMDiS
template<typename VectorType>
int BiCGStab2<VectorType>::solveSystem(MatVecMultiplier<VectorType> *mv,
- VectorType *x, VectorType *b)
+ VectorType *x, VectorType *b, bool reuseMatrix)
{
FUNCNAME("BiCGStab2::solveSystem()");
diff --git a/AMDiS/src/BiCGStab_M.h b/AMDiS/src/BiCGStab_M.h
index 1cada652..c9f36d0f 100644
--- a/AMDiS/src/BiCGStab_M.h
+++ b/AMDiS/src/BiCGStab_M.h
@@ -84,7 +84,7 @@ namespace AMDiS {
/** \brief
* Implements OEMSolver<VectorType>::solve().
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b, bool reuseMatrix);
private:
// Pointers to memory needed for solveSystem
diff --git a/AMDiS/src/BiCGStab_M.hh b/AMDiS/src/BiCGStab_M.hh
index aee65d69..57789557 100644
--- a/AMDiS/src/BiCGStab_M.hh
+++ b/AMDiS/src/BiCGStab_M.hh
@@ -34,7 +34,7 @@ namespace AMDiS {
template<typename VectorType>
int BiCGStab_M<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x, VectorType *b)
+ VectorType *x, VectorType *b, bool reuseMatrix)
{
FUNCNAME("BiCGStab_M::solveSystem");
double old_res = -1.0;
diff --git a/AMDiS/src/CGSolver.h b/AMDiS/src/CGSolver.h
index 81ece304..f39a9690 100644
--- a/AMDiS/src/CGSolver.h
+++ b/AMDiS/src/CGSolver.h
@@ -102,7 +102,7 @@ namespace AMDiS {
/** \brief
* Implements OEMSolver<VectorType>::solve().
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b, bool reuseMatrix);
private:
VectorType *p, *r, *v, *z;
diff --git a/AMDiS/src/CGSolver.hh b/AMDiS/src/CGSolver.hh
index ed43b1df..7b1fa1e7 100644
--- a/AMDiS/src/CGSolver.hh
+++ b/AMDiS/src/CGSolver.hh
@@ -16,7 +16,7 @@ namespace AMDiS {
template<typename VectorType>
int CGSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x, VectorType *b)
+ VectorType *x, VectorType *b, bool reuseMatrix)
{
FUNCNAME("CGSolver::solve()");
diff --git a/AMDiS/src/Flag.h b/AMDiS/src/Flag.h
index 40464381..d393a94c 100644
--- a/AMDiS/src/Flag.h
+++ b/AMDiS/src/Flag.h
@@ -41,36 +41,36 @@ namespace AMDiS {
/** \brief
* Constructs a unset Flag
*/
- inline Flag() : flags(0) {};
+ inline Flag() : flags(0) {}
/** \brief
* Constructs a Flag initialized by f
*/
- inline Flag(const unsigned long f) : flags(f) {};
+ inline Flag(const unsigned long f) : flags(f) {}
/** \brief
* Copy constructor
*/
- inline Flag(const Flag& f) : flags(f.flags) {};
+ inline Flag(const Flag& f) : flags(f.flags) {}
/** \brief
* Destructor
*/
- inline ~Flag() {};
+ inline ~Flag() {}
/** \brief
* Compares two Flags
*/
inline bool operator==(const Flag& f) const {
return (flags == f.flags);
- };
+ }
/** \brief
* Compares two Flags
*/
inline bool operator!=(const Flag& f) const {
return !(f == *this);
- };
+ }
/** \brief
* Assignment operator
@@ -79,60 +79,60 @@ namespace AMDiS {
if (this != &f)
flags = f.flags;
return *this;
- };
+ }
/** \brief
* Typecast
*/
inline operator bool() const {
return isAnySet();
- };
+ }
/** \brief
* Set \ref flags
*/
inline void setFlags(const unsigned long f) {
flags = f;
- };
+ }
/** \brief
* Set \ref flags
*/
inline void setFlags(const Flag& f) {
flags = f.flags;
- };
+ }
/** \brief
* Sets \ref flags to \ref flags | f
*/
inline void setFlag(const unsigned long f) {
flags |= f;
- };
+ }
/** \brief
* Sets \ref flags to \ref flags | f.flags
*/
inline void setFlag(const Flag& f) {
flags |= f.flags;
- };
+ }
/** \brief
* Sets \ref flags to \ref flags & ~f
*/
inline void unsetFlag(const unsigned long f) {
flags &= ~f;
- };
+ }
/** \brief
* Sets \ref flags to \ref flags & ~f.flags
*/
inline void unsetFlag(const Flag& f) {
flags &= ~f.flags;
- };
+ }
inline const unsigned long getFlags() const {
return flags;
- };
+ }
/** \brief
* Returns \ref flags | f.flags
@@ -141,7 +141,7 @@ namespace AMDiS {
Flag r(flags);
r.setFlag(f);
return r;
- };
+ }
/** \brief
* Returns \ref flags & ~f.flags
@@ -150,7 +150,7 @@ namespace AMDiS {
Flag r(flags);
r.unsetFlag(f);
return r;
- };
+ }
/** \brief
* Returns \ref flags | f.flags
@@ -159,7 +159,7 @@ namespace AMDiS {
Flag r(flags);
r.setFlag(f);
return r;
- };
+ }
/** \brief
* Returns \ref flags & f.flags
@@ -168,7 +168,7 @@ namespace AMDiS {
Flag r(flags);
r.flags &= f.flags;
return r;
- };
+ }
/** \brief
* Sets \ref flags to \ref flags &= f.flags
@@ -176,7 +176,7 @@ namespace AMDiS {
inline Flag operator&=(const Flag& f) {
flags &= f.flags;
return *this;
- };
+ }
/** \brief
* Returns \ref flags ^ f.flags
@@ -185,7 +185,7 @@ namespace AMDiS {
Flag r(flags);
r.flags ^=f.flags;
return r;
- };
+ }
/** \brief
* Sets \ref flags to \ref flags & f.flags
@@ -195,7 +195,7 @@ namespace AMDiS {
flags |= f.flags;
};
return *this;
- };
+ }
/** \brief
* Returns ~\ref flags
@@ -204,28 +204,28 @@ namespace AMDiS {
Flag r;
r.flags = ~flags;
return r;
- };
+ }
/** \brief
* Checks whether all set bits of f.flags are set in \ref flags too.
*/
inline bool isSet(const Flag& f) const {
return ((flags&f.flags) == f.flags);
- };
+ }
/** \brief
* Returns !\ref isSet(f)
*/
inline bool isUnset(const Flag& f) const {
return !isSet(f);
- };
+ }
/** \brief
* Returns true if \ref flags != 0
*/
inline bool isAnySet() const {
return (flags != 0);
- };
+ }
protected:
/** \brief
diff --git a/AMDiS/src/GMResSolver.h b/AMDiS/src/GMResSolver.h
index 6cbf0e18..9a03ed34 100644
--- a/AMDiS/src/GMResSolver.h
+++ b/AMDiS/src/GMResSolver.h
@@ -86,7 +86,8 @@ namespace AMDiS {
/** \brief
* realisation of OEMSolver::solveSystem
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv,
+ VectorType *x, VectorType *b, bool reuseMatrix);
/** \brief
* realisation of OEMSolver::init
diff --git a/AMDiS/src/GMResSolver.hh b/AMDiS/src/GMResSolver.hh
index e35fc959..cf442a3e 100644
--- a/AMDiS/src/GMResSolver.hh
+++ b/AMDiS/src/GMResSolver.hh
@@ -157,8 +157,7 @@ namespace AMDiS {
template<typename VectorType>
int GMResSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x,
- VectorType *b)
+ VectorType *x, VectorType *b, bool reuseMatrix)
{
FUNCNAME("GMResSolver::solveSystem()");
diff --git a/AMDiS/src/GMResSolver2.h b/AMDiS/src/GMResSolver2.h
index ace80972..6e820643 100644
--- a/AMDiS/src/GMResSolver2.h
+++ b/AMDiS/src/GMResSolver2.h
@@ -85,7 +85,8 @@ namespace AMDiS {
/** \brief
* realisation of OEMSolver::solveSystem
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv,
+ VectorType *x, VectorType *b, bool reuseMatrix);
/** \brief
* realisation of OEMSolver::init
diff --git a/AMDiS/src/GMResSolver2.hh b/AMDiS/src/GMResSolver2.hh
index 4ffee64e..63fb90c2 100644
--- a/AMDiS/src/GMResSolver2.hh
+++ b/AMDiS/src/GMResSolver2.hh
@@ -157,8 +157,8 @@ namespace AMDiS {
template<typename VectorType>
int GMResSolver2<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x,
- VectorType *b)
+ VectorType *x, VectorType *b,
+ bool reuseMatrix)
{
FUNCNAME("GMResSolver2::solveSystem()");
diff --git a/AMDiS/src/MultiGridWrapper.h b/AMDiS/src/MultiGridWrapper.h
index 841d52ac..a0f6dcf1 100644
--- a/AMDiS/src/MultiGridWrapper.h
+++ b/AMDiS/src/MultiGridWrapper.h
@@ -73,11 +73,12 @@ namespace AMDiS {
int solveSystem(MatVecMultiplier<DOFVector<double> > *mv,
DOFVector<double> *x,
- DOFVector<double> *b)
+ DOFVector<double> *b,
+ bool reuseMatrix)
{
DOFMatrix *matrix =
dynamic_cast<StandardMatVec<DOFMatrix, DOFVector<double> >*>(mv)->getMatrix();
- return mgSolver_->solve(*matrix, *x, *b);
+ return mgSolver_->solve(*matrix, *x, *b, reuseMatrix);
};
protected:
@@ -128,11 +129,12 @@ namespace AMDiS {
int solveSystem(MatVecMultiplier<SystemVector> *mv,
SystemVector *x,
- SystemVector *b)
+ SystemVector *b,
+ bool reuseMatrix)
{
Matrix<DOFMatrix*> *matrix =
dynamic_cast<StandardMatVec<Matrix<DOFMatrix*>, SystemVector>*>(mv)->getMatrix();
- return mgSolver_->solve(*matrix, *x, *b);
+ return mgSolver_->solve(*matrix, *x, *b, reuseMatrix);
};
protected:
diff --git a/AMDiS/src/ODirSolver.h b/AMDiS/src/ODirSolver.h
index 9fab35ad..15ba9c97 100644
--- a/AMDiS/src/ODirSolver.h
+++ b/AMDiS/src/ODirSolver.h
@@ -76,7 +76,8 @@ namespace AMDiS {
/** \brief
* realisation of OEMSolver::solveSystem
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv,
+ VectorType *x, VectorType *b, bool reuseMatrix);
/** \brief
* realisation of OEMSolver::init
diff --git a/AMDiS/src/ODirSolver.hh b/AMDiS/src/ODirSolver.hh
index a7142a31..11fcdaa2 100644
--- a/AMDiS/src/ODirSolver.hh
+++ b/AMDiS/src/ODirSolver.hh
@@ -35,7 +35,7 @@ namespace AMDiS {
template<typename VectorType>
int ODirSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x, VectorType *b)
+ VectorType *x, VectorType *b, bool reuseMatrix)
{
FUNCNAME("ODirSolver::solveSystem");
diff --git a/AMDiS/src/OEMSolver.h b/AMDiS/src/OEMSolver.h
index f16053b7..19a158fd 100644
--- a/AMDiS/src/OEMSolver.h
+++ b/AMDiS/src/OEMSolver.h
@@ -96,7 +96,8 @@ namespace AMDiS {
int solve(MatVecMultiplier<VectorType> *mv,
VectorType *x, VectorType *b,
Preconditioner<VectorType> *lPrecon = NULL,
- Preconditioner<VectorType> *rPrecon = NULL)
+ Preconditioner<VectorType> *rPrecon = NULL,
+ bool reuseMatrix = false)
{
FUNCNAME("OEMSolver::solve()");
TEST_EXIT(mv)("no matVec\n");
@@ -113,7 +114,7 @@ namespace AMDiS {
if (rightPrecon)
rightPrecon->init();
- int iter = solveSystem(mv, x, b);
+ int iter = solveSystem(mv, x, b, reuseMatrix);
if (leftPrecon)
leftPrecon->exit();
@@ -124,7 +125,7 @@ namespace AMDiS {
exit();
return iter;
- };
+ }
// ===== getting-methods ======================================================
@@ -137,28 +138,28 @@ namespace AMDiS {
*/
inline const std::string& getName() {
return name;
- };
+ }
/** \brief
* Returns \ref tolerance
*/
inline double getTolerance() {
return tolerance;
- };
+ }
/** \brief
* Returns \ref max_iter
*/
inline int getMaxIterations() {
return max_iter;
- };
+ }
/** \brief
* Returns \ref residual
*/
inline double getResidual() {
return residual;
- };
+ }
/** \} */
@@ -173,37 +174,36 @@ namespace AMDiS {
*/
inline void setTolerance(double tol) {
tolerance = tol;
- };
+ }
/** \brief
* Sets \ref relative
*/
inline void setRelative(bool rel) {
relative = rel;
- };
+ }
/** \brief
* Sets \ref max_iter
*/
inline void setMaxIterations(int i) {
max_iter = i;
- };
+ }
/** \brief
* Sets \ref info
*/
inline void setInfo(int i) {
info = i;
- };
+ }
inline void setVectorCreator(CreatorInterface<VectorType> *creator) {
vectorCreator = creator;
- };
+ }
- inline CreatorInterface<VectorType>* getVectorCreator()
- {
+ inline CreatorInterface<VectorType>* getVectorCreator() {
return vectorCreator;
- };
+ }
/** \} */
@@ -222,8 +222,9 @@ namespace AMDiS {
* To be overloaded by concrete solvers.
*/
virtual int solveSystem(MatVecMultiplier<VectorType> *mv,
- VectorType *x,
- VectorType *b) = 0;
+ VectorType *x,
+ VectorType *b,
+ bool reuseMatrix) = 0;
/** \brief
@@ -294,14 +295,14 @@ namespace AMDiS {
class OEMSolverCreator : public CreatorInterface<OEMSolver<T> >
{
public:
- virtual ~OEMSolverCreator() {};
+ virtual ~OEMSolverCreator() {}
/** \brief
* Sets \ref problem
*/
void setName(std::string name_) {
name = name_;
- };
+ }
protected:
/** \brief
diff --git a/AMDiS/src/OResSolver.h b/AMDiS/src/OResSolver.h
index ddeadce0..8d731b37 100644
--- a/AMDiS/src/OResSolver.h
+++ b/AMDiS/src/OResSolver.h
@@ -76,7 +76,8 @@ namespace AMDiS {
/** \brief
* realisation of OEMSolver::solveSystem
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv,
+ VectorType *x, VectorType *b, bool reuseMatrix);
/** \brief
* realisation of OEMSolver::init
diff --git a/AMDiS/src/OResSolver.hh b/AMDiS/src/OResSolver.hh
index e3fca6d6..74e95865 100644
--- a/AMDiS/src/OResSolver.hh
+++ b/AMDiS/src/OResSolver.hh
@@ -39,7 +39,8 @@ namespace AMDiS {
template<typename VectorType>
int OResSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x, VectorType *b)
+ VectorType *x, VectorType *b,
+ bool reuseMatrix)
{
FUNCNAME("OResSolver::solveSystem");
diff --git a/AMDiS/src/PardisoSolver.cc b/AMDiS/src/PardisoSolver.cc
index f4ee3a8f..bf395892 100644
--- a/AMDiS/src/PardisoSolver.cc
+++ b/AMDiS/src/PardisoSolver.cc
@@ -11,7 +11,8 @@ namespace AMDiS {
template<>
int PardisoSolver<SystemVector>::solveSystem(MatVecMultiplier<SystemVector> *matVec,
- SystemVector *x, SystemVector *b)
+ SystemVector *x, SystemVector *b,
+ bool reuseMatrix)
{
FUNCNAME("PardisoSolver::solveSystem()");
diff --git a/AMDiS/src/PardisoSolver.h b/AMDiS/src/PardisoSolver.h
index 11af655d..e4141546 100644
--- a/AMDiS/src/PardisoSolver.h
+++ b/AMDiS/src/PardisoSolver.h
@@ -102,7 +102,8 @@ namespace AMDiS {
/** \brief
* Implements OEMSolver<VectorType>::solve().
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv,
+ VectorType *x, VectorType *b, bool reuseMatrix);
private:
/** \brief
diff --git a/AMDiS/src/ProblemInstat.h b/AMDiS/src/ProblemInstat.h
index e3bd1d1a..e1621aaa 100644
--- a/AMDiS/src/ProblemInstat.h
+++ b/AMDiS/src/ProblemInstat.h
@@ -51,7 +51,7 @@ namespace AMDiS {
ProblemStatBase *initialProb)
: name(name_),
initialProblem(initialProb ? initialProb : this)
- {};
+ {}
/** \brief
* Destructor.
@@ -68,53 +68,53 @@ namespace AMDiS {
/** \brief
*/
- virtual void solve(AdaptInfo *adaptInfo) {};
+ virtual void solve(AdaptInfo *adaptInfo, bool fixedMatrix) {}
/** \brief
*/
- virtual void estimate(AdaptInfo *adaptInfo) {};
+ virtual void estimate(AdaptInfo *adaptInfo) {}
/** \brief
*/
- virtual void buildBeforeRefine(AdaptInfo *adaptInfo, Flag) {};
+ virtual void buildBeforeRefine(AdaptInfo *adaptInfo, Flag) {}
/** \brief
*/
- virtual void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag) {};
+ virtual void buildBeforeCoarsen(AdaptInfo *adaptInfo, Flag) {}
/** \brief
*/
- virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag) {};
+ virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag, bool, bool) {}
/** \brief
*/
virtual Flag markElements(AdaptInfo *adaptInfo) {
return 0;
- };
+ }
/** \brief
*/
virtual Flag refineMesh(AdaptInfo *adaptInfo) {
return 0;
- };
+ }
/** \brief
*/
virtual Flag coarsenMesh(AdaptInfo *adaptInfo) {
return 0;
- };
+ }
/** \brief
* Implementation of ProblemTimeInterface::closeTimestep().
*/
- virtual void closeTimestep() {};
+ virtual void closeTimestep() {}
/** \brief
* Returns \ref name.
*/
inline const std::string& getName() {
return name;
- };
+ }
/** \brief
* Used by \ref problemInitial
diff --git a/AMDiS/src/ProblemIterationInterface.h b/AMDiS/src/ProblemIterationInterface.h
index d2f6e7e3..742f9be8 100644
--- a/AMDiS/src/ProblemIterationInterface.h
+++ b/AMDiS/src/ProblemIterationInterface.h
@@ -29,11 +29,13 @@ namespace AMDiS {
class AdaptInfo;
class ProblemStatBase;
- const Flag BUILD = 0x01L;
- const Flag ADAPT = 0x02L;
- const Flag SOLVE = 0x04L;
- const Flag ESTIMATE = 0x08L;
- const Flag MARK = 0x10L;
+ const Flag BUILD = 1; // Assemble vectors and matrices
+ const Flag BUILD_RHS = 2; // Assemble rhs vectors only
+ const Flag ADAPT = 4; // Run adaption procedure
+ const Flag SOLVE = 8; // Solve system
+ const Flag SOLVE_RHS = 16; // Solve system, where only rhs vectors have changed
+ const Flag ESTIMATE = 32; // Estimate error
+ const Flag MARK = 64; // Mark elements
const Flag FULL_ITERATION = BUILD | ADAPT | SOLVE | ESTIMATE | MARK;
const Flag NO_ADAPTION = BUILD | SOLVE | ESTIMATE;
diff --git a/AMDiS/src/ProblemScal.cc b/AMDiS/src/ProblemScal.cc
index cf959792..ede6b31c 100644
--- a/AMDiS/src/ProblemScal.cc
+++ b/AMDiS/src/ProblemScal.cc
@@ -192,7 +192,7 @@ namespace AMDiS {
}
}
- void ProblemScal::solve(AdaptInfo *adaptInfo)
+ void ProblemScal::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
{
FUNCNAME("Problem::solve()");
if (!solver_) {
@@ -542,7 +542,8 @@ namespace AMDiS {
}
}
- void ProblemScal::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag)
+ void ProblemScal::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
+ bool assembleMatrix, bool assembleVector)
{
FUNCNAME("ProblemScal::buildAfterCoarsen()");
diff --git a/AMDiS/src/ProblemScal.h b/AMDiS/src/ProblemScal.h
index d60bf76e..fb6e52cc 100644
--- a/AMDiS/src/ProblemScal.h
+++ b/AMDiS/src/ProblemScal.h
@@ -130,7 +130,7 @@ namespace AMDiS {
* Implementation of ProblemStatBase::solve(). Deligates the solving
* of problems system to \ref solver.
*/
- virtual void solve(AdaptInfo *adaptInfo);
+ virtual void solve(AdaptInfo *adaptInfo, bool fixedMatrix);
/** \brief
* Implementation of ProblemStatBase::estimate(). Deligates the estimation
@@ -172,21 +172,23 @@ namespace AMDiS {
* Implementation of ProblemStatBase::buildAfterCoarsen().
* Assembles \ref A and \ref rhs.
*/
- virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag);
+ virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
+ bool assembleMatrix = true,
+ bool assembleVector = true);
/** \brief
* Returns number of managed problems
*/
virtual int getNumProblems() {
return 1;
- };
+ }
/** \brief
* Implementation of ProblemStatBase::getNumComponents()
*/
virtual int getNumComponents() {
return 1;
- };
+ }
/** \brief
* Returns the problem with the given number. If only one problem
@@ -194,7 +196,7 @@ namespace AMDiS {
*/
virtual ProblemStatBase *getProblem(int number = 0) {
return this;
- };
+ }
/** \brief
* Writes output files.
diff --git a/AMDiS/src/ProblemStatBase.h b/AMDiS/src/ProblemStatBase.h
index 264de49f..12d9633c 100644
--- a/AMDiS/src/ProblemStatBase.h
+++ b/AMDiS/src/ProblemStatBase.h
@@ -102,8 +102,11 @@ namespace AMDiS {
/** \brief
* Assembling of system matrices and vectors after coarsening.
+ * By the last two parameters, assembling can be restricted to either
+ * matrices or vectors only.
*/
- virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag) = 0;
+ virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
+ bool assembleMatrix, bool assembleVector) = 0;
/** \brief
* Refinement of the mesh.
@@ -118,7 +121,7 @@ namespace AMDiS {
/** \brief
* Solves the assembled system. The result is an approximative solution.
*/
- virtual void solve(AdaptInfo *adaptInfo) = 0;
+ virtual void solve(AdaptInfo *adaptInfo, bool fixedMatrix) = 0;
/** \brief
* A posteriori error estimation of the calculated solution. Should store
diff --git a/AMDiS/src/ProblemVec.cc b/AMDiS/src/ProblemVec.cc
index ea6c7e70..89f7b4da 100644
--- a/AMDiS/src/ProblemVec.cc
+++ b/AMDiS/src/ProblemVec.cc
@@ -506,7 +506,7 @@ namespace AMDiS {
{
}
- void ProblemVec::solve(AdaptInfo *adaptInfo)
+ void ProblemVec::solve(AdaptInfo *adaptInfo, bool fixedMatrix)
{
FUNCNAME("Problem::solve()");
@@ -520,7 +520,8 @@ namespace AMDiS {
#endif
clock_t first = clock();
- int iter = solver_->solve(matVec_, solution_, rhs_, leftPrecon_, rightPrecon_);
+ int iter = solver_->solve(matVec_, solution_, rhs_,
+ leftPrecon_, rightPrecon_, fixedMatrix);
#ifdef _OPENMP
INFO(info_, 8)("solution of discrete system needed %.5f seconds system time / %.5f seconds wallclock time\n",
@@ -648,7 +649,8 @@ namespace AMDiS {
return StandardProblemIteration::oneIteration(adaptInfo, toDo);
}
- void ProblemVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag)
+ void ProblemVec::buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
+ bool asmMatrix, bool asmVector)
{
FUNCNAME("ProblemVec::buildAfterCoarsen()");
@@ -687,7 +689,7 @@ namespace AMDiS {
if ((*systemMatrix_)[i][j]) {
// The matrix should not be deleted, if it was assembled before
// and it is marked to be assembled only once.
- if (!(assembleMatrixOnlyOnce_[i][j] && assembledMatrix_[i][j])) {
+ if (!(assembleMatrixOnlyOnce_[i][j] && assembledMatrix_[i][j]) && asmMatrix) {
(*systemMatrix_)[i][j]->clear();
}
}
@@ -712,6 +714,10 @@ namespace AMDiS {
assembleMatrix = false;
}
+ if (!asmMatrix) {
+ assembleMatrix = false;
+ }
+
// If the matrix should not be assembled, the rhs vector has to be considered.
// This will be only done, if i == j. So, if both is not true, we can jump
// to the next matrix.
@@ -727,7 +733,7 @@ namespace AMDiS {
assembleOnOneMesh(componentSpaces[i],
assembleFlag,
assembleMatrix ? matrix : NULL,
- (i == j) ? rhs_->getDOFVector(i) : NULL);
+ ((i == j) && asmVector) ? rhs_->getDOFVector(i) : NULL);
// std::cout << "--- Finished ---\n";
// std::cout << "--- Dim: " << matrix->getUsedSize() << "x" << matrix->getNumCols() << "---\n";
} else {
@@ -735,7 +741,7 @@ namespace AMDiS {
assembleOnDifMeshes(componentSpaces[i], componentSpaces[j],
assembleFlag,
assembleMatrix ? matrix : NULL,
- (i == j) ? rhs_->getDOFVector(i) : NULL);
+ ((i == j) && asmVector) ? rhs_->getDOFVector(i) : NULL);
// std::cout << "--- Finished ---\n";
// std::cout << "--- Dim: " << matrix->getUsedSize() << "x" << matrix->getNumCols() << "---\n";
}
diff --git a/AMDiS/src/ProblemVec.h b/AMDiS/src/ProblemVec.h
index e3c6b650..e840b33d 100644
--- a/AMDiS/src/ProblemVec.h
+++ b/AMDiS/src/ProblemVec.h
@@ -151,7 +151,7 @@ namespace AMDiS {
* Implementation of ProblemStatBase::solve(). Deligates the solving
* of problems system to \ref solver.
*/
- virtual void solve(AdaptInfo *adaptInfo);
+ virtual void solve(AdaptInfo *adaptInfo, bool fixedMatrix = false);
/** \brief
* Implementation of ProblemStatBase::estimate(). Deligates the estimation
@@ -191,9 +191,12 @@ namespace AMDiS {
/** \brief
* Implementation of ProblemStatBase::buildAfterCoarsen().
- * Assembles \ref A and \ref rhs.
+ * Assembles \ref A and \ref rhs. With the last two parameters, assembling
+ * can be restricted to matrices or vectors only.
*/
- virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag);
+ virtual void buildAfterCoarsen(AdaptInfo *adaptInfo, Flag flag,
+ bool assembleMatrix = true,
+ bool assembleVector = true);
/** \brief
* Determines the execution order of the single adaption steps. If adapt is
diff --git a/AMDiS/src/SolutionDataStorage.h b/AMDiS/src/SolutionDataStorage.h
index 15319ac1..4cd0b0eb 100644
--- a/AMDiS/src/SolutionDataStorage.h
+++ b/AMDiS/src/SolutionDataStorage.h
@@ -82,8 +82,9 @@ namespace AMDiS {
}
void deleteFeSpace(std::vector<FiniteElemSpace*> feSpaces) {
- for (int i = 0; i < static_cast<int>(feSpaces.size()); i++)
+ for (int i = 0; i < static_cast<int>(feSpaces.size()); i++) {
DELETE feSpaces[i];
+ }
feSpaces.empty();
}
diff --git a/AMDiS/src/SolutionDataStorage.hh b/AMDiS/src/SolutionDataStorage.hh
index 62e810f1..1bf73e5c 100644
--- a/AMDiS/src/SolutionDataStorage.hh
+++ b/AMDiS/src/SolutionDataStorage.hh
@@ -61,9 +61,13 @@ namespace AMDiS {
double timestamp,
typename SolutionHelper<T>::type feSpace)
{
+ FUNCNAME("SolutionDataStorage<T>::push()");
+
push(solution, timestamp);
// Store fe space only, if we do not have a fixed fe space.
+ TEST_EXIT(!fixedFESpace)("push wit fe space not possible!\n");
+
if (!fixedFESpace) {
feSpaces.push_back(feSpace);
}
@@ -111,8 +115,10 @@ namespace AMDiS {
DELETE solutions[i];
}
- for (int i = 0; i < static_cast<int>(feSpaces.size()); i++) {
- deleteFeSpace(feSpaces[i]);
+ if (!fixedFESpace) {
+ for (int i = 0; i < static_cast<int>(feSpaces.size()); i++) {
+ deleteFeSpace(feSpaces[i]);
+ }
}
solutions.empty();
diff --git a/AMDiS/src/StandardProblemIteration.cc b/AMDiS/src/StandardProblemIteration.cc
index 706b4bea..ac9742c9 100644
--- a/AMDiS/src/StandardProblemIteration.cc
+++ b/AMDiS/src/StandardProblemIteration.cc
@@ -24,7 +24,10 @@ namespace AMDiS {
Flag flag = buildAndAdapt(adaptInfo, toDo);
if (toDo.isSet(SOLVE))
- problem_->solve(adaptInfo);
+ problem_->solve(adaptInfo, false);
+
+ if (toDo.isSet(SOLVE_RHS))
+ problem_->solve(adaptInfo, true);
if (toDo.isSet(ESTIMATE))
problem_->estimate(adaptInfo);
@@ -72,7 +75,11 @@ namespace AMDiS {
}
if (toDo.isSet(BUILD)) {
- problem_->buildAfterCoarsen(adaptInfo, markFlag);
+ problem_->buildAfterCoarsen(adaptInfo, markFlag, true, true);
+ }
+
+ if (toDo.isSet(BUILD_RHS)) {
+ problem_->buildAfterCoarsen(adaptInfo, markFlag, false, true);
}
return flag;
diff --git a/AMDiS/src/StandardProblemIteration.h b/AMDiS/src/StandardProblemIteration.h
index bb9e8d58..b1b6faf9 100644
--- a/AMDiS/src/StandardProblemIteration.h
+++ b/AMDiS/src/StandardProblemIteration.h
@@ -71,14 +71,14 @@ namespace AMDiS {
*/
int getNumProblems() {
return 1;
- };
+ }
/** \brief
* implementation of \ref ProblemIterationInterface::getProblem(int)
*/
ProblemStatBase *getProblem(int number = 0) {
return problem_;
- };
+ }
/** \brief
* Returns the name of the problem.
diff --git a/AMDiS/src/TFQMR.h b/AMDiS/src/TFQMR.h
index 67087e1d..e127e267 100644
--- a/AMDiS/src/TFQMR.h
+++ b/AMDiS/src/TFQMR.h
@@ -72,7 +72,8 @@ namespace AMDiS {
/** \brief
* realisation of OEMSolver::solveSystem
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv,
+ VectorType *x, VectorType *b, bool reuseMatrix);
/** \brief
* realisation of OEMSolver::init
diff --git a/AMDiS/src/TFQMR.hh b/AMDiS/src/TFQMR.hh
index a3773e75..d2ae5e92 100644
--- a/AMDiS/src/TFQMR.hh
+++ b/AMDiS/src/TFQMR.hh
@@ -43,8 +43,9 @@ namespace AMDiS {
template<typename VectorType>
int TFQMR<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x,
- VectorType *b)
+ VectorType *x,
+ VectorType *b,
+ bool reuseMatrix)
{
DOFMatrix *m = (dynamic_cast<StandardMatVec<DOFMatrix, VectorType> *>(matVec))->getMatrix();
diff --git a/AMDiS/src/UmfPackSolver.h b/AMDiS/src/UmfPackSolver.h
index fcbc07c2..6bc57c89 100644
--- a/AMDiS/src/UmfPackSolver.h
+++ b/AMDiS/src/UmfPackSolver.h
@@ -82,7 +82,7 @@ namespace AMDiS {
void init() {
p = this->vectorCreator->create();
r = this->vectorCreator->create();
- };
+ }
/** \brief
* Implements OEMSolver<VectorType>::exit().
@@ -90,12 +90,12 @@ namespace AMDiS {
void exit() {
this->vectorCreator->free(p);
this->vectorCreator->free(r);
- };
+ }
/** \brief
* Implements OEMSolver<VectorType>::solve().
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b, bool reuseMatrix);
private:
/** \brief
@@ -109,6 +109,11 @@ namespace AMDiS {
*/
void *symbolic_;
+ /** \brief
+ * Stores the result of umfpack_di_numeric, if multiple right hand sides will occure.
+ */
+ void *numeric;
+
/** \brief
* If the symbolic analysis should be done only once (for example, if the matrices
* to solve have all the same pattern because of no adaptivity), this variable
@@ -116,6 +121,12 @@ namespace AMDiS {
* of solveSystem().
*/
int storeSymbolic_;
+
+ /** \brief
+ * If not zero, Umfpack is prepared to solve multiple following systems with equal
+ * system matrix but different right hand side vectors.
+ */
+ int multipleRhs;
};
}
diff --git a/AMDiS/src/UmfPackSolver.hh b/AMDiS/src/UmfPackSolver.hh
index 42e47c8b..2d7655db 100644
--- a/AMDiS/src/UmfPackSolver.hh
+++ b/AMDiS/src/UmfPackSolver.hh
@@ -13,9 +13,11 @@ namespace AMDiS {
UmfPackSolver<VectorType>::UmfPackSolver(std::string name)
: OEMSolver<VectorType>(name),
symbolic_(NULL),
- storeSymbolic_(0)
+ storeSymbolic_(0),
+ multipleRhs(0)
{
GET_PARAMETER(0, name + "->store symbolic", "%d", &storeSymbolic_);
+ GET_PARAMETER(0, name + "->multiple rhs", "%d", &multipleRhs);
}
template<typename VectorType>
@@ -28,7 +30,8 @@ namespace AMDiS {
template<typename VectorType>
int UmfPackSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x, VectorType *b)
+ VectorType *x, VectorType *b,
+ bool reuseMatrix)
{
FUNCNAME("UmfPackSolver::solveSystem()");
@@ -114,7 +117,6 @@ namespace AMDiS {
}
- void *numeric;
double Control[UMFPACK_CONTROL];
double Info[UMFPACK_INFO];
int status = 0;
@@ -128,15 +130,24 @@ namespace AMDiS {
if (!(storeSymbolic_) ||
(storeSymbolic_ && (symbolic_ == NULL))) {
- status = umfpack_di_symbolic(newMatrixSize, newMatrixSize, Ap, Ai, Ax, &symbolic_, Control, Info);
+ status = umfpack_di_symbolic(newMatrixSize, newMatrixSize, Ap, Ai, Ax,
+ &symbolic_, Control, Info);
if (!status == UMFPACK_OK) {
ERROR_EXIT("UMFPACK Error in function umfpack_di_symbolic");
}
}
-
- status = umfpack_di_numeric(Ap, Ai, Ax, symbolic_, &numeric, Control, Info);
- if (!status == UMFPACK_OK) {
- ERROR_EXIT("UMFPACK Error in function umfpack_di_numeric");
+
+ if (!((multipleRhs != 0) && (reuseMatrix))) {
+ // With multiple rhs are allowed, but not this time, free the matrix
+ // factorization of last solution step,
+ if ((multipleRhs != 0) && !reuseMatrix) {
+ umfpack_di_free_numeric(&numeric);
+ }
+
+ status = umfpack_di_numeric(Ap, Ai, Ax, symbolic_, &numeric, Control, Info);
+ if (!status == UMFPACK_OK) {
+ ERROR_EXIT("UMFPACK Error in function umfpack_di_numeric");
+ }
}
status = umfpack_di_solve(UMFPACK_A, Ap, Ai, Ax, xvec, bvec, numeric, Control, Info);
@@ -147,7 +158,10 @@ namespace AMDiS {
if (!storeSymbolic_) {
umfpack_di_free_symbolic(&symbolic_);
}
- umfpack_di_free_numeric(&numeric);
+ if (multipleRhs == 0) {
+ std::cout << "FREE THE STUFF" << std::endl;
+ umfpack_di_free_numeric(&numeric);
+ }
// Copy and resort solution.
diff --git a/AMDiS/src/VecSymSolver.h b/AMDiS/src/VecSymSolver.h
index 632719ef..130d50b3 100644
--- a/AMDiS/src/VecSymSolver.h
+++ b/AMDiS/src/VecSymSolver.h
@@ -89,7 +89,8 @@ namespace AMDiS {
/** \brief
* Implements OEMSolver<VectorType>::solve().
*/
- int solveSystem(MatVecMultiplier<VectorType> *mv, VectorType *x, VectorType *b);
+ int solveSystem(MatVecMultiplier<VectorType> *mv,
+ VectorType *x, VectorType *b, bool reuseMatrix);
private:
VectorType *r, *p;
diff --git a/AMDiS/src/VecSymSolver.hh b/AMDiS/src/VecSymSolver.hh
index 3e5f522d..76a1cd6f 100644
--- a/AMDiS/src/VecSymSolver.hh
+++ b/AMDiS/src/VecSymSolver.hh
@@ -227,7 +227,7 @@ namespace AMDiS {
template<typename VectorType>
int VecSymSolver<VectorType>::solveSystem(MatVecMultiplier<VectorType> *matVec,
- VectorType *x, VectorType *b)
+ VectorType *x, VectorType *b, bool reuseMatrix)
{
FUNCNAME("VecSymSolver::solveSystem()");
--
GitLab