Actual source code: ex28.c
1: static char help[] = "Test procedural KSPSetFromOptions() or at runtime; Test PCREDUNDANT.\n\n";
3: #include <petscksp.h>
5: int main(int argc, char **args)
6: {
7: Vec x, b, u; /* approx solution, RHS, exact solution */
8: Mat A; /* linear system matrix */
9: KSP ksp; /* linear solver context */
10: PC pc; /* preconditioner context */
11: PetscReal norm; /* norm of solution error */
12: PetscInt i, n = 10, col[3], its, rstart, rend, nlocal;
13: PetscMPIInt size, rank, subsize;
14: PetscScalar one = 1.0, value[3];
15: PetscBool TEST_PROCEDURAL = PETSC_FALSE;
17: PetscFunctionBeginUser;
18: PetscCall(PetscInitialize(&argc, &args, NULL, help));
19: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
20: PetscCall(PetscOptionsGetBool(NULL, NULL, "-procedural", &TEST_PROCEDURAL, NULL));
22: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
23: Compute the matrix and right-hand-side vector that define
24: the linear system, Ax = b.
25: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
27: /*
28: Create vectors. Note that we form 1 vector from scratch and
29: then duplicate as needed. For this simple case let PETSc decide how
30: many elements of the vector are stored on each processor. The second
31: argument to VecSetSizes() below causes PETSc to decide.
32: */
33: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
34: PetscCall(VecSetSizes(x, PETSC_DECIDE, n));
35: PetscCall(VecSetFromOptions(x));
36: PetscCall(VecDuplicate(x, &b));
37: PetscCall(VecDuplicate(x, &u));
39: /* Identify the starting and ending mesh points on each
40: processor for the interior part of the mesh. We let PETSc decide
41: above. */
43: PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
44: PetscCall(VecGetLocalSize(x, &nlocal));
46: /* Create a tridiagonal matrix. See ../tutorials/ex23.c */
47: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
48: PetscCall(MatSetSizes(A, nlocal, nlocal, n, n));
49: PetscCall(MatSetFromOptions(A));
50: PetscCall(MatSetUp(A));
51: /* Assemble matrix */
52: if (!rstart) {
53: rstart = 1;
54: i = 0;
55: col[0] = 0;
56: col[1] = 1;
57: value[0] = 2.0;
58: value[1] = -1.0;
59: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
60: }
61: if (rend == n) {
62: rend = n - 1;
63: i = n - 1;
64: col[0] = n - 2;
65: col[1] = n - 1;
66: value[0] = -1.0;
67: value[1] = 2.0;
68: PetscCall(MatSetValues(A, 1, &i, 2, col, value, INSERT_VALUES));
69: }
71: /* Set entries corresponding to the mesh interior */
72: value[0] = -1.0;
73: value[1] = 2.0;
74: value[2] = -1.0;
75: for (i = rstart; i < rend; i++) {
76: col[0] = i - 1;
77: col[1] = i;
78: col[2] = i + 1;
79: PetscCall(MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES));
80: }
81: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
82: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
84: /* Set exact solution; then compute right-hand-side vector. */
85: PetscCall(VecSet(u, one));
86: PetscCall(MatMult(A, u, b));
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Create the linear solver and set various options
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
92: PetscCall(KSPSetOperators(ksp, A, A));
94: /*
95: Set linear solver defaults for this problem (optional).
96: - By extracting the KSP and PC contexts from the KSP context,
97: we can then directly call any KSP and PC routines to set
98: various options.
99: - The following statements are optional; all of these
100: parameters could alternatively be specified at runtime via
101: KSPSetFromOptions();
102: */
103: if (TEST_PROCEDURAL) {
104: /* Example of runtime options: '-pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi' */
105: PetscCall(KSPGetPC(ksp, &pc));
106: PetscCall(PCSetType(pc, PCREDUNDANT));
107: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
108: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
109: PetscCheck(size > 2, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, "Num of processes %d must greater than 2", size);
110: PetscCall(PCRedundantSetNumber(pc, size - 2));
111: }
112: PetscCall(KSPSetFromOptions(ksp));
113: if (TEST_PROCEDURAL) {
114: PetscBool flg;
116: PetscCall(KSPSetUp(ksp));
117: PetscCall(KSPGetPC(ksp, &pc));
118: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCREDUNDANT, &flg));
119: if (flg) {
120: KSP innerksp;
121: PC innerpc;
123: PetscCall(PCRedundantGetKSP(pc, &innerksp));
124: PetscCall(KSPSetFromOptions(innerksp));
125: PetscCall(KSPGetPC(innerksp, &innerpc));
126: PetscCall(PetscObjectTypeCompare((PetscObject)innerpc, PCILU, &flg));
127: if (flg) PetscCall(PCFactorSetLevels(innerpc, 1));
128: }
129: }
131: /* Solve linear system */
132: PetscCall(KSPSolve(ksp, b, x));
134: if (TEST_PROCEDURAL) {
135: Mat A_redundant;
136: KSP innerksp;
137: PC innerpc;
138: MPI_Comm subcomm;
140: /* Get subcommunicator and redundant matrix */
141: PetscCall(PCRedundantGetKSP(pc, &innerksp));
142: PetscCall(KSPGetPC(innerksp, &innerpc));
143: PetscCall(PCGetOperators(innerpc, NULL, &A_redundant));
144: PetscCall(PetscObjectGetComm((PetscObject)A_redundant, &subcomm));
145: PetscCallMPI(MPI_Comm_size(subcomm, &subsize));
146: if (subsize == 1 && rank == 0) {
147: PetscCall(PetscPrintf(PETSC_COMM_SELF, "A_redundant:\n"));
148: PetscCall(MatView(A_redundant, PETSC_VIEWER_STDOUT_SELF));
149: }
150: }
152: /* Check the error */
153: PetscCall(VecAXPY(x, -1.0, u));
154: PetscCall(VecNorm(x, NORM_2, &norm));
155: PetscCall(KSPGetIterationNumber(ksp, &its));
156: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its));
158: /* Free work space. */
159: PetscCall(VecDestroy(&x));
160: PetscCall(VecDestroy(&u));
161: PetscCall(VecDestroy(&b));
162: PetscCall(MatDestroy(&A));
163: PetscCall(KSPDestroy(&ksp));
164: PetscCall(PetscFinalize());
165: return 0;
166: }
168: /*TEST
170: test:
171: nsize: 3
172: output_file: output/ex28.out
174: test:
175: suffix: 2
176: args: -procedural -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi
177: nsize: 3
179: test:
180: suffix: 3
181: args: -procedural -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi
182: nsize: 5
184: test:
185: suffix: 4
186: requires: defined(PETSC_USE_LOG)
187: args: -procedural -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type ilu -log_view
188: filter: grep -o "MatLUFactorNum 1"
189: nsize: 3
191: TEST*/