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