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

Some performance issues in FETI-DP code.

parent e80c79ac
...@@ -27,7 +27,7 @@ namespace AMDiS { ...@@ -27,7 +27,7 @@ namespace AMDiS {
void *ctx; void *ctx;
MatShellGetContext(mat, &ctx); MatShellGetContext(mat, &ctx);
PetscSchurPrimalData* data = static_cast<PetscSchurPrimalData*>(ctx); SchurPrimalData* data = static_cast<SchurPrimalData*>(ctx);
MatMult(*(data->mat_b_primal), x, data->tmp_vec_b); MatMult(*(data->mat_b_primal), x, data->tmp_vec_b);
KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b); KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b);
...@@ -43,28 +43,25 @@ namespace AMDiS { ...@@ -43,28 +43,25 @@ namespace AMDiS {
// y = mat * x // y = mat * x
int petscMultMatFeti(Mat mat, Vec x, Vec y) int petscMultMatFeti(Mat mat, Vec x, Vec y)
{ {
// F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L) // F = L inv(K_BB) trans(L) + L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L)
// => F = L [I + inv(K_BB) K_BPi inv(S_PiPi) K_PiB] inv(K_BB) trans(L)
void *ctx; void *ctx;
MatShellGetContext(mat, &ctx); MatShellGetContext(mat, &ctx);
FetiData* data = static_cast<FetiData*>(ctx); FetiData* data = static_cast<FetiData*>(ctx);
// y = L inv(K_BB) trans(L) x // tmp_vec_b0 = inv(K_BB) trans(L) x
MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b); MatMultTranspose(*(data->mat_lagrange), x, data->tmp_vec_b0);
KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b); KSPSolve(*(data->ksp_b), data->tmp_vec_b0, data->tmp_vec_b0);
MatMult(*(data->mat_lagrange), data->tmp_vec_b, y);
// tmp_vec_primal = inv(S_PiPi) K_PiB inv(K_BB) trans(L) // tmp_vec_b1 = inv(K_BB) K_BPi inv(S_PiPi) K_PiB tmp_vec_b0
MatMult(*(data->mat_primal_b), data->tmp_vec_b, data->tmp_vec_primal); MatMult(*(data->mat_primal_b), data->tmp_vec_b0, data->tmp_vec_primal);
KSPSolve(*(data->ksp_schur_primal), data->tmp_vec_primal, 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_b1);
// tmp_vec_lagrange = L inv(K_BB) K_BPi tmp_vec_primal // y = L (tmp_vec_b0 + tmp_vec_b1)
// = L inv(K_BB) K_BPi inv(S_PiPi) K_PiB inv(K_BB) trans(L) VecAXPBY(data->tmp_vec_b0, 1.0, 1.0, data->tmp_vec_b1);
MatMult(*(data->mat_b_primal), data->tmp_vec_primal, data->tmp_vec_b); MatMult(*(data->mat_lagrange), data->tmp_vec_b0, y);
KSPSolve(*(data->ksp_b), data->tmp_vec_b, data->tmp_vec_b);
MatMult(*(data->mat_lagrange), data->tmp_vec_b, data->tmp_vec_lagrange);
VecAXPBY(y, 1.0, 1.0, data->tmp_vec_lagrange);
return 0; return 0;
} }
...@@ -593,18 +590,18 @@ namespace AMDiS { ...@@ -593,18 +590,18 @@ namespace AMDiS {
{ {
FUNCNAME("PetscSolverFeti::createSchurPrimal()"); FUNCNAME("PetscSolverFeti::createSchurPrimal()");
petscSchurPrimalData.mat_primal_primal = &mat_primal_primal; schurPrimalData.mat_primal_primal = &mat_primal_primal;
petscSchurPrimalData.mat_primal_b = &mat_primal_b; schurPrimalData.mat_primal_b = &mat_primal_b;
petscSchurPrimalData.mat_b_primal = &mat_b_primal; schurPrimalData.mat_b_primal = &mat_b_primal;
petscSchurPrimalData.ksp_b = &ksp_b; schurPrimalData.ksp_b = &ksp_b;
VecDuplicate(f_b, &(petscSchurPrimalData.tmp_vec_b)); VecDuplicate(f_b, &(schurPrimalData.tmp_vec_b));
VecDuplicate(f_primal, &(petscSchurPrimalData.tmp_vec_primal)); VecDuplicate(f_primal, &(schurPrimalData.tmp_vec_primal));
MatCreateShell(PETSC_COMM_WORLD, MatCreateShell(PETSC_COMM_WORLD,
nRankPrimals * nComponents, nRankPrimals * nComponents, nRankPrimals * nComponents, nRankPrimals * nComponents,
nOverallPrimals * nComponents, nOverallPrimals * nComponents, nOverallPrimals * nComponents, nOverallPrimals * nComponents,
&petscSchurPrimalData, &schurPrimalData,
&mat_schur_primal); &mat_schur_primal);
MatShellSetOperation(mat_schur_primal, MATOP_MULT, MatShellSetOperation(mat_schur_primal, MATOP_MULT,
(void(*)(void))petscMultMatSchurPrimal); (void(*)(void))petscMultMatSchurPrimal);
...@@ -620,13 +617,13 @@ namespace AMDiS { ...@@ -620,13 +617,13 @@ namespace AMDiS {
{ {
FUNCNAME("PetscSolverFeti::destroySchurPrimal()"); FUNCNAME("PetscSolverFeti::destroySchurPrimal()");
petscSchurPrimalData.mat_primal_primal = PETSC_NULL; schurPrimalData.mat_primal_primal = PETSC_NULL;
petscSchurPrimalData.mat_primal_b = PETSC_NULL; schurPrimalData.mat_primal_b = PETSC_NULL;
petscSchurPrimalData.mat_b_primal = PETSC_NULL; schurPrimalData.mat_b_primal = PETSC_NULL;
petscSchurPrimalData.ksp_b = PETSC_NULL; schurPrimalData.ksp_b = PETSC_NULL;
VecDestroy(&petscSchurPrimalData.tmp_vec_b); VecDestroy(&schurPrimalData.tmp_vec_b);
VecDestroy(&petscSchurPrimalData.tmp_vec_primal); VecDestroy(&schurPrimalData.tmp_vec_primal);
MatDestroy(&mat_schur_primal); MatDestroy(&mat_schur_primal);
KSPDestroy(&ksp_schur_primal); KSPDestroy(&ksp_schur_primal);
...@@ -646,11 +643,9 @@ namespace AMDiS { ...@@ -646,11 +643,9 @@ namespace AMDiS {
fetiData.ksp_b = &ksp_b; fetiData.ksp_b = &ksp_b;
fetiData.ksp_schur_primal = &ksp_schur_primal; fetiData.ksp_schur_primal = &ksp_schur_primal;
VecDuplicate(f_b, &(fetiData.tmp_vec_b0));
VecDuplicate(f_b, &(fetiData.tmp_vec_b)); VecDuplicate(f_b, &(fetiData.tmp_vec_b1));
VecDuplicate(f_primal, &(fetiData.tmp_vec_primal)); VecDuplicate(f_primal, &(fetiData.tmp_vec_primal));
MatGetVecs(mat_lagrange, PETSC_NULL, &(fetiData.tmp_vec_lagrange));
MatCreateShell(PETSC_COMM_WORLD, MatCreateShell(PETSC_COMM_WORLD,
nRankLagrange * nComponents, nRankLagrange * nComponents, nRankLagrange * nComponents, nRankLagrange * nComponents,
...@@ -729,10 +724,9 @@ namespace AMDiS { ...@@ -729,10 +724,9 @@ namespace AMDiS {
fetiData.ksp_b = PETSC_NULL; fetiData.ksp_b = PETSC_NULL;
fetiData.ksp_schur_primal = PETSC_NULL; fetiData.ksp_schur_primal = PETSC_NULL;
VecDestroy(&fetiData.tmp_vec_b); VecDestroy(&fetiData.tmp_vec_b0);
VecDestroy(&fetiData.tmp_vec_b1);
VecDestroy(&fetiData.tmp_vec_primal); VecDestroy(&fetiData.tmp_vec_primal);
VecDestroy(&fetiData.tmp_vec_lagrange);
MatDestroy(&mat_feti); MatDestroy(&mat_feti);
KSPDestroy(&ksp_feti); KSPDestroy(&ksp_feti);
...@@ -1391,17 +1385,19 @@ namespace AMDiS { ...@@ -1391,17 +1385,19 @@ namespace AMDiS {
// === Create new rhs === // === Create new rhs ===
// vec_rhs = L * inv(K_BB) * f_b // 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
KSPSolve(ksp_b, f_b, tmp_b0); KSPSolve(ksp_b, f_b, tmp_b0);
MatMult(mat_lagrange, tmp_b0, vec_rhs); MatMult(mat_lagrange, tmp_b0, vec_rhs);
// tmp_primal0 = M_PiB * inv(K_BB) * f_b // tmp_primal0 = M_PiB * inv(K_BB) * f_B
MatMult(mat_primal_b, tmp_b0, tmp_primal0); MatMult(mat_primal_b, tmp_b0, tmp_primal0);
// tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_b // tmp_primal0 = f_Pi - M_PiB * inv(K_BB) * f_B
VecAXPBY(tmp_primal0, -1.0, 1.0, f_primal); VecAXPBY(tmp_primal0, 1.0, -1.0, f_primal);
// tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_b) // tmp_primal0 = inv(S_PiPi) (f_Pi - M_PiB * inv(K_BB) * f_B)
KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0); KSPSolve(ksp_schur_primal, tmp_primal0, tmp_primal0);
// //
...@@ -1410,7 +1406,7 @@ namespace AMDiS { ...@@ -1410,7 +1406,7 @@ namespace AMDiS {
MatMult(mat_lagrange, tmp_b0, tmp_lagrange0); MatMult(mat_lagrange, tmp_b0, tmp_lagrange0);
// //
VecAXPBY(vec_rhs, 1.0, 1.0, tmp_lagrange0); VecAXPBY(vec_rhs, -1.0, 1.0, tmp_lagrange0);
// === Solve with FETI-DP operator. === // === Solve with FETI-DP operator. ===
......
...@@ -37,7 +37,7 @@ namespace AMDiS { ...@@ -37,7 +37,7 @@ namespace AMDiS {
* This structure is used when defining the MatShell operation for solving * This structure is used when defining the MatShell operation for solving
* primal schur complement. \ref petscMultMatSchurPrimal * primal schur complement. \ref petscMultMatSchurPrimal
*/ */
struct PetscSchurPrimalData { struct SchurPrimalData {
/// Pointers to the matrix containing the primal variables. /// Pointers to the matrix containing the primal variables.
Mat *mat_primal_primal; Mat *mat_primal_primal;
...@@ -70,15 +70,12 @@ namespace AMDiS { ...@@ -70,15 +70,12 @@ namespace AMDiS {
/// Matrix of Lagrange variables. /// Matrix of Lagrange variables.
Mat *mat_lagrange; Mat *mat_lagrange;
/// Temporal vector on the B variables. /// Temporal vectors on the B variables.
Vec tmp_vec_b; Vec tmp_vec_b0, tmp_vec_b1;
/// Temporal vector on the primal variables. /// Temporal vector on the primal variables.
Vec tmp_vec_primal; Vec tmp_vec_primal;
/// Temporal vector on the Lagrange variables.
Vec tmp_vec_lagrange;
/// Pointer to the solver for \ref PetscSolverFeti::mat_bb. /// Pointer to the solver for \ref PetscSolverFeti::mat_bb.
KSP *ksp_b; KSP *ksp_b;
...@@ -293,7 +290,7 @@ namespace AMDiS { ...@@ -293,7 +290,7 @@ namespace AMDiS {
Mat mat_schur_primal; Mat mat_schur_primal;
/// Data for MatMult operation in matrix \ref mat_schur_primal /// Data for MatMult operation in matrix \ref mat_schur_primal
PetscSchurPrimalData petscSchurPrimalData; SchurPrimalData schurPrimalData;
/// PETSc solver object to solve a system with FETI-DP. /// PETSc solver object to solve a system with FETI-DP.
KSP ksp_feti; KSP ksp_feti;
......
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