Actual source code: ex28.c
2: static char help[] = "Test procedural KSPSetFromOptions() or at runtime; Test PCREDUNDANT.\n\n";
4: #include <petscksp.h>
6: int main(int argc, char **args)
7: {
8: Vec x, b, u; /* approx solution, RHS, exact solution */
9: Mat A; /* linear system matrix */
10: KSP ksp; /* linear solver context */
11: PC pc; /* preconditioner context */
12: PetscReal norm; /* norm of solution error */
13: PetscInt i, n = 10, col[3], its, rstart, rend, nlocal;
14: PetscScalar one = 1.0, value[3];
15: PetscBool TEST_PROCEDURAL = PETSC_FALSE;
18: PetscInitialize(&argc, &args, (char *)0, help);
19: PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL);
20: 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: VecCreate(PETSC_COMM_WORLD, &x);
34: VecSetSizes(x, PETSC_DECIDE, n);
35: VecSetFromOptions(x);
36: VecDuplicate(x, &b);
37: 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: VecGetOwnershipRange(x, &rstart, &rend);
44: VecGetLocalSize(x, &nlocal);
46: /* Create a tridiagonal matrix. See ../tutorials/ex23.c */
47: MatCreate(PETSC_COMM_WORLD, &A);
48: MatSetSizes(A, nlocal, nlocal, n, n);
49: MatSetFromOptions(A);
50: 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: 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: 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: MatSetValues(A, 1, &i, 3, col, value, INSERT_VALUES);
80: }
81: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
82: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
84: /* Set exact solution; then compute right-hand-side vector. */
85: VecSet(u, one);
86: MatMult(A, u, b);
88: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
89: Create the linear solver and set various options
90: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
91: KSPCreate(PETSC_COMM_WORLD, &ksp);
92: 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: PetscMPIInt size, rank, subsize;
106: Mat A_redundant;
107: KSP innerksp;
108: PC innerpc;
109: MPI_Comm subcomm;
111: KSPGetPC(ksp, &pc);
112: PCSetType(pc, PCREDUNDANT);
113: MPI_Comm_size(PETSC_COMM_WORLD, &size);
114: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
116: PCRedundantSetNumber(pc, size - 2);
117: KSPSetFromOptions(ksp);
119: /* Get subcommunicator and redundant matrix */
120: KSPSetUp(ksp);
121: PCRedundantGetKSP(pc, &innerksp);
122: KSPGetPC(innerksp, &innerpc);
123: PCGetOperators(innerpc, NULL, &A_redundant);
124: PetscObjectGetComm((PetscObject)A_redundant, &subcomm);
125: MPI_Comm_size(subcomm, &subsize);
126: if (subsize == 1 && rank == 0) {
127: PetscPrintf(PETSC_COMM_SELF, "A_redundant:\n");
128: MatView(A_redundant, PETSC_VIEWER_STDOUT_SELF);
129: }
130: } else {
131: KSPSetFromOptions(ksp);
132: }
134: /* Solve linear system */
135: KSPSolve(ksp, b, x);
137: /* Check the error */
138: VecAXPY(x, -1.0, u);
139: VecNorm(x, NORM_2, &norm);
140: KSPGetIterationNumber(ksp, &its);
141: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g, Iterations %" PetscInt_FMT "\n", (double)norm, its);
143: /* Free work space. */
144: VecDestroy(&x);
145: VecDestroy(&u);
146: VecDestroy(&b);
147: MatDestroy(&A);
148: KSPDestroy(&ksp);
149: PetscFinalize();
150: return 0;
151: }
153: /*TEST
155: test:
156: nsize: 3
157: output_file: output/ex28.out
159: test:
160: suffix: 2
161: args: -procedural -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi
162: nsize: 3
164: test:
165: suffix: 3
166: args: -procedural -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type bjacobi
167: nsize: 5
169: TEST*/