Actual source code: ex23.c
1: const char help[] = "Test whether KSPSolve() can be executed without zeroing output for KSPPREONLY";
3: #include <petscksp.h>
5: static PetscErrorCode PCApply_scale(PC pc, Vec x, Vec y)
6: {
7: PetscFunctionBegin;
8: PetscCall(VecCopy(x, y));
9: PetscCall(VecScale(y, 0.5));
10: PetscFunctionReturn(PETSC_SUCCESS);
11: }
13: static PetscErrorCode VecSet_Error(Vec x, PetscReal _alpha)
14: {
15: PetscFunctionBegin;
16: SETERRQ(PetscObjectComm((PetscObject)x), PETSC_ERR_SUP, "Cannot zero this vector");
17: PetscFunctionReturn(PETSC_SUCCESS);
18: }
20: int main(int argc, char **argv)
21: {
22: Mat mat;
23: PetscInt M = 10;
24: MPI_Comm comm;
25: Vec x, y;
26: KSP ksp;
27: PC pc;
29: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
30: comm = PETSC_COMM_WORLD;
31: PetscCall(MatCreateConstantDiagonal(comm, PETSC_DETERMINE, PETSC_DETERMINE, M, M, 2.0, &mat));
32: PetscCall(MatCreateVecs(mat, &x, &y));
33: PetscCall(VecSet(x, 1.0));
34: PetscCall(VecSet(y, 3.0));
35: PetscCall(VecSetOperation(y, VECOP_SET, (void (*)(void))VecSet_Error));
37: PetscCall(KSPCreate(comm, &ksp));
38: PetscCall(KSPSetOperators(ksp, mat, mat));
39: PetscCall(KSPSetType(ksp, KSPPREONLY));
40: PetscCall(KSPGetPC(ksp, &pc));
41: PetscCall(PCSetType(pc, PCSHELL));
43: PetscCall(PCShellSetApply(pc, PCApply_scale));
44: PetscCall(KSPSolve(ksp, x, y));
46: PetscCall(KSPDestroy(&ksp));
47: PetscCall(VecDestroy(&y));
48: PetscCall(VecDestroy(&x));
49: PetscCall(MatDestroy(&mat));
50: PetscCall(PetscFinalize());
51: return 0;
52: }
54: /*TEST
56: test:
57: suffix: 0
59: TEST*/