Commit 5e529972 authored by Thomas Witkowski's avatar Thomas Witkowski

Fixed UMFPACK problem for solving with multiple RHS vectors.

parent 169b9d9a
......@@ -51,12 +51,14 @@ namespace AMDiS {
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)
VectorialMapper& m,
bool createMatrixData,
bool storeMatrixData)
{
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, createMatrixData, storeMatrixData);
}
......@@ -96,9 +98,11 @@ namespace AMDiS {
int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
SystemVector& b,
VectorialMapper& m)
VectorialMapper& m,
bool createMatrixData,
bool storeMatrixData)
{
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, createMatrixData, storeMatrixData);
}
~ITL_OEMSolver_para() {}
......
......@@ -42,7 +42,7 @@ namespace AMDiS {
template< typename Mapper, bool IsDist>
void init(Mapper& mapper, boost::mpl::bool_<IsDist>)
{
matrix.change_dim(mapper.getNumRows(), mapper.getNumCols());
matrix.change_dim(mapper.getNumRows(), mapper.getNumCols());
}
#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(HAVE_PARALLEL_MTL4)
......@@ -66,16 +66,18 @@ namespace AMDiS {
}
template< typename Matrix, typename Vector, typename Mapper >
int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper)
int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper,
bool createMatrixData,
bool storeMatrixData)
{
FUNCNAME("MTL4Solver::solve()");
Timer t;
if (num_rows(matrix) == 0 || !getMultipleRhs()) {
if (createMatrixData) {
init(mapper, mtl::traits::is_distributed<MTLMatrix>());
set_to_zero(matrix);
MatMap< const Matrix, Mapper > matMap(A,mapper);
MatMap< const Matrix, Mapper > matMap(A, mapper);
matrix << matMap;
worker.init(matrix);
}
......@@ -120,9 +122,11 @@ namespace AMDiS {
virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
SystemVector& b,
VectorialMapper& m)
VectorialMapper& m,
bool createMatrixData,
bool storeMatrixData)
{
return solve(A,x,b,m);
return solve(A, x, b, m, createMatrixData, storeMatrixData);
}
......
......@@ -71,8 +71,7 @@ namespace AMDiS {
residual(0),
print_cycle(100),
iterations(-1),
error(-1),
multipleRhs(false)
error(-1)
{}
///
......@@ -95,7 +94,9 @@ namespace AMDiS {
virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
SystemVector& b,
VectorialMapper&)
VectorialMapper&,
bool createMatrixData,
bool storeMatrixData)
{
FUNCNAME("OEMSolver::solveSystem()");
TEST_EXIT(false)
......@@ -105,10 +106,12 @@ namespace AMDiS {
inline int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
SystemVector& b)
SystemVector& b,
bool createMatrixData,
bool storeMatrixData)
{
VectorialMapper mapper(A);
return solveSystem(A, x, b, mapper);
return solveSystem(A, x, b, mapper, createMatrixData, storeMatrixData);
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
......@@ -182,12 +185,6 @@ namespace AMDiS {
return relative;
}
/// Returns \ref multipleRhs
inline bool getMultipleRhs() const
{
return multipleRhs;
}
/** \} */
......@@ -235,11 +232,6 @@ namespace AMDiS {
info = i;
}
inline void setMultipleRhs(bool b)
{
multipleRhs = b;
}
/** \} */
protected:
......@@ -269,11 +261,6 @@ namespace AMDiS {
/// Error code in last solver (not set by UmfPack)
int error;
/// 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;
};
/**
......
......@@ -201,9 +201,9 @@ namespace AMDiS {
if (solver) {
WARNING("solver already created\n");
} else {
if (initFlag.isSet(INIT_SOLVER)) {
if (initFlag.isSet(INIT_SOLVER))
createSolver();
}
if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
TEST_EXIT(!solver)("solver already created\n");
solver = adoptProblem->getSolver();
......@@ -576,11 +576,12 @@ namespace AMDiS {
void ProblemStatSeq::doOtherStuff()
{
}
{}
void ProblemStatSeq::solve(AdaptInfo *adaptInfo, bool, bool)
void ProblemStatSeq::solve(AdaptInfo *adaptInfo,
bool createMatrixData,
bool storeMatrixData)
{
FUNCNAME("Problem::solve()");
......@@ -590,8 +591,9 @@ namespace AMDiS {
}
clock_t first = clock();
solver->solveSystem(solverMatrix, *solution, *rhs);
solver->solveSystem(solverMatrix, *solution, *rhs,
createMatrixData, storeMatrixData);
INFO(info, 8)("solution of discrete system needed %.5f seconds\n",
TIME_USED(first, clock()));
......@@ -709,9 +711,6 @@ namespace AMDiS {
{
FUNCNAME("ProblemStat::buildAfterCoarsen()");
// buildAfterCoarsen_sebastianMode(adaptInfo, flag);
// return;
if (dualMeshTraverseRequired()) {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
ERROR_EXIT("Dual mesh assemble does not work in parallel code!\n");
......@@ -858,13 +857,10 @@ namespace AMDiS {
}
if (asmMatrix) {
solver->setMultipleRhs(false);
solverMatrix.setMatrix(*systemMatrix);
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
} else {
solver->setMultipleRhs(true);
}
}
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
......@@ -878,187 +874,6 @@ namespace AMDiS {
}
void ProblemStatSeq::buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag)
{
FUNCNAME("ProblemStat::buildAfterCoarsen()");
clock_t first = clock();
for (int i = 0; i < static_cast<int>(meshes.size()); i++)
meshes[i]->dofCompress();
Flag assembleFlag =
flag |
(*systemMatrix)[0][0]->getAssembleFlag() |
rhs->getDOFVector(0)->getAssembleFlag() |
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS |
Mesh::FILL_DET |
Mesh::FILL_GRD_LAMBDA |
Mesh::FILL_NEIGH;
if (useGetBound)
assembleFlag |= Mesh::FILL_BOUND;
traverseInfo.updateStatus();
// Used to calculate the overall number of non zero entries.
int nnz = 0;
/// === INITIALIZE ===
for (int i = 0; i < nComponents; i++) {
MSG("%d DOFs for %s\n",
componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[i]->getName().c_str());
rhs->getDOFVector(i)->set(0.0);
for (int j = 0; j < nComponents; j++) {
// Only if this variable is true, the current matrix will be assembled.
bool assembleMatrix = true;
// The DOFMatrix which should be assembled (or not, if assembleMatrix
// will be set to false).
DOFMatrix *matrix = (*systemMatrix)[i][j];
if (matrix)
matrix->calculateNnz();
// If the matrix was assembled before and it is marked to be assembled
// only once, it will not be assembled.
if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j]) {
assembleMatrix = false;
} else if (matrix) {
matrix->getBaseMatrix().
change_dim(componentSpaces[i]->getAdmin()->getUsedSize(),
componentSpaces[j]->getAdmin()->getUsedSize());
set_to_zero(matrix->getBaseMatrix());
}
// If there is no DOFMatrix, e.g., if it is completly 0, do not assemble.
if (!matrix || !assembleMatrix)
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.
if (!assembleMatrix && i != j) {
if (matrix)
nnz += matrix->getBaseMatrix().nnz();
continue;
}
if (assembleMatrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->initMatrix(matrix);
if (matrix && assembleMatrix)
matrix->startInsertion(matrix->getNnz());
}
}
// === TRAVERSE ===
Mesh *mesh = componentMeshes[0];
const FiniteElemSpace *feSpace = componentSpaces[0];
const BasisFunction *basisFcts = feSpace->getBasisFcts();
ElementMatrix elMat(basisFcts->getNumber(), basisFcts->getNumber());
ElementMatrix tmpElMat(elMat);
ElementVector elVec(basisFcts->getNumber());
ElementVector tmpElVec(elVec);
TraverseStack stack;
BoundaryType *bound =
useGetBound ? new BoundaryType[basisFcts->getNumber()] : NULL;
ElInfo *elInfo = stack.traverseFirst(mesh, -1, assembleFlag);
while (elInfo) {
if (useGetBound)
basisFcts->getBound(elInfo, bound);
for (map<Operator*, vector<OperatorPos> >::iterator opIt = operators.begin();
opIt != operators.end(); ++opIt) {
if (opIt->first->getNeedDualTraverse() == true)
continue;
if (opFlags[opIt->first].isSet(Operator::MATRIX_OPERATOR)) {
set_to_zero(elMat);
opIt->first->getElementMatrix(elInfo, elMat, 1.0);
}
if (opFlags[opIt->first].isSet(Operator::VECTOR_OPERATOR)) {
set_to_zero(elVec);
opIt->first->getElementVector(elInfo, elVec, 1.0);
}
for (vector<OperatorPos>::iterator posIt = opIt->second.begin();
posIt != opIt->second.end(); ++posIt) {
if (posIt->operatorType.isSet(Operator::MATRIX_OPERATOR)) {
if (!posIt->factor || *(posIt->factor) == 1.0) {
(*systemMatrix)[posIt->row][posIt->col]->addElementMatrix(elMat, bound, elInfo, NULL);
} else {
tmpElMat = *(posIt->factor) * elMat;
(*systemMatrix)[posIt->row][posIt->col]->addElementMatrix(tmpElMat, bound, elInfo, NULL);
}
}
if (posIt->operatorType.isSet(Operator::VECTOR_OPERATOR)) {
if (!posIt->factor || *(posIt->factor) == 1.0) {
rhs->getDOFVector(posIt->row)->addElementVector(1.0, elVec, bound, elInfo);
} else {
tmpElVec = *(posIt->factor) * elVec;
rhs->getDOFVector(posIt->row)->addElementVector(1.0, tmpElVec, bound, elInfo);
}
}
}
}
elInfo = stack.traverseNext(elInfo);
}
if (useGetBound)
delete [] bound;
// === FINELIZE ===
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
bool assembleMatrix = true;
DOFMatrix *matrix = (*systemMatrix)[i][j];
if (assembleMatrixOnlyOnce[i][j] && assembledMatrix[i][j])
assembleMatrix = false;
if (!matrix || !assembleMatrix)
assembleMatrix = false;
if (!assembleMatrix && i != j)
continue;
assembledMatrix[i][j] = true;
if (assembleMatrix) {
matrix->clearDirichletRows();
matrix->finishInsertion();
}
if (assembleMatrix && matrix->getBoundaryManager())
matrix->getBoundaryManager()->exitMatrix(matrix);
if (matrix)
nnz += matrix->getBaseMatrix().nnz();
}
assembleBoundaryConditions(rhs->getDOFVector(i),
solution->getDOFVector(i),
componentMeshes[i],
assembleFlag);
}
solverMatrix.setMatrix(*systemMatrix);
INFO(info, 8)("fillin of assembled matrix: %d\n", nnz);
INFO(info, 8)("buildAfterCoarsen needed %.5f seconds\n",
TIME_USED(first, clock()));
}
bool ProblemStatSeq::dualMeshTraverseRequired()
{
FUNCNAME("ProblemStat::dualMeshTraverseRequired()");
......
......@@ -140,8 +140,6 @@ namespace AMDiS {
bool assembleMatrix = true,
bool assembleVector = true);
void buildAfterCoarsen_sebastianMode(AdaptInfo *adaptInfo, Flag flag);
bool dualMeshTraverseRequired();
void dualAssemble(AdaptInfo *adaptInfo, Flag flag,
......
......@@ -64,8 +64,12 @@ namespace AMDiS {
void init(const Matrix& A)
{
FUNCNAME("Umfpack_runner::init()")
if(solver != NULL)
delete solver;
if (solver != NULL) {
delete solver;
solver = NULL;
}
solver = new mtl::matrix::umfpack::solver<matrix_type>(A, symmetric_strategy, alloc_init);
}
......@@ -75,12 +79,10 @@ namespace AMDiS {
FUNCNAME("Umfpack_runner::solve()");
TEST_EXIT(solver != NULL)("The umfpack solver was not initialized\n");
#if 0
if (!oem.getMultipleRhs()) {
if (store_symbolic)
solver->update_numeric();
else
solver->update();
}
if (store_symbolic)
solver->update_numeric();
else
solver->update();
#endif
int code = (*solver)(x, b);
......@@ -101,8 +103,10 @@ namespace AMDiS {
~Umfpack_runner()
{
if (solver != NULL)
if (solver != NULL) {
delete solver;
solver == NULL;
}
}
private:
......@@ -142,12 +146,12 @@ namespace AMDiS {
int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
SystemVector& x,
SystemVector& b,
VectorialMapper& m)
VectorialMapper& m,
bool createMatrixData,
bool storeMatrixData)
{
return solve(A,x,b,m);
}
private:
return solve(A, x, b, m, createMatrixData, storeMatrixData);
}
};
}
......
......@@ -84,14 +84,11 @@ namespace AMDiS {
MSG("iter. | this->residual | red. | n |\n");
for (iter = 1; iter <= this->maxIter; iter++) {
if (iter == 1 || (buildCycle > 0 && (iter-1) % buildCycle == 0)) {
this->getLinearSolver()->setMultipleRhs(false);
// Assemble DF(x) and F(x)
// Assemble DF(x) and F(x)
if (iter == 1 || (buildCycle > 0 && (iter-1) % buildCycle == 0))
prob->buildAfterCoarsen(adaptInfo, 0, true, true);
} else {
this->getLinearSolver()->setMultipleRhs(true);
else
prob->buildAfterCoarsen(adaptInfo, 0, false, true);
}
// Initial guess is zero
b->set(0.0);
......
......@@ -91,14 +91,11 @@ namespace AMDiS {
MSG("iter. | this->residual | red. | n | lambda |\n");
for (iter = 1; iter <= this->maxIter; iter++) {
if (iter == 1 || (buildCycle > 0 && (iter-1) % buildCycle == 0)) {
this->getLinearSolver()->setMultipleRhs(false);
// Assemble DF(x) and F(x)
// Assemble DF(x) and F(x)
if (iter == 1 || (buildCycle > 0 && (iter-1) % buildCycle == 0))
prob->buildAfterCoarsen(adaptInfo, 0, true, true);
} else {
this->getLinearSolver()->setMultipleRhs(true);
prob->buildAfterCoarsen(adaptInfo, 0, false, true);
}
else
prob->buildAfterCoarsen(adaptInfo, 0, false, true);
// Initial guess is zero
b->set(0.0);
......
......@@ -111,7 +111,8 @@ namespace AMDiS {
/// Solves the non linear system. Must be overriden in sub classes.
virtual int nlsolve(SolverMatrix<Matrix<DOFMatrix*> >& matVec,
SystemVector& x, SystemVector& rhs,
SystemVector& x,
SystemVector& rhs,
AdaptInfo *adaptInfo,
ProblemStat *prob) = 0;
......@@ -119,11 +120,12 @@ namespace AMDiS {
virtual void exit() = 0;
virtual int solveLinearSystem(SolverMatrix<Matrix<DOFMatrix*> >& mat,
SystemVector &x, SystemVector &b)
SystemVector &x,
SystemVector &b)
{
FUNCNAME("NonLinSolver::solveLinearSystem()");
TEST_EXIT(linSolver)("no solver\n");
return linSolver->solveSystem(mat, x, b);
return linSolver->solveSystem(mat, x, b, true, false);
}
protected:
......
......@@ -762,24 +762,50 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::createFetiKsp()");
MSG("START A\n");
// === Create FETI-DP solver object. ===
fetiData.mat_lagrange = &mat_lagrange;
fetiData.subSolver = subdomain;
fetiData.ksp_schur_primal = &ksp_schur_primal;
FetiData &data = (!enableStokesMode ? fetiData : fetiDataInterface);
data.mat_lagrange = &mat_lagrange;
data.subSolver = subdomain;
data.ksp_schur_primal = &ksp_schur_primal;
localDofMap.createVec(data.tmp_vec_b, nGlobalOverallInterior);
lagrangeMap.createVec(data.tmp_vec_lagrange);
primalDofMap.createVec(data.tmp_vec_primal);
localDofMap.createVec(fetiData.tmp_vec_b, nGlobalOverallInterior);
lagrangeMap.createVec(fetiData.tmp_vec_lagrange);
primalDofMap.createVec(fetiData.tmp_vec_primal);
if (enableStokesMode == false) {
MatCreateShell(mpiCommGlobal,
lagrangeMap.getRankDofs(),
lagrangeMap.getRankDofs(),
lagrangeMap.getOverallDofs(),
lagrangeMap.getOverallDofs(),
&fetiData, &mat_feti);
MatShellSetOperation(mat_feti, MATOP_MULT,
(void(*)(void))petscMultMatFeti);
} else {
MSG("TEST A0\n");
MatCreateShell(mpiCommGlobal,
lagrangeMap.getRankDofs(),
lagrangeMap.getRankDofs(),
lagrangeMap.getOverallDofs(),
lagrangeMap.getOverallDofs(),
&fetiData, &mat_feti);
MatShellSetOperation(mat_feti, MATOP_MULT, (void(*)(void))petscMultMatFeti);
fetiDataInterface.mat_interior_interface =
getMatInteriorCoarseByComponent(2);
fetiDataInterface.mat_interface_interior =
getMatCoarseInteriorByComponent(2);
fetiDataInterface.mat_primal_interface = getMatCoarse(0, 1);
fetiDataInterface.mat_interface_primal = getMatCoarse(1, 0);
MSG("TEST A1\n");
MatCreateShell(mpiCommGlobal,
interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
interfaceDofMap.getRankDofs() + lagrangeMap.getRankDofs(),
interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
interfaceDofMap.getOverallDofs() + lagrangeMap.getOverallDofs(),
&fetiDataInterface, &mat_feti);
MatShellSetOperation(mat_feti, MATOP_MULT,
(void(*)(void))petscMultMatFetiInterface);
MSG("TEST A2\n");
}
KSPCreate(mpiCommGlobal, &ksp_feti);
KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
......@@ -889,6 +915,8 @@ namespace AMDiS {
default:
break;
}
MSG("END A\n");
}
......@@ -1414,43 +1442,53 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::solveReducedFetiMatrix()");
// Some temporary vectors.
Vec tmp_b0, tmp_b1, tmp_lagrange0, tmp_primal0, tmp_primal1;
// === Some temporary vectors. ===
Vec tmp_b0, tmp_b1;
localDofMap.createVec(tmp_b0, nGlobalOverallInterior);
localDofMap.createVec(tmp_b1, nGlobalOverallInterior);
Vec tmp_primal0, tmp_primal1;
primalDofMap.createVec(tmp_primal0);
primalDofMap.createVec(tmp_primal1);
// RHS vector.
Vec vecRhs, vecSol;
MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange0);
MatGetVecs(mat_lagrange, PETSC_NULL, &vecRhs);
MatGetVecs(mat_lagrange, PETSC_NULL, &vecSol);
Vec tmp_lagrange;
MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange);
// === Create RHS and solution vectors. ===
Vec vecRhs, vecSol;
Vec vecRhsLagrange, vecRhsInterface, vecSolLagrange, vecSolInterface;