diff --git a/AMDiS/src/parallel/PetscSolverFeti.cc b/AMDiS/src/parallel/PetscSolverFeti.cc
index e7241aaf32085cd31b04a5943eb119612df3ff60..8aaf8d05007f6b39e7ca22cbc63b80b61476221d 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.cc
+++ b/AMDiS/src/parallel/PetscSolverFeti.cc
@@ -70,6 +70,62 @@ namespace AMDiS {
}
+ // y = PC * x
+ PetscErrorCode petscApplyFetiPrecon(PC pc, Vec x, Vec y)
+ {
+ void *ctx;
+ PCShellGetContext(pc, &ctx);
+ PetscFetiPreconData* data = static_cast<PetscFetiPreconData*>(ctx);
+
+ MatMultTranspose(*(data->mat_lagrange_scaled), x, data->tmp_vec_b);
+
+ int sizeB;
+ int sizeBound;
+ VecGetLocalSize(data->tmp_vec_b, &sizeB);
+ VecGetLocalSize(data->tmp_vec_bound0, &sizeBound);
+
+ PetscScalar *local_b;
+ VecGetArray(data->tmp_vec_b, &local_b);
+
+ PetscScalar *local_bound;
+ VecGetArray(data->tmp_vec_bound0, &local_bound);
+
+ for (int i = sizeB - sizeBound, j = 0; i < sizeB; i++, j++)
+ local_bound[j] = local_b[i];
+
+ VecRestoreArray(data->tmp_vec_b, &local_b);
+ VecRestoreArray(data->tmp_vec_bound0, &local_bound);
+
+
+
+ MatMult(*(data->mat_bound_bound), data->tmp_vec_bound0, data->tmp_vec_bound1);
+
+ MatMult(*(data->mat_interior_bound), data->tmp_vec_bound0, data->tmp_vec_interior);
+ KSPSolve(*(data->ksp_interior), data->tmp_vec_interior, data->tmp_vec_interior);
+ MatMult(*(data->mat_bound_interior), data->tmp_vec_interior, data->tmp_vec_bound0);
+
+ VecAXPBY(data->tmp_vec_bound0, 1.0, -1.0, data->tmp_vec_bound1);
+
+
+
+ VecGetArray(data->tmp_vec_b, &local_b);
+ VecGetArray(data->tmp_vec_bound0, &local_bound);
+
+ for (int i = sizeB - sizeBound, j = 0; i < sizeB; i++, j++)
+ local_b[i] = local_bound[j];
+
+ VecRestoreArray(data->tmp_vec_b, &local_b);
+ VecRestoreArray(data->tmp_vec_bound0, &local_bound);
+
+
+
+
+ MatMult(*(data->mat_lagrange_scaled), data->tmp_vec_b, y);
+
+ return 0;
+ }
+
+
void PetscSolverFeti::updateDofData()
{
FUNCNAME("PetscSolverFeti::updateDofData()");
@@ -379,14 +435,15 @@ namespace AMDiS {
// === Get correct index for all interior nodes. ===
- int nInterior = 0;
+ nLocalInterior = 0;
for (DofMapping::iterator it = globalIndexB.begin();
it != globalIndexB.end(); ++it) {
- it->second = nInterior + rStartB;
- nInterior++;
+ it->second = nLocalInterior + rStartB;
+ nLocalInterior++;
}
+ nLocalBound = duals.size();
- TEST_EXIT_DBG(nInterior + primals.size() + duals.size() ==
+ TEST_EXIT_DBG(nLocalInterior + primals.size() + duals.size() ==
static_cast<unsigned int>(admin->getUsedDofs()))
("Should not happen!\n");
@@ -505,6 +562,8 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::createFetiKsp()");
+ // === Create FETI-DP solver object. ===
+
petscFetiData.mat_primal_primal = &mat_primal_primal;
petscFetiData.mat_primal_b = &mat_primal_b;
petscFetiData.mat_b_primal = &mat_b_primal;
@@ -529,6 +588,38 @@ namespace AMDiS {
KSPSetOperators(ksp_feti, mat_feti, mat_feti, SAME_NONZERO_PATTERN);
KSPSetOptionsPrefix(ksp_feti, "solver_feti_");
KSPSetFromOptions(ksp_feti);
+
+
+ // === Create FETI-DP Dirichlet preconditioner object. ===
+
+ KSPCreate(PETSC_COMM_SELF, &ksp_interior);
+ KSPSetOperators(ksp_interior, mat_interior_interior, mat_interior_interior, SAME_NONZERO_PATTERN);
+ KSPSetOptionsPrefix(ksp_interior, "solver_interior_");
+ KSPSetFromOptions(ksp_interior);
+
+
+ MatDuplicate(mat_lagrange, MAT_COPY_VALUES, &mat_lagrange_scaled);
+ MatScale(mat_lagrange_scaled, 0.5);
+
+ petscFetiPreconData.mat_lagrange_scaled = &mat_lagrange_scaled;
+ petscFetiPreconData.mat_interior_interior = &mat_interior_interior;
+ petscFetiPreconData.mat_bound_bound = &mat_bound_bound;
+ petscFetiPreconData.mat_interior_bound = &mat_interior_bound;
+ petscFetiPreconData.mat_bound_interior = &mat_bound_interior;
+ petscFetiPreconData.ksp_interior = &ksp_interior;
+ petscFetiPreconData.nInterior = nRankB - duals.size();
+
+ VecDuplicate(f_b, &(petscFetiPreconData.tmp_vec_b));
+
+ MatGetVecs(mat_bound_bound, PETSC_NULL, &(petscFetiPreconData.tmp_vec_bound0));
+ MatGetVecs(mat_bound_bound, PETSC_NULL, &(petscFetiPreconData.tmp_vec_bound1));
+ MatGetVecs(mat_interior_interior, PETSC_NULL, &(petscFetiPreconData.tmp_vec_interior));
+
+
+ KSPGetPC(ksp_feti, &precon_feti);
+ PCSetType(precon_feti, PCSHELL);
+ PCShellSetContext(precon_feti, static_cast<void*>(&petscFetiPreconData));
+ PCShellSetApply(precon_feti, petscApplyFetiPrecon);
}
@@ -536,6 +627,8 @@ namespace AMDiS {
{
FUNCNAME("PetscSolverFeti::destroyFetiKsp()");
+ // === Destroy FETI-DP solver object. ===
+
petscFetiData.mat_primal_primal = PETSC_NULL;
petscFetiData.mat_primal_b = PETSC_NULL;
petscFetiData.mat_b_primal = PETSC_NULL;
@@ -549,6 +642,24 @@ namespace AMDiS {
MatDestroy(mat_feti);
KSPDestroy(ksp_feti);
+
+
+ // === Destroy FETI-DP Dirichlet preconditioner object. ===
+
+ KSPDestroy(ksp_interior);
+
+ petscFetiPreconData.mat_lagrange_scaled = NULL;
+ petscFetiPreconData.mat_interior_interior = NULL;
+ petscFetiPreconData.mat_bound_bound = NULL;
+ petscFetiPreconData.mat_interior_bound = NULL;
+ petscFetiPreconData.mat_bound_interior = NULL;
+ petscFetiPreconData.ksp_interior = NULL;
+
+ VecDestroy(petscFetiPreconData.tmp_vec_b);
+ VecDestroy(petscFetiPreconData.tmp_vec_bound0);
+ VecDestroy(petscFetiPreconData.tmp_vec_bound1);
+ VecDestroy(petscFetiPreconData.tmp_vec_interior);
+ MatDestroy(mat_lagrange_scaled);
}
@@ -658,6 +769,8 @@ namespace AMDiS {
int nRowsOverallB = nOverallB * nComponents;
int nRowsRankPrimal = nRankPrimals * nComponents;
int nRowsOverallPrimal = nOverallPrimals * nComponents;
+ int nRowsInterior = nLocalInterior * nComponents;
+ int nRowsBound = nLocalBound * nComponents;
MatCreateMPIAIJ(PETSC_COMM_WORLD,
nRowsRankB, nRowsRankB, nRowsOverallB, nRowsOverallB,
@@ -678,6 +791,25 @@ namespace AMDiS {
nRowsOverallPrimal, nRowsOverallB,
100, PETSC_NULL, 100, PETSC_NULL, &mat_primal_b);
+
+ // === Create matrices for Dirichlet FETI-DP preconditioner. ===
+
+ MatCreateSeqAIJ(PETSC_COMM_SELF,
+ nRowsInterior, nRowsInterior, 100, PETSC_NULL,
+ &mat_interior_interior);
+
+ MatCreateSeqAIJ(PETSC_COMM_SELF,
+ nRowsBound, nRowsBound, 100, PETSC_NULL,
+ &mat_bound_bound);
+
+ MatCreateSeqAIJ(PETSC_COMM_SELF,
+ nRowsInterior, nRowsBound, 100, PETSC_NULL,
+ &mat_interior_bound);
+
+ MatCreateSeqAIJ(PETSC_COMM_SELF,
+ nRowsBound, nRowsInterior, 100, PETSC_NULL,
+ &mat_bound_interior);
+
// === Prepare traverse of sequentially created matrices. ===
@@ -695,6 +827,13 @@ namespace AMDiS {
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. ===
@@ -710,19 +849,27 @@ namespace AMDiS {
// 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 = primals.count(*cursor) != 0;
-
+
cols.clear();
- values.clear();
-
colsOther.clear();
+ values.clear();
valuesOther.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) {
- if (primals.count(col(*icursor)) != 0) {
+ bool colPrimal = primals.count(col(*icursor)) != 0;
+
+ if (colPrimal) {
// Column is a primal variable.
TEST_EXIT_DBG(globalPrimalIndex.count(col(*icursor)))
@@ -753,6 +900,47 @@ namespace AMDiS {
values.push_back(value(*icursor));
}
}
+
+
+
+ // === For preconditioner ===
+
+ if (!rowPrimal && !colPrimal) {
+ int rowIndex = globalIndexB[*cursor] - rStartB;
+ int colIndex = globalIndexB[col(*icursor)] - rStartB;
+
+ if (rowIndex < nLocalInterior) {
+ if (colIndex < nLocalInterior) {
+ int colIndex2 =
+ (globalIndexB[col(*icursor)] - rStartB) * nComponents + j;
+
+ colsLocal.push_back(colIndex2);
+ valuesLocal.push_back(value(*icursor));
+ } else {
+ int colIndex2 =
+ (globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;
+
+ colsLocalOther.push_back(colIndex2);
+ valuesLocalOther.push_back(value(*icursor));
+ }
+ } else {
+ if (colIndex < nLocalInterior) {
+ int colIndex2 =
+ (globalIndexB[col(*icursor)] - rStartB) * nComponents + j;
+
+ colsLocalOther.push_back(colIndex2);
+ valuesLocalOther.push_back(value(*icursor));
+ } else {
+ int colIndex2 =
+ (globalIndexB[col(*icursor)] - rStartB - nLocalInterior) * nComponents + j;
+
+ colsLocal.push_back(colIndex2);
+ valuesLocal.push_back(value(*icursor));
+ }
+ }
+ }
+
+
}
if (rowPrimal) {
@@ -778,6 +966,37 @@ namespace AMDiS {
MatSetValues(mat_b_primal, 1, &rowIndex, colsOther.size(),
&(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
}
+
+ // === For preconditioner ===
+
+ if (!rowPrimal) {
+ int rowIndex = globalIndexB[*cursor] - rStartB;
+
+ if (rowIndex < nLocalInterior) {
+ int rowIndex2 =
+ (globalIndexB[*cursor] - rStartB) * nComponents + i;
+
+ MatSetValues(mat_interior_interior, 1, &rowIndex2, colsLocal.size(),
+ &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
+
+ if (colsLocalOther.size())
+ MatSetValues(mat_interior_bound, 1, &rowIndex2, colsLocalOther.size(),
+ &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
+ } else {
+ int rowIndex2 =
+ (globalIndexB[*cursor] - rStartB - nLocalInterior) * nComponents + i;
+
+ MatSetValues(mat_bound_bound, 1, &rowIndex2, colsLocal.size(),
+ &(colsLocal[0]), &(valuesLocal[0]), INSERT_VALUES);
+
+ if (colsLocalOther.size())
+ MatSetValues(mat_bound_interior, 1, &rowIndex2, colsLocalOther.size(),
+ &(colsLocalOther[0]), &(valuesLocalOther[0]), INSERT_VALUES);
+
+ }
+ }
+
+
}
}
}
@@ -798,6 +1017,19 @@ namespace AMDiS {
MatAssemblyEnd(mat_primal_b, MAT_FINAL_ASSEMBLY);
+ MatAssemblyBegin(mat_interior_interior, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mat_interior_interior, MAT_FINAL_ASSEMBLY);
+
+ MatAssemblyBegin(mat_bound_bound, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mat_bound_bound, MAT_FINAL_ASSEMBLY);
+
+ MatAssemblyBegin(mat_interior_bound, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mat_interior_bound, MAT_FINAL_ASSEMBLY);
+
+ MatAssemblyBegin(mat_bound_interior, MAT_FINAL_ASSEMBLY);
+ MatAssemblyEnd(mat_bound_interior, MAT_FINAL_ASSEMBLY);
+
+
// === Create and fill PETSc's right hand side vectors. ===
VecCreate(PETSC_COMM_WORLD, &f_b);
@@ -1100,7 +1332,15 @@ namespace AMDiS {
destroySchurPrimalKsp();
- destroyFetiKsp();
+ destroyFetiKsp();
+
+
+ // === Destroy preconditioner data structures. ===
+
+ MatDestroy(mat_interior_interior);
+ MatDestroy(mat_bound_bound);
+ MatDestroy(mat_interior_bound);
+ MatDestroy(mat_bound_interior);
}
@@ -1117,9 +1357,9 @@ namespace AMDiS {
solveFetiMatrix(vec);
} else {
solveReducedFetiMatrix(vec);
- }
-
+ }
}
+
#endif
}
diff --git a/AMDiS/src/parallel/PetscSolverFeti.h b/AMDiS/src/parallel/PetscSolverFeti.h
index 5f6e96ff37be58ec3ac953ed99fc610dacedd049..d63f0031b22c22aaa12757e532813ac144bbbdaa 100644
--- a/AMDiS/src/parallel/PetscSolverFeti.h
+++ b/AMDiS/src/parallel/PetscSolverFeti.h
@@ -87,6 +87,26 @@ namespace AMDiS {
};
+ struct PetscFetiPreconData {
+ /// Matrix of scaled Lagrange variables.
+ Mat *mat_lagrange_scaled;
+
+ Mat *mat_interior_interior, *mat_bound_bound, *mat_interior_bound, *mat_bound_interior;
+
+ /// Pointer to the solver for \ref PetscSolverFeti::mat_bb.
+ KSP *ksp_interior;
+
+ /// Temporal vector on the B variables.
+ Vec tmp_vec_b;
+
+ Vec tmp_vec_bound0, tmp_vec_bound1;
+
+ Vec tmp_vec_interior;
+
+
+ int nInterior;
+ };
+
/** \brief
* FETI-DP implementation based on PETSc.
@@ -267,6 +287,21 @@ namespace AMDiS {
/// Data for MatMult operation in matrix \ref mat_feti
PetscFetiData petscFetiData;
+
+
+ Mat mat_lagrange_scaled;
+
+ PC precon_feti;
+
+ PetscFetiPreconData petscFetiPreconData;
+
+ Mat mat_interior_interior, mat_bound_bound, mat_interior_bound, mat_bound_interior;
+
+ KSP ksp_interior;
+
+ int nLocalInterior;
+
+ int nLocalBound;
};
#endif