Skip to content
GitLab
Commit bc23a38c authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Changed augmented coarse grid problem to its transpose.

parent 01656427
Hide whitespace changes
Inline Side-by-side
AMDiS/src/parallel/PetscSolverFeti.cc
View file @ bc23a38c
...@@ -683,13 +683,13 @@ namespace AMDiS { ...@@ -683,13 +683,13 @@ namespace AMDiS {
nOverallEdges *= feSpaces.size(); nOverallEdges *= feSpaces.size();
MatCreateAIJ(mpiCommGlobal, MatCreateAIJ(mpiCommGlobal,
lagrangeMap.getRankDofs(), nRankEdges, nRankEdges, lagrangeMap.getRankDofs(),
lagrangeMap.getOverallDofs(), nOverallEdges, nOverallEdges, lagrangeMap.getOverallDofs(),
1, PETSC_NULL, 1, PETSC_NULL, 1, PETSC_NULL, 1, PETSC_NULL,
&mat_augmented_lagrange); &mat_augmented_lagrange);
MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); MatSetOption(mat_augmented_lagrange, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
int colCounter = rStartEdges; int rowCounter = rStartEdges;
for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) { for (InteriorBoundary::iterator it(intBound.getOwn()); !it.end(); ++it) {
if (it->rankObj.subObj == EDGE) { if (it->rankObj.subObj == EDGE) {
for (int i = 0; i < feSpaces.size(); i++) { for (int i = 0; i < feSpaces.size(); i++) {
...@@ -703,12 +703,12 @@ namespace AMDiS { ...@@ -703,12 +703,12 @@ namespace AMDiS {
TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false) TEST_EXIT_DBG(isPrimal(feSpaces[i], **it) == false)
("Should not be primal!\n"); ("Should not be primal!\n");
int row = lagrangeMap.getMatIndex(i, **it); int col = lagrangeMap.getMatIndex(i, **it);
double value = 1.0; double value = 1.0;
MatSetValue(mat_augmented_lagrange, row, colCounter, value, INSERT_VALUES); MatSetValue(mat_augmented_lagrange, rowCounter, col, value, INSERT_VALUES);
} }
colCounter++; rowCounter++;
} }
} }
} }
...@@ -1928,10 +1928,14 @@ namespace AMDiS { ...@@ -1928,10 +1928,14 @@ namespace AMDiS {
MatMult(mat_lagrange, tmp_b0, tmp_lagrange); MatMult(mat_lagrange, tmp_b0, tmp_lagrange);
} else { } else {
Vec tmp_mu; Vec tmp_mu;
MatGetVecs(mat_augmented_lagrange, &tmp_mu, PETSC_NULL); MatGetVecs(mat_augmented_lagrange, PETSC_NULL, &tmp_mu);
MatMult(mat_augmented_lagrange, vecRhsLagrange, tmp_mu);
VecScale(tmp_mu, -1.0);
Vec vec_array[2] = {tmp_primal0, tmp_mu}; Vec vec_array[2] = {tmp_primal0, tmp_mu};
Vec vec_nest; Vec vec_nest;
VecCreateNest(PETSC_COMM_WORLD, 2, PETSC_NULL, vec_array, &vec_nest); VecCreateNest(PETSC_COMM_WORLD, 2, PETSC_NULL, vec_array, &vec_nest);
KSPSolve(ksp_schur_primal, vec_nest, vec_nest); KSPSolve(ksp_schur_primal, vec_nest, vec_nest);
...@@ -1943,7 +1947,7 @@ namespace AMDiS { ...@@ -1943,7 +1947,7 @@ namespace AMDiS {
// 2 Step // 2 Step
Vec tmp_lagrange1; Vec tmp_lagrange1;
VecDuplicate(tmp_lagrange, &tmp_lagrange1); VecDuplicate(tmp_lagrange, &tmp_lagrange1);
MatMult(mat_augmented_lagrange, tmp_mu, tmp_lagrange1); MatMultTranspose(mat_augmented_lagrange, tmp_mu, tmp_lagrange1);
MatMultTranspose(mat_lagrange, tmp_lagrange1, tmp_b0); MatMultTranspose(mat_lagrange, tmp_lagrange1, tmp_b0);
subdomain->solveGlobal(tmp_b0, tmp_b0); subdomain->solveGlobal(tmp_b0, tmp_b0);
MatMult(mat_lagrange, tmp_b0, tmp_lagrange1); MatMult(mat_lagrange, tmp_b0, tmp_lagrange1);
...@@ -1955,7 +1959,7 @@ namespace AMDiS { ...@@ -1955,7 +1959,7 @@ namespace AMDiS {
VecDestroy(&vec_nest); VecDestroy(&vec_nest);
} }
VecAXPBY(vecRhsLagrange, -1.0, 1.0, tmp_lagrange); VecAXPY(vecRhsLagrange, -1.0, tmp_lagrange);
} else { } else {
TEST_EXIT(augmentedLagrange)("Not yet implemented!\n"); TEST_EXIT(augmentedLagrange)("Not yet implemented!\n");
......
AMDiS/src/parallel/PetscSolverFetiOperators.cc
View file @ bc23a38c
...@@ -51,8 +51,8 @@ namespace AMDiS { ...@@ -51,8 +51,8 @@ namespace AMDiS {
MatMult(data->subSolver->getMatInteriorCoarse(), x_primal, data->tmp_vec_b0); MatMult(data->subSolver->getMatInteriorCoarse(), x_primal, data->tmp_vec_b0);
data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0); data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
// inv(K_BB) J^T Q x_mu // inv(K_BB) trans(J) trans(Q) x_mu
MatMult(*(data->mat_augmented_lagrange), x_mu, data->tmp_vec_lagrange); MatMultTranspose(*(data->mat_augmented_lagrange), x_mu, data->tmp_vec_lagrange);
MatMultTranspose(*(data->mat_lagrange), data->tmp_vec_lagrange, data->tmp_vec_b1); MatMultTranspose(*(data->mat_lagrange), data->tmp_vec_lagrange, data->tmp_vec_b1);
data->subSolver->solveGlobal(data->tmp_vec_b1, data->tmp_vec_b1); data->subSolver->solveGlobal(data->tmp_vec_b1, data->tmp_vec_b1);
...@@ -62,16 +62,16 @@ namespace AMDiS { ...@@ -62,16 +62,16 @@ namespace AMDiS {
MatMult(data->subSolver->getMatCoarse(), x_primal, y_primal); MatMult(data->subSolver->getMatCoarse(), x_primal, y_primal);
VecAXPY(y_primal, -1.0, data->tmp_vec_primal); VecAXPY(y_primal, -1.0, data->tmp_vec_primal);
// y_Pi += (-K_PiB inv(K_BB) J^T Q) x_mu // y_Pi += (-K_PiB inv(K_BB) J^T Q^T) x_mu
MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b1, MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b1,
data->tmp_vec_primal); data->tmp_vec_primal);
VecAXPY(y_primal, -1.0, data->tmp_vec_primal); VecAXPY(y_primal, -1.0, data->tmp_vec_primal);
// y_mu = (-Q^T J inv(K_BB) K_BPi) x_Pi + (-Q^T J inv(K_BB) J^T Q) x_mu // y_mu = (-Q J inv(K_BB) K_BPi) x_Pi + (-Q J inv(K_BB) J^T Q) x_mu
// = -Q^T J (inv(K_BB) K_BPi x_Pi + inv(K_BB) J^T Q x_mu) // = -Q J (inv(K_BB) K_BPi x_Pi + inv(K_BB) J^T Q x_mu)
VecAXPY(data->tmp_vec_b0, 1.0, data->tmp_vec_b1); VecAXPY(data->tmp_vec_b0, 1.0, data->tmp_vec_b1);
MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange); MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange);
MatMultTranspose(*(data->mat_augmented_lagrange), data->tmp_vec_lagrange, y_mu); MatMult(*(data->mat_augmented_lagrange), data->tmp_vec_lagrange, y_mu);
VecScale(y_mu, -1.0); VecScale(y_mu, -1.0);
return 0; return 0;
...@@ -135,7 +135,7 @@ namespace AMDiS { ...@@ -135,7 +135,7 @@ namespace AMDiS {
FetiData* data = static_cast<FetiData*>(ctx); FetiData* data = static_cast<FetiData*>(ctx);
Vec vec_mu0, vec_mu1; Vec vec_mu0, vec_mu1;
MatGetVecs(*(data->mat_augmented_lagrange), &vec_mu0, PETSC_NULL); MatGetVecs(*(data->mat_augmented_lagrange), PETSC_NULL, &vec_mu0);
VecDuplicate(vec_mu0, &vec_mu1); VecDuplicate(vec_mu0, &vec_mu1);
MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0); MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0);
...@@ -143,7 +143,7 @@ namespace AMDiS { ...@@ -143,7 +143,7 @@ namespace AMDiS {
MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b0, data->tmp_vec_primal0); MatMult(data->subSolver->getMatCoarseInterior(), data->tmp_vec_b0, data->tmp_vec_primal0);
MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange); MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange);
MatMultTranspose(*(data->mat_augmented_lagrange), data->tmp_vec_lagrange, vec_mu0); MatMult(*(data->mat_augmented_lagrange), data->tmp_vec_lagrange, vec_mu0);
Vec vec_array0[2] = {data->tmp_vec_primal0, vec_mu0}; Vec vec_array0[2] = {data->tmp_vec_primal0, vec_mu0};
Vec vec_array1[2] = {data->tmp_vec_primal1, vec_mu1}; Vec vec_array1[2] = {data->tmp_vec_primal1, vec_mu1};
...@@ -163,7 +163,7 @@ namespace AMDiS { ...@@ -163,7 +163,7 @@ namespace AMDiS {
VecAXPY(y, 1.0, data->tmp_vec_lagrange); VecAXPY(y, 1.0, data->tmp_vec_lagrange);
// Step 3 // Step 3
MatMult(*(data->mat_augmented_lagrange), vec_mu1, data->tmp_vec_lagrange); MatMultTranspose(*(data->mat_augmented_lagrange), vec_mu1, data->tmp_vec_lagrange);
MatMultTranspose(*(data->mat_lagrange), data->tmp_vec_lagrange, data->tmp_vec_b0); MatMultTranspose(*(data->mat_lagrange), data->tmp_vec_lagrange, data->tmp_vec_b0);
data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0); data->subSolver->solveGlobal(data->tmp_vec_b0, data->tmp_vec_b0);
MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange); MatMult(*(data->mat_lagrange), data->tmp_vec_b0, data->tmp_vec_lagrange);
......
Supports Markdown
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