Skip to content

Commit

Permalink
PCApplyTranspose_FieldSplit(): fix Schur complement implementation
Browse files Browse the repository at this point in the history 
Reported-by: Chris Douglas <[email protected]>

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Mat object's type is not set: Argument # 1
[0]PETSC ERROR: #1 MatGetFactorAvailable() at src/mat/interface/matrix.c:4820
[0]PETSC ERROR: #2 PCGetDefaultType_Private() at src/ksp/pc/interface/precon.c:25
[0]PETSC ERROR: #3 PCSetUp() at src/ksp/pc/interface/precon.c:1109
[0]PETSC ERROR: #4 KSPSetUp() at src/ksp/ksp/interface/itfunc.c:406
[0]PETSC ERROR: #5 KSPSolve_Private() at src/ksp/ksp/interface/itfunc.c:826
[0]PETSC ERROR: #6 KSPSolveTranspose() at src/ksp/ksp/interface/itfunc.c:1122
[0]PETSC ERROR: #7 PCApplyTranspose_FieldSplit() at src/ksp/pc/impls/fieldsplit/fieldsplit.c:1487
[0]PETSC ERROR: #8 PCApplyTranspose() at src/ksp/pc/interface/precon.c:715
[0]PETSC ERROR: #9 KSP_PCApply() at include/petsc/private/kspimpl.h:380
[0]PETSC ERROR: #10 KSPInitialResidual() at src/ksp/ksp/interface/itres.c:63
[0]PETSC ERROR: #11 KSPSolve_GMRES() at src/ksp/ksp/impls/gmres/gmres.c:226
[0]PETSC ERROR: #12 KSPSolve_Private() at src/ksp/ksp/interface/itfunc.c:900
[0]PETSC ERROR: #13 KSPSolveTranspose() at src/ksp/ksp/interface/itfunc.c:1122
[0]PETSC ERROR: #14 main() at ex9.c:80
[0]PETSC ERROR: Reached the main program with an out-of-range error code 1. This should never happen
  • Loading branch information
@prj-
prj- committed Jul 21, 2023
1 parent 6dde82a commit 7b66572
Show file tree
Hide file tree
Showing 4 changed files with 117 additions and 2 deletions.
1 change: 1 addition & 0 deletions src/ksp/ksp/tests/ex9.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,6 +77,7 @@ int main(int argc, char *argv[])
PetscCall(replace_submats(A));
PetscCall(replace_submats(P));
PetscCall(KSPSolve(ksp, b, x));
PetscCall(KSPSolveTranspose(ksp, b, x));

