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;
 47:   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;
 85:     i = Ii / n;
 86:     j = Ii - i * n;
 87:     if (i > 0) {
 88:       J = Ii - n;
 89:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 90:     }
 91:     if (i < m - 1) {
 92:       J = Ii + n;
 93:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 94:     }
 95:     if (j > 0) {
 96:       J = Ii - 1;
 97:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
 98:     }
 99:     if (j < n - 1) {
100:       J = Ii + 1;
101:       MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
102:     }
103:     v = 4.0;
104:     MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
105:   }

107:   /*
108:      Assemble matrix, using the 2-step process:
109:        MatAssemblyBegin(), MatAssemblyEnd()
110:      Computations can be done while messages are in transition
111:      by placing code between these two statements.
112:   */
113:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
114:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

116:   /*
117:      Create parallel vectors.
118:       - When using VecCreate() VecSetSizes() and VecSetFromOptions(),
119:         we specify only the vector's global
120:         dimension; the parallel partitioning is determined at runtime.
121:       - Note: We form 1 vector from scratch and then duplicate as needed.
122:   */
123:   VecCreate(PETSC_COMM_WORLD, &u);
124:   VecSetSizes(u, PETSC_DECIDE, m * n);
125:   VecSetFromOptions(u);
126:   VecDuplicate(u, &b);
127:   VecDuplicate(b, &x);

129:   /*
130:      Set exact solution; then compute right-hand-side vector.
131:   */
132:   VecSet(u, one);
133:   MatMult(A, u, b);

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:                 Create the linear solver and set various options
137:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

139:   /*
140:      Create linear solver context
141:   */
142:   KSPCreate(PETSC_COMM_WORLD, &ksp);

144:   /*
145:      Set operators. Here the matrix that defines the linear system
146:      also serves as the preconditioning matrix.
147:   */
148:   KSPSetOperators(ksp, A, A);

150:   /*
151:      Set linear solver defaults for this problem (optional).
152:      - By extracting the KSP and PC contexts from the KSP context,
153:        we can then directly call any KSP and PC routines
154:        to set various options.
155:   */
156:   KSPGetPC(ksp, &pc);
157:   KSPSetTolerances(ksp, 1.e-7, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);

159:   /*
160:      Set a user-defined "shell" preconditioner if desired
161:   */
162:   PetscOptionsGetBool(NULL, NULL, "-user_defined_pc", &user_defined_pc, NULL);
163:   if (user_defined_pc) {
164:     /* (Required) Indicate to PETSc that we're using a "shell" preconditioner */
165:     PCSetType(pc, PCSHELL);

167:     /* (Optional) Create a context for the user-defined preconditioner; this
168:        context can be used to contain any application-specific data. */
169:     SampleShellPCCreate(&shell);

171:     /* (Required) Set the user-defined routine for applying the preconditioner */
172:     PCShellSetApply(pc, SampleShellPCApply);
173:     PCShellSetContext(pc, shell);

175:     /* (Optional) Set user-defined function to free objects used by custom preconditioner */
176:     PCShellSetDestroy(pc, SampleShellPCDestroy);

178:     /* (Optional) Set a name for the preconditioner, used for PCView() */
179:     PCShellSetName(pc, "MyPreconditioner");

181:     /* (Optional) Do any setup required for the preconditioner */
182:     /* Note: This function could be set with PCShellSetSetUp and it would be called when necessary */
183:     SampleShellPCSetUp(pc, A, x);

185:   } else {
186:     PCSetType(pc, PCJACOBI);
187:   }

189:   /*
190:     Set runtime options, e.g.,
191:         -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
192:     These options will override those specified above as long as
193:     KSPSetFromOptions() is called _after_ any other customization
194:     routines.
195:   */
196:   KSPSetFromOptions(ksp);

198:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199:                       Solve the linear system
200:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

202:   KSPSolve(ksp, b, x);

204:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205:                       Check solution and clean up
206:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

208:   /*
209:      Check the error
210:   */
211:   VecAXPY(x, none, u);
212:   VecNorm(x, NORM_2, &norm);
213:   KSPGetIterationNumber(ksp, &its);
214:   PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g iterations %" PetscInt_FMT "\n", (double)norm, its);

216:   /*
217:      Free work space.  All PETSc objects should be destroyed when they
218:      are no longer needed.
219:   */
220:   KSPDestroy(&ksp);
221:   VecDestroy(&u);
222:   VecDestroy(&x);
223:   VecDestroy(&b);
224:   MatDestroy(&A);

226:   PetscFinalize();
227:   return 0;
228: }

230: /***********************************************************************/
231: /*          Routines for a user-defined shell preconditioner           */
232: /***********************************************************************/

234: /*
235:    SampleShellPCCreate - This routine creates a user-defined
236:    preconditioner context.

238:    Output Parameter:
239: .  shell - user-defined preconditioner context
240: */
241: PetscErrorCode SampleShellPCCreate(SampleShellPC **shell)
242: {
243:   SampleShellPC *newctx;

245:   PetscNew(&newctx);
246:   newctx->diag = 0;
247:   *shell       = newctx;
248:   return 0;
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:   PCShellGetContext(pc, &shell);
275:   VecDuplicate(x, &diag);
276:   MatGetDiagonal(pmat, diag);
277:   VecReciprocal(diag);

279:   shell->diag = diag;
280:   return 0;
281: }
282: /* ------------------------------------------------------------------- */
283: /*
284:    SampleShellPCApply - This routine demonstrates the use of a
285:    user-provided preconditioner.

287:    Input Parameters:
288: +  pc - preconditioner object
289: -  x - input vector

291:    Output Parameter:
292: .  y - preconditioned vector

294:    Notes:
295:    This code implements the Jacobi preconditioner, merely as an
296:    example of working with a PCSHELL.  Note that the Jacobi method
297:    is already provided within PETSc.
298: */
299: PetscErrorCode SampleShellPCApply(PC pc, Vec x, Vec y)
300: {
301:   SampleShellPC *shell;

303:   PCShellGetContext(pc, &shell);
304:   VecPointwiseMult(y, x, shell->diag);

306:   return 0;
307: }
308: /* ------------------------------------------------------------------- */
309: /*
310:    SampleShellPCDestroy - This routine destroys a user-defined
311:    preconditioner context.

313:    Input Parameter:
314: .  shell - user-defined preconditioner context
315: */
316: PetscErrorCode SampleShellPCDestroy(PC pc)
317: {
318:   SampleShellPC *shell;

320:   PCShellGetContext(pc, &shell);
321:   VecDestroy(&shell->diag);
322:   PetscFree(shell);

324:   return 0;
325: }

327: /*TEST

329:    build:
330:       requires: !complex !single

332:    test:
333:       nsize: 2
334:       args: -ksp_view -user_defined_pc -ksp_gmres_cgs_refinement_type refine_always

336:    test:
337:       suffix: tsirm
338:       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
339:       timeoutfactor: 4

341: TEST*/