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