PetscCall(KSPDestroy(&ksp));
PetscCall(VecDestroy(&x));
Expand Down
1 change: 1 addition & 0 deletions src/ksp/ksp/tests/output/ex9_1.out
View file
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
Linear solve converged due to CONVERGED_ATOL iterations 1
Linear solve converged due to CONVERGED_ATOL iterations 1
Linear solve converged due to CONVERGED_ATOL iterations 1
1 change: 1 addition & 0 deletions src/ksp/ksp/tests/output/ex9_schur.out
View file
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
Linear solve converged due to CONVERGED_ATOL iterations 1
Linear solve converged due to CONVERGED_ATOL iterations 1
Linear solve converged due to CONVERGED_ATOL iterations 1
116 changes: 114 additions & 2 deletions src/ksp/pc/impls/fieldsplit/fieldsplit.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1204,6 +1204,117 @@ static PetscErrorCode PCApply_FieldSplit_Schur(PC pc, Vec x, Vec y)
PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApplyTranspose_FieldSplit_Schur(PC pc, Vec x, Vec y)
{
PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
PC_FieldSplitLink ilinkA = jac->head, ilinkD = ilinkA->next;
KSP kspA = ilinkA->ksp, kspLower = kspA, kspUpper = jac->kspupper;

PetscFunctionBegin;
switch (jac->schurfactorization) {
case PC_FIELDSPLIT_SCHUR_FACT_DIAG:
/* [A00 0; 0 -S], positive definite, suitable for MINRES */
PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
PetscCall(VecScale(ilinkD->y, jac->schurscale));
PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
break;
case PC_FIELDSPLIT_SCHUR_FACT_UPPER:
PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
PetscCall(VecScale(ilinkD->x, -1.));
PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
break;
case PC_FIELDSPLIT_SCHUR_FACT_LOWER:
PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
PetscCall(VecScale(ilinkA->x, -1.));
PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, ADD_VALUES, SCATTER_FORWARD));
PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
break;
case PC_FIELDSPLIT_SCHUR_FACT_FULL:
PetscCall(VecScatterBegin(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(VecScatterEnd(ilinkA->sctx, x, ilinkA->x, INSERT_VALUES, SCATTER_FORWARD));
PetscCall(PetscLogEventBegin(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
PetscCall(KSPSolveTranspose(kspUpper, ilinkA->x, ilinkA->y));
PetscCall(KSPCheckSolve(kspUpper, pc, ilinkA->y));
PetscCall(PetscLogEventEnd(KSP_Solve_FS_U, kspUpper, ilinkA->x, ilinkA->y, NULL));
PetscCall(MatMultTranspose(jac->B, ilinkA->y, ilinkD->x));
PetscCall(VecScale(ilinkD->x, -1.0));
PetscCall(VecScatterBegin(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));
PetscCall(VecScatterEnd(ilinkD->sctx, x, ilinkD->x, ADD_VALUES, SCATTER_FORWARD));

PetscCall(PetscLogEventBegin(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
PetscCall(KSPSolveTranspose(jac->kspschur, ilinkD->x, ilinkD->y));
PetscCall(KSPCheckSolve(jac->kspschur, pc, ilinkD->y));
PetscCall(PetscLogEventEnd(KSP_Solve_FS_S, jac->kspschur, ilinkD->x, ilinkD->y, NULL));
PetscCall(VecScatterBegin(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));

if (kspLower == kspA) {
PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->y));
PetscCall(VecAXPY(ilinkA->x, -1.0, ilinkA->y));
PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
PetscCall(PetscLogEventEnd(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
} else {
PetscCall(PetscLogEventBegin(ilinkA->event, kspA, ilinkA->x, ilinkA->y, NULL));
PetscCall(KSPSolveTranspose(kspA, ilinkA->x, ilinkA->y));
PetscCall(KSPCheckSolve(kspA, pc, ilinkA->y));
PetscCall(MatMultTranspose(jac->C, ilinkD->y, ilinkA->x));
PetscCall(PetscLogEventBegin(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
PetscCall(KSPSolveTranspose(kspLower, ilinkA->x, ilinkA->z));
PetscCall(KSPCheckSolve(kspLower, pc, ilinkA->z));
PetscCall(PetscLogEventEnd(KSP_Solve_FS_L, kspLower, ilinkA->x, ilinkA->z, NULL));
PetscCall(VecAXPY(ilinkA->y, -1.0, ilinkA->z));
}
PetscCall(VecScatterEnd(ilinkD->sctx, ilinkD->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterBegin(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
PetscCall(VecScatterEnd(ilinkA->sctx, ilinkA->y, y, INSERT_VALUES, SCATTER_REVERSE));
}
PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode PCApply_FieldSplit(PC pc, Vec x, Vec y)
{
PC_FieldSplit *jac = (PC_FieldSplit *)pc->data;
Expand Down Expand Up @@ -2751,8 +2862,9 @@ static PetscErrorCode PCFieldSplitSetType_FieldSplit(PC pc, PCCompositeType type
PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetGKBDelay_C", NULL));

if (type == PC_COMPOSITE_SCHUR) {
pc->ops->apply = PCApply_FieldSplit_Schur;
pc->ops->view = PCView_FieldSplit_Schur;
pc->ops->apply = PCApply_FieldSplit_Schur;
pc->ops->applytranspose = PCApplyTranspose_FieldSplit_Schur;
pc->ops->view = PCView_FieldSplit_Schur;

PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitGetSubKSP_C", PCFieldSplitGetSubKSP_FieldSplit_Schur));
PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCFieldSplitSetSchurPre_C", PCFieldSplitSetSchurPre_FieldSplit));
Expand Down

0 comments on commit 7b66572

Please sign in to comment.