Actual source code: ex55.c
1: static const char help[] = "Example demonstrating PCCOMPOSITE where one of the inner PCs uses a different operator\n\
2: \n";
4: #include <petscksp.h>
6: int main(int argc, char **argv)
7: {
8: PetscInt n = 10, i, col[3];
9: Vec x, b;
10: Mat A, B;
11: KSP ksp;
12: PC pc, subpc;
13: PetscScalar value[3];
15: PetscFunctionBeginUser;
16: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
18: /* Create a diagonal matrix with a given distribution of diagonal elements */
19: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
20: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
21: PetscCall(MatSetFromOptions(A));
22: PetscCall(MatSetUp(A));
23: /*
24: Assemble matrix
25: */
26: value[0] = -1.0;
27: value[1] = 2.0;
28: value[2] = -1.0;
29: for (i = 1; i < n - 1; i++) {
30: col[0] = i - 1;
31: col[1] = i;
32: col[2] = i + 1;
33: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
34: }
35: i = n - 1;
36: col[0] = n - 2;
37: col[1] = n - 1;
38: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
39: i = 0;
40: col[0] = 0;
41: col[1] = 1;
42: value[0] = 2.0;
43: value[1] = -1.0;
44: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
45: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
46: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
48: PetscCall(MatCreateVecs(A, &x, &b));
50: /* Create a KSP object */
51: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
52: PetscCall(KSPSetOperators(ksp, A, A));
54: /* Set up a composite preconditioner */
55: PetscCall(KSPGetPC(ksp, &pc));
56: PetscCall(PCSetType(pc, PCCOMPOSITE)); /* default composite with single Identity PC */
57: PetscCall(PCCompositeSetType(pc, PC_COMPOSITE_ADDITIVE));
58: PetscCall(PCCompositeAddPCType(pc, PCLU));
59: PetscCall(PCCompositeGetPC(pc, 0, &subpc));
60: /* B is set to the diagonal of A; this demonstrates that setting the operator for a subpc changes the preconditioning */
61: PetscCall(MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &B));
62: PetscCall(MatGetDiagonal(A, b));
63: PetscCall(MatDiagonalSet(B, b, ADD_VALUES));
64: PetscCall(PCSetOperators(subpc, B, B));
65: PetscCall(PCCompositeAddPCType(pc, PCNONE));
67: PetscCall(KSPSetFromOptions(ksp));
68: PetscCall(KSPSolve(ksp, b, x));
70: PetscCall(KSPDestroy(&ksp));
71: PetscCall(MatDestroy(&A));
72: PetscCall(MatDestroy(&B));
73: PetscCall(VecDestroy(&x));
74: PetscCall(VecDestroy(&b));
75: PetscCall(PetscFinalize());
76: return 0;
77: }
79: /*TEST
81: test:
82: args: -ksp_monitor -pc_composite_type multiplicative
84: TEST*/