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