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