Actual source code: ex15.c

  1: static char help[] = "Solves a linear system in parallel with KSP.  Also\n\
  2: illustrates setting a user-defined shell preconditioner and using the\n\
  3: Input parameters include:\n\
  4:   -user_defined_pc : Activate a user-defined preconditioner\n\n";

  6: /*
  7:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  8:   automatically includes:
  9:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 10:      petscmat.h - matrices
 11:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 12:      petscviewer.h - viewers               petscpc.h  - preconditioners
 13: */
 14: #include <petscksp.h>

 16: /* Define context for user-provided preconditioner */
 17: typedef struct {
 18:   Vec diag;
 19: } SampleShellPC;

 21: /* Declare routines for user-provided preconditioner */
 22: extern PetscErrorCode SampleShellPCCreate(SampleShellPC **);
 23: extern PetscErrorCode SampleShellPCSetUp(PC, Mat, Vec);
 24: extern PetscErrorCode SampleShellPCApply(PC, Vec x, Vec y);
 25: extern PetscErrorCode SampleShellPCDestroy(PC);

 27: /*
 28:    User-defined routines.  Note that immediately before each routine below,
 29:    If defined, this macro is used in the PETSc error handlers to provide a
 30:    complete traceback of routine names.  All PETSc library routines use this
 31:    macro, and users can optionally employ it as well in their application
 32:    codes.  Note that users can get a traceback of PETSc errors regardless of
 33:    provides the added traceback detail of the application routine names.
 34: */

 36: int main(int argc, char **args)
 37: {
 38:   Vec            x, b, u; /* approx solution, RHS, exact solution */
 39:   Mat            A;       /* linear system matrix */
 40:   KSP            ksp;     /* linear solver context */
 41:   PC             pc;      /* preconditioner context */
 42:   PetscReal      norm;    /* norm of solution error */
 43:   SampleShellPC *shell;   /* user-defined preconditioner context */
 44:   PetscScalar    v, one = 1.0, none = -1.0;
 45:   PetscInt       i, j, Ii, J, Istart, Iend, m = 8, n = 7, its;
 46:   PetscBool      user_defined_pc = PETSC_FALSE;

 48:   PetscFunctionBeginUser;
 49:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
 50:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &m, NULL));
 51:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));

 53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 54:          Compute the matrix and right-hand-side vector that define
 55:          the linear system, Ax = b.
 56:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 57:   /*
 58:      Create parallel matrix, specifying only its global dimensions.
 59:      When using MatCreate(), the matrix format can be specified at
 60:      runtime. Also, the parallel partitioning of the matrix is
 61:      determined by PETSc at runtime.
 62:   */
 63:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 64:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m * n, m * n));
 65:   PetscCall(MatSetFromOptions(A));
 66:   PetscCall(MatSetUp(A));

 68:   /*
 69:      Currently, all PETSc parallel matrix formats are partitioned by
 70:      contiguous chunks of rows across the processors.  Determine which
 71:      rows of the matrix are locally owned.
 72:   */
 73:   PetscCall(MatGetOwnershipRange(A, &Istart, &Iend));

 75:   /*
 76:      Set matrix elements for the 2-D, five-point stencil in parallel.
 77:       - Each processor needs to insert only elements that it owns
 78:         locally (but any non-local elements will be sent to the
 79:         appropriate processor during matrix assembly).
 80:       - Always specify global rows and columns of matrix entries.
 81:    */
 82:   for (Ii = Istart; Ii < Iend; Ii++) {
 83:     v = -1.0;
 84:     i = Ii / n;
 85:     j = Ii - i * n;
 86:     if (i > 0) {
 87:       J = Ii - n;
 88:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 89:     }
 90:     if (i < m - 1) {
 91:       J = Ii + n;
 92:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 93:     }
 94:     if (j > 0) {
 95:       J = Ii - 1;
 96:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
 97:     }
 98:     if (j < n - 1) {
 99:       J = Ii + 1;
100:       PetscCall(MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES));
101:     }
102:     v = 4.0;
103:     PetscCall(MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES));
104:   }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

201:   PetscCall(KSPSolve(ksp, b, x));

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

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

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

225:   PetscCall(PetscFinalize());
226:   return 0;
227: }

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

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;

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

280:   shell->diag = diag;
281:   PetscFunctionReturn(PETSC_SUCCESS);
282: }
283: /* ------------------------------------------------------------------- */
284: /*
285:    SampleShellPCApply - This routine demonstrates the use of a
286:    user-provided preconditioner.

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

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

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

304:   PetscFunctionBeginUser;
305:   PetscCall(PCShellGetContext(pc, &shell));
306:   PetscCall(VecPointwiseMult(y, x, shell->diag));
307:   PetscFunctionReturn(PETSC_SUCCESS);
308: }
309: /* ------------------------------------------------------------------- */
310: /*
311:    SampleShellPCDestroy - This routine destroys a user-defined
312:    preconditioner context.

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

321:   PetscFunctionBeginUser;
322:   PetscCall(PCShellGetContext(pc, &shell));
323:   PetscCall(VecDestroy(&shell->diag));
324:   PetscCall(PetscFree(shell));
325:   PetscFunctionReturn(PETSC_SUCCESS);
326: }

328: /*TEST

330:    build:
331:       requires: !complex !single

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

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

342: TEST*/