Commit bb1e643b authored by Thomas Witkowski's avatar Thomas Witkowski

SubDomainSolver interface for sequentiel sub solvers in FETI-DP method is...

SubDomainSolver interface for sequentiel sub solvers in FETI-DP method is fixed and should work now.
parent bc02aae5
......@@ -95,6 +95,8 @@ namespace AMDiS {
petscSolver->solvePetscMatrix(*solution, adaptInfo);
petscSolver->destroyVectorData();
if (!storeMatrixData)
petscSolver->destroyMatrixData();
......
......@@ -79,6 +79,9 @@ namespace AMDiS {
/// Destroys all matrix data structures.
virtual void destroyMatrixData() = 0;
/// Detroys all vector data structures.
virtual void destroyVectorData() = 0;
virtual Flag getBoundaryDofRequirement()
{
return 0;
......
......@@ -30,10 +30,10 @@ namespace AMDiS {
MatShellGetContext(mat, &ctx);
SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
MatMult(*(data->mat_primal_primal), x, y);
MatMult(data->subSolver->getMatIntCoarse(), x, data->tmp_vec_b);
data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
MatMult(data->subSolver->getMatCoarseCoarse(), x, y);
VecAXPBY(y, -1.0, 1.0, data->tmp_vec_primal);
return 0;
......@@ -51,13 +51,13 @@ namespace AMDiS {
FetiData* data = static_cast<FetiData*>(ctx);
MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b);
data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal);
MatMult(data->subSolver->getMatCoarseInt(), data->tmp_vec_b, data->tmp_vec_primal);
KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, data->tmp_vec_primal);
MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b);
data->fetiSolver->solveLocalProblem(data->tmp_vec_b, data->tmp_vec_b);
MatMult(data->subSolver->getMatIntCoarse(), data->tmp_vec_primal, data->tmp_vec_b);
data->subSolver->solveGlobal(data->tmp_vec_b, data->tmp_vec_b);
MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
......@@ -487,10 +487,7 @@ namespace AMDiS {
if (schurPrimalSolver == 0) {
MSG("Create iterative schur primal solver!\n");
schurPrimalData.mat_primal_primal = &(subDomainSolver->getMatCoarseCoarse());
schurPrimalData.mat_primal_b = &(subDomainSolver->getMatCoarseInt());
schurPrimalData.mat_b_primal = &(subDomainSolver->getMatIntCoarse());
schurPrimalData.fetiSolver = this;
schurPrimalData.subSolver = subDomainSolver;
VecCreateMPI(mpiComm,
localDofMap.getRankDofs(), localDofMap.getOverallDofs(),
......@@ -579,20 +576,20 @@ namespace AMDiS {
MatAssemblyBegin(matBPi, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(matBPi, MAT_FINAL_ASSEMBLY);
MatDuplicate(subDomainSolver->getMatCoarseCoarse(), MAT_COPY_VALUES, &mat_schur_primal);
MatMatMult(subDomainSolver->getMatCoarseInt(), matBPi, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matPrimal);
MatAXPY(subDomainSolver->getMatCoarseCoarse(), -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
MatAXPY(mat_schur_primal, -1.0, matPrimal, DIFFERENT_NONZERO_PATTERN);
MatDestroy(&matPrimal);
MatDestroy(&matBPi);
MatInfo minfo;
MatGetInfo(subDomainSolver->getMatCoarseCoarse(), MAT_GLOBAL_SUM, &minfo);
MatGetInfo(mat_schur_primal, MAT_GLOBAL_SUM, &minfo);
MSG("Schur primal matrix nnz = %f\n", minfo.nz_used);
KSPCreate(mpiComm, &ksp_schur_primal);
KSPSetOperators(ksp_schur_primal,
subDomainSolver->getMatCoarseCoarse(),
subDomainSolver->getMatCoarseCoarse(),
KSPSetOperators(ksp_schur_primal, mat_schur_primal, mat_schur_primal,
SAME_NONZERO_PATTERN);
KSPSetOptionsPrefix(ksp_schur_primal, "schur_primal_");
KSPSetType(ksp_schur_primal, KSPPREONLY);
......@@ -613,19 +610,14 @@ namespace AMDiS {
FUNCNAME("PetscSolverFeti::destroySchurPrimal()");
if (schurPrimalSolver == 0) {
schurPrimalData.mat_primal_primal = PETSC_NULL;
schurPrimalData.mat_primal_b = PETSC_NULL;
schurPrimalData.mat_b_primal = PETSC_NULL;
schurPrimalData.fetiSolver = NULL;
schurPrimalData.subSolver = NULL;
VecDestroy(&schurPrimalData.tmp_vec_b);
VecDestroy(&schurPrimalData.tmp_vec_primal);
MatDestroy(&mat_schur_primal);
KSPDestroy(&ksp_schur_primal);
} else {
KSPDestroy(&ksp_schur_primal);
}
MatDestroy(&mat_schur_primal);
KSPDestroy(&ksp_schur_primal);
}
......@@ -635,10 +627,8 @@ namespace AMDiS {
// === Create FETI-DP solver object. ===
fetiData.mat_primal_b = &(subDomainSolver->getMatCoarseInt());
fetiData.mat_b_primal = &(subDomainSolver->getMatIntCoarse());
fetiData.mat_lagrange = &mat_lagrange;
fetiData.fetiSolver = this;
fetiData.subSolver = subDomainSolver;
fetiData.ksp_schur_primal = &ksp_schur_primal;
VecCreateMPI(mpiComm,
......@@ -765,10 +755,8 @@ namespace AMDiS {
// === Destroy FETI-DP solver object. ===
fetiData.mat_primal_b = PETSC_NULL;
fetiData.mat_b_primal = PETSC_NULL;
fetiData.mat_lagrange = PETSC_NULL;
fetiData.fetiSolver = NULL;
fetiData.subSolver = NULL;
fetiData.ksp_schur_primal = PETSC_NULL;
VecDestroy(&fetiData.tmp_vec_b);
......@@ -933,26 +921,6 @@ namespace AMDiS {
int nRowsRankPrimal = primalDofMap.getRankDofs();
int nRowsOverallPrimal = primalDofMap.getOverallDofs();
#if 0
MatCreateSeqAIJ(PETSC_COMM_SELF, nRowsRankB, nRowsRankB, 60, PETSC_NULL,
&mat_b_b);
MatCreateMPIAIJ(mpiComm,
nRowsRankPrimal, nRowsRankPrimal,
nRowsOverallPrimal, nRowsOverallPrimal,
60, PETSC_NULL, 60, PETSC_NULL, &mat_primal_primal);
MatCreateMPIAIJ(mpiComm,
nRowsRankB, nRowsRankPrimal,
nRowsOverallB, nRowsOverallPrimal,
60, PETSC_NULL, 60, PETSC_NULL, &mat_b_primal);
MatCreateMPIAIJ(mpiComm,
nRowsRankPrimal, nRowsRankB,
nRowsOverallPrimal, nRowsOverallB,
30, PETSC_NULL, 30, PETSC_NULL, &mat_primal_b);
#endif
subDomainSolver->fillPetscMatrix(mat);
......@@ -982,88 +950,56 @@ namespace AMDiS {
}
}
// === Prepare traverse of sequentially created matrices. ===
using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits = mtl::traits;
typedef DOFMatrix::base_matrix_type Matrix;
typedef traits::range_generator<row, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
vector<int> cols, colsOther;
vector<double> values, valuesOther;
cols.reserve(300);
colsOther.reserve(300);
values.reserve(300);
valuesOther.reserve(300);
vector<int> colsLocal, colsLocalOther;
vector<double> valuesLocal, valuesLocalOther;
colsLocal.reserve(300);
colsLocalOther.reserve(300);
valuesLocal.reserve(300);
valuesLocalOther.reserve(300);
// === Traverse all sequentially created matrices and add the values to ===
// === the global PETSc matrices. ===
if (fetiPreconditioner != FETI_NONE) {
// === Prepare traverse of sequentially created matrices. ===
using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
namespace traits = mtl::traits;
typedef DOFMatrix::base_matrix_type Matrix;
typedef traits::range_generator<row, Matrix>::type cursor_type;
typedef traits::range_generator<nz, cursor_type>::type icursor_type;
vector<int> colsLocal, colsLocalOther;
vector<double> valuesLocal, valuesLocalOther;
colsLocal.reserve(300);
colsLocalOther.reserve(300);
valuesLocal.reserve(300);
valuesLocalOther.reserve(300);
// === Traverse all sequentially created matrices and add the values to ===
// === the global PETSc matrices. ===
int nComponents = mat->getSize();
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
if (!(*mat)[i][j])
continue;
int nComponents = mat->getSize();
for (int i = 0; i < nComponents; i++) {
for (int j = 0; j < nComponents; j++) {
if (!(*mat)[i][j])
continue;
traits::col<Matrix>::type col((*mat)[i][j]->getBaseMatrix());
traits::const_value<Matrix>::type value((*mat)[i][j]->getBaseMatrix());
traits::col<Matrix>::type col((*mat)[i][j]->getBaseMatrix());
traits::const_value<Matrix>::type value((*mat)[i][j]->getBaseMatrix());
// Traverse all rows.
for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()),
cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
// Traverse all rows.
for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()),
cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
bool rowPrimal = isPrimal(feSpaces[i], *cursor);
bool rowPrimal = isPrimal(feSpaces[i], *cursor);
cols.clear();
colsOther.clear();
values.clear();
valuesOther.clear();
colsLocal.clear();
colsLocalOther.clear();
valuesLocal.clear();
valuesLocalOther.clear();
colsLocal.clear();
colsLocalOther.clear();
valuesLocal.clear();
valuesLocalOther.clear();
// Traverse all columns.
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
bool colPrimal = isPrimal(feSpaces[j], col(*icursor));
if (colPrimal) {
if (rowPrimal) {
cols.push_back(col(*icursor));
values.push_back(value(*icursor));
} else {
colsOther.push_back(col(*icursor));
valuesOther.push_back(value(*icursor));
}
} else {
if (rowPrimal) {
colsOther.push_back(col(*icursor));
valuesOther.push_back(value(*icursor));
} else {
cols.push_back(col(*icursor));
values.push_back(value(*icursor));
}
}
// Traverse all columns.
for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor);
icursor != icend; ++icursor) {
// === For preconditioner ===
bool colPrimal = isPrimal(feSpaces[j], col(*icursor));
if (fetiPreconditioner != FETI_NONE) {
if (!rowPrimal && !colPrimal) {
if (!isDual(feSpaces[i], *cursor)) {
if (!isDual(feSpaces[j], col(*icursor))) {
......@@ -1087,129 +1023,65 @@ namespace AMDiS {
}
}
}
}
} // for each nnz in row
// === Set matrix values. ===
if (rowPrimal) {
int rowIndex = primalDofMap.getMatIndex(i, *cursor);
for (unsigned int k = 0; k < cols.size(); k++)
cols[k] = primalDofMap.getMatIndex(j, cols[k]);
} // for each nnz in row
#if 0
MatSetValues(mat_primal_primal, 1, &rowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES);
#endif
if (colsOther.size()) {
for (unsigned int k = 0; k < colsOther.size(); k++)
colsOther[k] = localDofMap.getMatIndex(j, colsOther[k]);
// === Set matrix values for preconditioner ===
#if 0
MatSetValues(mat_primal_b, 1, &rowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
#endif
}
} else {
int localRowIndex = localDofMap.getLocalMatIndex(i, *cursor);
for (unsigned int k = 0; k < cols.size(); k++)
cols[k] = localDofMap.getLocalMatIndex(j, cols[k]);
#if 0
MatSetValues(mat_b_b, 1, &localRowIndex, cols.size(),
&(cols[0]), &(values[0]), ADD_VALUES);
#endif
if (colsOther.size()) {
int globalRowIndex = localDofMap.getMatIndex(i, *cursor);
for (unsigned int k = 0; k < colsOther.size(); k++)
colsOther[k] = primalDofMap.getMatIndex(j, colsOther[k]);
#if 0
MatSetValues(mat_b_primal, 1, &globalRowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
#endif
if (!rowPrimal) {
switch (fetiPreconditioner) {
case FETI_DIRICHLET:
if (!isDual(feSpaces[i], *cursor)) {
int rowIndex = interiorDofMap.getLocalMatIndex(i, *cursor);
MatSetValues(mat_interior_interior, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_interior_duals, 1, &rowIndex, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
} else {
int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);
MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_duals_interior, 1, &rowIndex, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
}
break;
case FETI_LUMPED:
if (isDual(feSpaces[i], *cursor)) {
int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);
MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
}
break;
default:
break;
}
}
}
// === Set matrix values for preconditioner ===
if (!rowPrimal) {
switch (fetiPreconditioner) {
case FETI_DIRICHLET:
if (!isDual(feSpaces[i], *cursor)) {
int rowIndex = interiorDofMap.getLocalMatIndex(i, *cursor);
MatSetValues(mat_interior_interior, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_interior_duals, 1, &rowIndex, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
} else {
int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);
MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
if (colsLocalOther.size())
MatSetValues(mat_duals_interior, 1, &rowIndex, colsLocalOther.size(),
&(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
}
break;
case FETI_LUMPED:
if (isDual(feSpaces[i], *cursor)) {
int rowIndex = dualDofMap.getLocalMatIndex(i, *cursor);
MatSetValues(mat_duals_duals, 1, &rowIndex, colsLocal.size(),
&(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
}
break;
default:
break;
}
}
}
}
}
}
}
#if 0
// === Start global assembly procedure. ===
MatAssemblyBegin(mat_b_b, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_b_b, MAT_FINAL_ASSEMBLY);
// === Start global assembly procedure for preconditioner matrices. ===
MatAssemblyBegin(mat_primal_primal, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_primal_primal, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_b_primal, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_b_primal, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_primal_b, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_primal_b, MAT_FINAL_ASSEMBLY);
#endif
// === Start global assembly procedure for preconditioner matrices. ===
if (fetiPreconditioner != FETI_NONE) {
MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY);
}
if (fetiPreconditioner != FETI_NONE) {
MatAssemblyBegin(mat_duals_duals, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_duals_duals, MAT_FINAL_ASSEMBLY);
}
if (fetiPreconditioner == FETI_DIRICHLET) {
MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);
if (fetiPreconditioner == FETI_DIRICHLET) {
MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_interior_duals, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_interior_duals, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_interior_duals, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_interior_duals, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_duals_interior, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_duals_interior, MAT_FINAL_ASSEMBLY);
MatAssemblyBegin(mat_duals_interior, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(mat_duals_interior, MAT_FINAL_ASSEMBLY);
}
}
// === Create and fill PETSc matrix for Lagrange constraints. ===
createMatLagrange(feSpaces);
......@@ -1230,72 +1102,10 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::fillPetscRhs()");
#if 0
vector<const FiniteElemSpace*> feSpaces = getFeSpaces(vec);
VecCreateMPI(mpiComm,
localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &f_b);
VecCreateMPI(mpiComm,
primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(),
&f_primal);
for (unsigned int i = 0; i < feSpaces.size(); i++) {
DOFVector<double>::Iterator dofIt(vec->getDOFVector(i), USED_DOFS);
for (dofIt.reset(); !dofIt.end(); ++dofIt) {
int index = dofIt.getDOFIndex();
if (isPrimal(feSpaces[i], index)) {
index = primalDofMap.getMatIndex(i, index);
VecSetValue(f_primal, index, *dofIt, ADD_VALUES);
} else {
index = localDofMap.getMatIndex(i, index);
VecSetValue(f_b, index, *dofIt, INSERT_VALUES);
}
}
}
VecAssemblyBegin(f_b);
VecAssemblyEnd(f_b);
VecAssemblyBegin(f_primal);
VecAssemblyEnd(f_primal);
#else
subDomainSolver->fillPetscRhs(vec);
#endif
}
void PetscSolverFeti::solveLocalProblem(Vec &rhs, Vec &sol)
{
FUNCNAME("PetscSolverFeti::solveLocalProblem()");
Vec tmp;
VecCreateSeq(PETSC_COMM_SELF, localDofMap.getRankDofs(), &tmp);
PetscScalar *tmpValues, *rhsValues;
VecGetArray(tmp, &tmpValues);
VecGetArray(rhs, &rhsValues);
for (int i = 0; i < localDofMap.getRankDofs(); i++)
tmpValues[i] = rhsValues[i];
VecRestoreArray(rhs, &rhsValues);
VecRestoreArray(tmp, &tmpValues);
subDomainSolver->solve(tmp, tmp);
VecGetArray(tmp, &tmpValues);
VecGetArray(sol, &rhsValues);
for (int i = 0; i < localDofMap.getRankDofs(); i++)
rhsValues[i] = tmpValues[i];
VecRestoreArray(sol, &rhsValues);
VecRestoreArray(tmp, &tmpValues);
VecDestroy(&tmp);
}
void PetscSolverFeti::solveFetiMatrix(SystemVector &vec)
{
FUNCNAME("PetscSolverFeti::solveFetiMatrix()");
......@@ -1320,18 +1130,19 @@ namespace AMDiS {
localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &tmp_b0);
VecCreateMPI(mpiComm,
localDofMap.getRankDofs(), localDofMap.getOverallDofs(), &tmp_b1);
VecCreateMPI(mpiComm,
primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(), &tmp_primal0);
VecCreateMPI(mpiComm,
primalDofMap.getRankDofs(), primalDofMap.getOverallDofs(), &tmp_primal1);
MatGetVecs(mat_lagrange, PETSC_NULL, &tmp_lagrange0);
MatGetVecs(mat_lagrange, PETSC_NULL, &vec_rhs);
MatGetVecs(subDomainSolver->getMatCoarseCoarse(), PETSC_NULL, &tmp_primal0);
MatGetVecs(subDomainSolver->getMatCoarseCoarse(), PETSC_NULL, &tmp_primal1);
// === Create new rhs ===
// d = L inv(K_BB) f_B - L inv(K_BB) K_BPi inv(S_PiPi) [f_Pi - K_PiB inv(K_BB) f_B]
// vec_rhs = L * inv(K_BB) * f_B
solveLocalProblem(subDomainSolver->getRhsInterior(), tmp_b0);
subDomainSolver->solveGlobal(subDomainSolver->getRhsInterior(), tmp_b0);
MatMult(mat_lagrange, tmp_b0, vec_rhs);
// tmp_primal0 = M_PiB * inv(K_BB) * f_B
......@@ -1345,7 +1156,7 @@ namespace AMDiS {
//
MatMult(subDomainSolver->getMatIntCoarse(), tmp_primal0, tmp_b0);
solveLocalProblem(tmp_b0, tmp_b0);
subDomainSolver->solveGlobal(tmp_b0, tmp_b0);
MatMult(mat_lagrange, tmp_b0, tmp_lagrange0);
//
......@@ -1358,11 +1169,11 @@ namespace AMDiS {
// === Solve for u_primals. ===
VecCopy(subDomainSolver->getRhsCoarseSpace(), tmp_primal0);
solveLocalProblem(subDomainSolver->getRhsInterior(), tmp_b0);
subDomainSolver->solveGlobal(subDomainSolver->getRhsInterior(), tmp_b0);
MatMult(subDomainSolver->getMatCoarseInt(), tmp_b0, tmp_primal1);
VecAXPBY(tmp_primal0, -1.0, 1.0, tmp_primal1);
MatMultTranspose(mat_lagrange, vec_rhs, tmp_b0);
solveLocalProblem(tmp_b0, tmp_b0);
subDomainSolver->solveGlobal(tmp_b0, tmp_b0);
MatMult(subDomainSolver->getMatCoarseInt(), tmp_b0, tmp_primal1);
VecAXPBY(tmp_primal0, 1.0, 1.0, tmp_primal1);
......@@ -1377,28 +1188,18 @@ namespace AMDiS {
MatMult(subDomainSolver->getMatIntCoarse(), tmp_primal0, tmp_b1);
VecAXPBY(tmp_b0, -1.0, 1.0, tmp_b1);
solveLocalProblem(tmp_b0, tmp_b0);
subDomainSolver->solveGlobal(tmp_b0, tmp_b0);
// === And recover AMDiS solution vectors. ===
recoverSolution(tmp_b0, tmp_primal0, vec);
// === Destroy all data structures. ===
VecDestroy(&vec_rhs);
VecDestroy(&tmp_b0);
VecDestroy(&tmp_b1);
VecDestroy(&tmp_lagrange0);
VecDestroy(&tmp_primal0);
VecDestroy(&tmp_primal1);
#if 0
VecDestroy(&f_b);
VecDestroy(&f_primal);
#else
subDomainSolver->solvePetscMatrix(vec, NULL);
#endif
}
......@@ -1406,13 +1207,6 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::destroyMatrixData()");