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: macro __FUNCT__ to define routine names for use in error handling.\n\
5: Input parameters include:\n\
6: -user_defined_pc : Activate a user-defined preconditioner\n\n";
8: /*T
9: Concepts: KSP^basic parallel example
10: Concepts: PC^setting a user-defined shell preconditioner
11: Concepts: error handling^Using the macro __FUNCT__ to define routine names;
12: Processors: n
13: 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: we define the macro __FUNCT__ to be a string containing the routine name.
39: If defined, this macro is used in the PETSc error handlers to provide a
40: complete traceback of routine names. All PETSc library routines use this
41: macro, and users can optionally employ it as well in their application
42: codes. Note that users can get a traceback of PETSc errors regardless of
43: whether they define __FUNCT__ in application codes; this macro merely
44: provides the added traceback detail of the application routine names.
45: */
49: int main(int argc,char **args) 50: {
51: Vec x,b,u; /* approx solution, RHS, exact solution */
52: Mat A; /* linear system matrix */
53: KSP ksp; /* linear solver context */
54: PC pc; /* preconditioner context */
55: PetscReal norm; /* norm of solution error */
56: SampleShellPC *shell; /* user-defined preconditioner context */
57: PetscScalar v,one = 1.0,none = -1.0;
58: PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
60: PetscBool user_defined_pc = PETSC_FALSE;
62: PetscInitialize(&argc,&args,(char*)0,help);
63: PetscOptionsGetInt(NULL,"-m",&m,NULL);
64: PetscOptionsGetInt(NULL,"-n",&n,NULL);
66: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
67: Compute the matrix and right-hand-side vector that define
68: the linear system, Ax = b.
69: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70: /*
71: Create parallel matrix, specifying only its global dimensions.
72: When using MatCreate(), the matrix format can be specified at
73: runtime. Also, the parallel partioning of the matrix is
74: determined by PETSc at runtime.
75: */
76: MatCreate(PETSC_COMM_WORLD,&A);
77: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
78: MatSetFromOptions(A);
79: MatSetUp(A);
81: /*
82: Currently, all PETSc parallel matrix formats are partitioned by
83: contiguous chunks of rows across the processors. Determine which
84: rows of the matrix are locally owned.
85: */
86: MatGetOwnershipRange(A,&Istart,&Iend);
88: /*
89: Set matrix elements for the 2-D, five-point stencil in parallel.
90: - Each processor needs to insert only elements that it owns
91: locally (but any non-local elements will be sent to the
92: appropriate processor during matrix assembly).
93: - Always specify global rows and columns of matrix entries.
94: */
95: for (Ii=Istart; Ii<Iend; Ii++) {
96: v = -1.0; i = Ii/n; j = Ii - i*n;
97: if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
98: if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
99: if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
100: if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);}
101: v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);
102: }
104: /*
105: Assemble matrix, using the 2-step process:
106: MatAssemblyBegin(), MatAssemblyEnd()
107: Computations can be done while messages are in transition
108: by placing code between these two statements.
109: */
110: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
111: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
113: /*
114: Create parallel vectors.
115: - When using VecCreate() VecSetSizes() and VecSetFromOptions(),
116: we specify only the vector's global
117: dimension; the parallel partitioning is determined at runtime.
118: - Note: We form 1 vector from scratch and then duplicate as needed.
119: */
120: VecCreate(PETSC_COMM_WORLD,&u);
121: VecSetSizes(u,PETSC_DECIDE,m*n);
122: VecSetFromOptions(u);
123: VecDuplicate(u,&b);
124: VecDuplicate(b,&x);
126: /*
127: Set exact solution; then compute right-hand-side vector.
128: */
129: VecSet(u,one);
130: MatMult(A,u,b);
132: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: Create the linear solver and set various options
134: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: /*
137: Create linear solver context
138: */
139: KSPCreate(PETSC_COMM_WORLD,&ksp);
141: /*
142: Set operators. Here the matrix that defines the linear system
143: also serves as the preconditioning matrix.
144: */
145: KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
147: /*
148: Set linear solver defaults for this problem (optional).
149: - By extracting the KSP and PC contexts from the KSP context,
150: we can then directly call any KSP and PC routines
151: to set various options.
152: */
153: KSPGetPC(ksp,&pc);
154: KSPSetTolerances(ksp,1.e-7,PETSC_DEFAULT,PETSC_DEFAULT,
155: PETSC_DEFAULT);
157: /*
158: Set a user-defined "shell" preconditioner if desired
159: */
160: PetscOptionsGetBool(NULL,"-user_defined_pc",&user_defined_pc,NULL);
161: if (user_defined_pc) {
162: /* (Required) Indicate to PETSc that we're using a "shell" preconditioner */
163: PCSetType(pc,PCSHELL);
165: /* (Optional) Create a context for the user-defined preconditioner; this
166: context can be used to contain any application-specific data. */
167: SampleShellPCCreate(&shell);
169: /* (Required) Set the user-defined routine for applying the preconditioner */
170: PCShellSetApply(pc,SampleShellPCApply);
171: PCShellSetContext(pc,shell);
173: /* (Optional) Set user-defined function to free objects used by custom preconditioner */
174: PCShellSetDestroy(pc,SampleShellPCDestroy);
176: /* (Optional) Set a name for the preconditioner, used for PCView() */
177: PCShellSetName(pc,"MyPreconditioner");
179: /* (Optional) Do any setup required for the preconditioner */
180: /* Note: This function could be set with PCShellSetSetUp and it would be called when necessary */
181: SampleShellPCSetUp(pc,A,x);
183: } else {
184: PCSetType(pc,PCJACOBI);
185: }
187: /*
188: Set runtime options, e.g.,
189: -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
190: These options will override those specified above as long as
191: KSPSetFromOptions() is called _after_ any other customization
192: routines.
193: */
194: KSPSetFromOptions(ksp);
196: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197: Solve the linear system
198: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
200: KSPSolve(ksp,b,x);
202: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: Check solution and clean up
204: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206: /*
207: Check the error
208: */
209: VecAXPY(x,none,u);
210: VecNorm(x,NORM_2,&norm);
211: KSPGetIterationNumber(ksp,&its);
212: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G iterations %D\n",norm,its);
214: /*
215: Free work space. All PETSc objects should be destroyed when they
216: are no longer needed.
217: */
218: KSPDestroy(&ksp);
219: VecDestroy(&u); VecDestroy(&x);
220: VecDestroy(&b); MatDestroy(&A);
222: PetscFinalize();
223: return 0;
225: }
227: /***********************************************************************/
228: /* Routines for a user-defined shell preconditioner */
229: /***********************************************************************/
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;
245: PetscNew(SampleShellPC,&newctx);
246: newctx->diag = 0;
247: *shell = newctx;
248: return 0;
249: }
250: /* ------------------------------------------------------------------- */
253: /*
254: SampleShellPCSetUp - This routine sets up a user-defined
255: preconditioner context.
257: Input Parameters:
258: . pc - preconditioner object
259: . pmat - preconditioner matrix
260: . x - vector
262: Output Parameter:
263: . shell - fully set up user-defined preconditioner context
265: Notes:
266: In this example, we define the shell preconditioner to be Jacobi's
267: method. Thus, here we create a work vector for storing the reciprocal
268: of the diagonal of the preconditioner matrix; this vector is then
269: used within the routine SampleShellPCApply().
270: */
271: PetscErrorCode SampleShellPCSetUp(PC pc,Mat pmat,Vec x)272: {
273: SampleShellPC *shell;
274: Vec diag;
277: PCShellGetContext(pc,(void**)&shell);
278: VecDuplicate(x,&diag);
279: MatGetDiagonal(pmat,diag);
280: VecReciprocal(diag);
282: shell->diag = diag;
283: return 0;
284: }
285: /* ------------------------------------------------------------------- */
288: /*
289: SampleShellPCApply - This routine demonstrates the use of a
290: user-provided preconditioner.
292: Input Parameters:
293: + pc - preconditioner object
294: - x - input vector
296: Output Parameter:
297: . y - preconditioned vector
299: Notes:
300: This code implements the Jacobi preconditioner, merely as an
301: example of working with a PCSHELL. Note that the Jacobi method
302: is already provided within PETSc.
303: */
304: PetscErrorCode SampleShellPCApply(PC pc,Vec x,Vec y)305: {
306: SampleShellPC *shell;
309: PCShellGetContext(pc,(void**)&shell);
310: VecPointwiseMult(y,x,shell->diag);
312: return 0;
313: }
314: /* ------------------------------------------------------------------- */
317: /*
318: SampleShellPCDestroy - This routine destroys a user-defined
319: preconditioner context.
321: Input Parameter:
322: . shell - user-defined preconditioner context
323: */
324: PetscErrorCode SampleShellPCDestroy(PC pc)325: {
326: SampleShellPC *shell;
329: PCShellGetContext(pc,(void**)&shell);
330: VecDestroy(&shell->diag);
331: PetscFree(shell);
333: return 0;
334: }