Commit f0f8d33c authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* Several changes to make Rosenbrock method possible to solve systems with...

* Several changes to make Rosenbrock method possible to solve systems with equal matrix but different right hand side vectors

parent ecc316ac
......@@ -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");
......
......@@ -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()");
......
......@@ -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
......
......@@ -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
......
......@@ -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;
......
......@@ -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()");
......
......@@ -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.
......
......@@ -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
......
......@@ -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";
}
......
......@@ -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
......
......@@ -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();
}
......
......@@ -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();
......
......@@ -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;
......
......@@ -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.
......
......@@ -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
......
......@@ -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();
......
......@@ -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;
};
}
......
......@@ -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.
......
......@@ -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;
......
......@@ -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()");
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment