Actual source code: ex15.c
1: static char help[] = "Solves a linear system in parallel with KSP. Also\n\
2: illustrates setting a user-defined shell preconditioner and using the\n\
3: Input parameters include:\n\
4: -user_defined_pc : Activate a user-defined preconditioner\n\n";
6: /*
7: Include "petscksp.h" so that we can use KSP solvers. Note that this file
8: automatically includes:
9: petscsys.h - base PETSc routines petscvec.h - vectors
10: petscmat.h - matrices
11: petscis.h - index sets petscksp.h - Krylov subspace methods
12: petscviewer.h - viewers petscpc.h - preconditioners
13: */
14: #include <petscksp.h>
16: /* Define context for user-provided preconditioner */
17: typedef struct {
18: Vec diag;
19: } SampleShellPC;
21: /* Declare routines for user-provided preconditioner */
22: extern PetscErrorCode SampleShellPCCreate(SampleShellPC **);
23: extern PetscErrorCode SampleShellPCSetUp(PC, Mat, Vec);
24: extern PetscErrorCode SampleShellPCApply(PC, Vec x, Vec y);
25: extern PetscErrorCode SampleShellPCDestroy(PC);
27: /*
28: User-defined routines. Note that immediately before each routine below,
29: If defined, this macro is used in the PETSc error handlers to provide a
30: complete traceback of routine names. All PETSc library routines use this
31: macro, and users can optionally employ it as well in their application
32: codes. Note that users can get a traceback of PETSc errors regardless of
33: provides the added traceback detail of the application routine names.
34: */
36: int main(int argc, char **args)
37: {
38: Vec x, b, u; /* approx solution, RHS, exact solution */
39: Mat A; /* linear system matrix */
40: KSP ksp; /* linear solver context */
41: PC pc; /* preconditioner context */
42: PetscReal norm; /* norm of solution error */
43: SampleShellPC *shell; /* user-defined preconditioner context */
44: PetscScalar v, one = 1.0, none = -1.0;
45: PetscInt i, j, Ii, J, Istart, Iend, m = 8, n = 7, its;
46: PetscBool user_defined_pc = PETSC_FALSE;
48: PetscFunctionBeginUser;
49: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
50: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
51: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
53: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: Compute the matrix and right-hand-side vector that define
55: the linear system, Ax = b.
56: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
57: /*
58: Create parallel matrix, specifying only its global dimensions.
59: When using MatCreate(), the matrix format can be specified at
60: runtime. Also, the parallel partitioning of the matrix is
61: determined by PETSc at runtime.
62: */
63: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
64: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
65: PetscCall(MatSetFromOptions(A));
66: PetscCall(MatSetUp(A));
68: /*
69: Currently, all PETSc parallel matrix formats are partitioned by
70: contiguous chunks of rows across the processors. Determine which
71: rows of the matrix are locally owned.
72: */
73: PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));
75: /*
76: Set matrix elements for the 2-D, five-point stencil in parallel.
77: - Each processor needs to insert only elements that it owns
78: locally (but any non-local elements will be sent to the
79: appropriate processor during matrix assembly).
80: - Always specify global rows and columns of matrix entries.
81: */
82: for (Ii = Istart; Ii < Iend; Ii++) {
83: v = -1.0;
84: i = Ii / n;
85: j = Ii - i * n;
86: if (i > 0) {
87: J = Ii - n;
88: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
89: }
90: if (i < m - 1) {
91: J = Ii + n;
92: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
93: }
94: if (j > 0) {
95: J = Ii - 1;
96: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
97: }
98: if (j < n - 1) {
99: J = Ii + 1;
100: PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
101: }
102: v = 4.0;
103: PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
104: }
106: /*
107: Assemble matrix, using the 2-step process:
108: MatAssemblyBegin(), MatAssemblyEnd()
109: Computations can be done while messages are in transition
110: by placing code between these two statements.
111: */
112: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
113: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
115: /*
116: Create parallel vectors.
117: - When using VecCreate() VecSetSizes() and VecSetFromOptions(),
118: we specify only the vector's global
119: dimension; the parallel partitioning is determined at runtime.
120: - Note: We form 1 vector from scratch and then duplicate as needed.
121: */
122: PetscCall(VecCreate(PETSC_COMM_WORLD, &u));
123: PetscCall(VecSetSizes(u, PETSC_DECIDE, m * n));
124: PetscCall(VecSetFromOptions(u));
125: PetscCall(VecDuplicate(u, &b));
126: PetscCall(VecDuplicate(b, &x));
128: /*
129: Set exact solution; then compute right-hand-side vector.
130: */
131: PetscCall(VecSet(u, one));
132: PetscCall(MatMult(A, u, b));
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Create the linear solver and set various options
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: /*
139: Create linear solver context
140: */
141: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
143: /*
144: Set operators. Here the matrix that defines the linear system
145: also serves as the preconditioning matrix.
146: */
147: PetscCall(KSPSetOperators(ksp, A, A));
149: /*
150: Set linear solver defaults for this problem (optional).
151: - By extracting the KSP and PC contexts from the KSP context,
152: we can then directly call any KSP and PC routines
153: to set various options.
154: */
155: PetscCall(KSPGetPC(ksp, &pc));
156: PetscCall(KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT));
158: /*
159: Set a user-defined "shell" preconditioner if desired
160: */
161: PetscCall(PetscOptionsGetBool(NULL, NULL, "-user_defined_pc", &user_defined_pc, NULL));
162: if (user_defined_pc) {
163: /* (Required) Indicate to PETSc that we're using a "shell" preconditioner */
164: PetscCall(PCSetType(pc, PCSHELL));
166: /* (Optional) Create a context for the user-defined preconditioner; this
167: context can be used to contain any application-specific data. */
168: PetscCall(SampleShellPCCreate(&shell));
170: /* (Required) Set the user-defined routine for applying the preconditioner */
171: PetscCall(PCShellSetApply(pc, SampleShellPCApply));
172: PetscCall(PCShellSetContext(pc, shell));
174: /* (Optional) Set user-defined function to free objects used by custom preconditioner */
175: PetscCall(PCShellSetDestroy(pc, SampleShellPCDestroy));
177: /* (Optional) Set a name for the preconditioner, used for PCView() */
178: PetscCall(PCShellSetName(pc, "MyPreconditioner"));
180: /* (Optional) Do any setup required for the preconditioner */
181: /* Note: This function could be set with PCShellSetSetUp and it would be called when necessary */
182: PetscCall(SampleShellPCSetUp(pc, A, x));
184: } else {
185: PetscCall(PCSetType(pc, PCJACOBI));
186: }
188: /*
189: Set runtime options, e.g.,
190: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
191: These options will override those specified above as long as
192: KSPSetFromOptions() is called _after_ any other customization
193: routines.
194: */
195: PetscCall(KSPSetFromOptions(ksp));
197: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: Solve the linear system
199: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201: PetscCall(KSPSolve(ksp, b, x));
203: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: Check solution and clean up
205: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
207: /*
208: Check the error
209: */
210: PetscCall(VecAXPY(x, none, u));
211: PetscCall(VecNorm(x, NORM_2, &norm));
212: PetscCall(KSPGetIterationNumber(ksp, &its));
213: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its));
215: /*
216: Free work space. All PETSc objects should be destroyed when they
217: are no longer needed.
218: */
219: PetscCall(KSPDestroy(&ksp));
220: PetscCall(VecDestroy(&u));
221: PetscCall(VecDestroy(&x));
222: PetscCall(VecDestroy(&b));
223: PetscCall(MatDestroy(&A));
225: PetscCall(PetscFinalize());
226: return 0;
227: }
229: /***********************************************************************/
230: /* Routines for a user-defined shell preconditioner */
231: /***********************************************************************/
233: /*
234: SampleShellPCCreate - This routine creates a user-defined
235: preconditioner context.
237: Output Parameter:
238: . shell - user-defined preconditioner context
239: */
240: PetscErrorCode SampleShellPCCreate(SampleShellPC **shell)
241: {
242: SampleShellPC *newctx;
244: PetscFunctionBeginUser;
245: PetscCall(PetscNew(&newctx));
246: newctx->diag = 0;
247: *shell = newctx;
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
250: /* ------------------------------------------------------------------- */
251: /*
252: SampleShellPCSetUp - This routine sets up a user-defined
253: preconditioner context.
255: Input Parameters:
256: . pc - preconditioner object
257: . pmat - preconditioner matrix
258: . x - vector
260: Output Parameter:
261: . shell - fully set up user-defined preconditioner context
263: Notes:
264: In this example, we define the shell preconditioner to be Jacobi's
265: method. Thus, here we create a work vector for storing the reciprocal
266: of the diagonal of the preconditioner matrix; this vector is then
267: used within the routine SampleShellPCApply().
268: */
269: PetscErrorCode SampleShellPCSetUp(PC pc, Mat pmat, Vec x)
270: {
271: SampleShellPC *shell;
272: Vec diag;
274: PetscFunctionBeginUser;
275: PetscCall(PCShellGetContext(pc, &shell));
276: PetscCall(VecDuplicate(x, &diag));
277: PetscCall(MatGetDiagonal(pmat, diag));
278: PetscCall(VecReciprocal(diag));
280: shell->diag = diag;
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
283: /* ------------------------------------------------------------------- */
284: /*
285: SampleShellPCApply - This routine demonstrates the use of a
286: user-provided preconditioner.
288: Input Parameters:
289: + pc - preconditioner object
290: - x - input vector
292: Output Parameter:
293: . y - preconditioned vector
295: Notes:
296: This code implements the Jacobi preconditioner, merely as an
297: example of working with a PCSHELL. Note that the Jacobi method
298: is already provided within PETSc.
299: */
300: PetscErrorCode SampleShellPCApply(PC pc, Vec x, Vec y)
301: {
302: SampleShellPC *shell;
304: PetscFunctionBeginUser;
305: PetscCall(PCShellGetContext(pc, &shell));
306: PetscCall(VecPointwiseMult(y, x, shell->diag));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
309: /* ------------------------------------------------------------------- */
310: /*
311: SampleShellPCDestroy - This routine destroys a user-defined
312: preconditioner context.
314: Input Parameter:
315: . shell - user-defined preconditioner context
316: */
317: PetscErrorCode SampleShellPCDestroy(PC pc)
318: {
319: SampleShellPC *shell;
321: PetscFunctionBeginUser;
322: PetscCall(PCShellGetContext(pc, &shell));
323: PetscCall(VecDestroy(&shell->diag));
324: PetscCall(PetscFree(shell));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: /*TEST
330: build:
331: requires: !complex !single
333: test:
334: nsize: 2
335: args: -ksp_view -user_defined_pc -ksp_gmres_cgs_refinement_type refine_always
337: test:
338: suffix: tsirm
339: args: -m 60 -n 60 -ksp_type tsirm -pc_type ksp -ksp_monitor_short -ksp_ksp_type fgmres -ksp_ksp_rtol 1e-10 -ksp_pc_type mg -ksp_ksp_max_it 30
340: timeoutfactor: 4
342: TEST*/