Actual source code: ex3.c
1: static char help[] = "Newton methods to solve u'' + u^{2} = f in parallel.\n\
2: This example employs a user-defined monitoring routine and optionally a user-defined\n\
3: routine to check candidate iterates produced by line search routines.\n\
4: The command line options include:\n\
5: -pre_check_iterates : activate checking of iterates\n\
6: -post_check_iterates : activate checking of iterates\n\
7: -check_tol <tol>: set tolerance for iterate checking\n\
8: -user_precond : activate a (trivial) user-defined preconditioner\n\n";
10: /*
11: Include "petscdm.h" so that we can use data management objects (DMs)
12: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
13: Include "petscsnes.h" so that we can use SNES solvers. Note that this
14: file automatically includes:
15: petscsys.h - base PETSc routines
16: petscvec.h - vectors
17: petscmat.h - matrices
18: petscis.h - index sets
19: petscksp.h - Krylov subspace methods
20: petscviewer.h - viewers
21: petscpc.h - preconditioners
22: petscksp.h - linear solvers
23: */
25: #include <petscdm.h>
26: #include <petscdmda.h>
27: #include <petscsnes.h>
29: /*
30: User-defined routines.
31: */
32: PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
33: PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
34: PetscErrorCode FormInitialGuess(Vec);
35: PetscErrorCode Monitor(SNES, PetscInt, PetscReal, void *);
36: PetscErrorCode PreCheck(SNESLineSearch, Vec, Vec, PetscBool *, void *);
37: PetscErrorCode PostCheck(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
38: PetscErrorCode PostSetSubKSP(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *);
39: PetscErrorCode MatrixFreePreconditioner(PC, Vec, Vec);
41: /*
42: User-defined application context
43: */
44: typedef struct {
45: DM da; /* distributed array */
46: Vec F; /* right-hand side of PDE */
47: PetscMPIInt rank; /* rank of processor */
48: PetscMPIInt size; /* size of communicator */
49: PetscReal h; /* mesh spacing */
50: PetscBool sjerr; /* if or not to test jacobian domain error */
51: } ApplicationCtx;
53: /*
54: User-defined context for monitoring
55: */
56: typedef struct {
57: PetscViewer viewer;
58: } MonitorCtx;
60: /*
61: User-defined context for checking candidate iterates that are
62: determined by line search methods
63: */
64: typedef struct {
65: Vec last_step; /* previous iterate */
66: PetscReal tolerance; /* tolerance for changes between successive iterates */
67: ApplicationCtx *user;
68: } StepCheckCtx;
70: typedef struct {
71: PetscInt its0; /* num of previous outer KSP iterations */
72: } SetSubKSPCtx;
74: int main(int argc, char **argv)
75: {
76: SNES snes; /* SNES context */
77: SNESLineSearch linesearch; /* SNESLineSearch context */
78: Mat J; /* Jacobian matrix */
79: ApplicationCtx ctx; /* user-defined context */
80: Vec x, r, U, F; /* vectors */
81: MonitorCtx monP; /* monitoring context */
82: StepCheckCtx checkP; /* step-checking context */
83: SetSubKSPCtx checkP1;
84: PetscBool pre_check, post_check, post_setsubksp; /* flag indicating whether we're checking candidate iterates */
85: PetscScalar xp, *FF, *UU, none = -1.0;
86: PetscInt its, N = 5, i, maxit, maxf, xs, xm;
87: PetscReal abstol, rtol, stol, norm;
88: PetscBool flg, viewinitial = PETSC_FALSE;
90: PetscFunctionBeginUser;
91: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
92: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank));
93: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size));
94: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL));
95: ctx.h = 1.0 / (N - 1);
96: ctx.sjerr = PETSC_FALSE;
97: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test_jacobian_domain_error", &ctx.sjerr, NULL));
98: PetscCall(PetscOptionsGetBool(NULL, NULL, "-view_initial", &viewinitial, NULL));
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create nonlinear solver context
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create vector data structures; set function evaluation routine
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: /*
111: Create distributed array (DMDA) to manage parallel grid and vectors
112: */
113: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, 1, 1, NULL, &ctx.da));
114: PetscCall(DMSetFromOptions(ctx.da));
115: PetscCall(DMSetUp(ctx.da));
117: /*
118: Extract global and local vectors from DMDA; then duplicate for remaining
119: vectors that are the same types
120: */
121: PetscCall(DMCreateGlobalVector(ctx.da, &x));
122: PetscCall(VecDuplicate(x, &r));
123: PetscCall(VecDuplicate(x, &F));
124: ctx.F = F;
125: PetscCall(VecDuplicate(x, &U));
127: /*
128: Set function evaluation routine and vector. Whenever the nonlinear
129: solver needs to compute the nonlinear function, it will call this
130: routine.
131: - Note that the final routine argument is the user-defined
132: context that provides application-specific data for the
133: function evaluation routine.
134: */
135: PetscCall(SNESSetFunction(snes, r, FormFunction, &ctx));
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Create matrix data structure; set Jacobian evaluation routine
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
142: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, N, N));
143: PetscCall(MatSetFromOptions(J));
144: PetscCall(MatSeqAIJSetPreallocation(J, 3, NULL));
145: PetscCall(MatMPIAIJSetPreallocation(J, 3, NULL, 3, NULL));
147: /*
148: Set Jacobian matrix data structure and default Jacobian evaluation
149: routine. Whenever the nonlinear solver needs to compute the
150: Jacobian matrix, it will call this routine.
151: - Note that the final routine argument is the user-defined
152: context that provides application-specific data for the
153: Jacobian evaluation routine.
154: */
155: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, &ctx));
157: /*
158: Optionally allow user-provided preconditioner
159: */
160: PetscCall(PetscOptionsHasName(NULL, NULL, "-user_precond", &flg));
161: if (flg) {
162: KSP ksp;
163: PC pc;
164: PetscCall(SNESGetKSP(snes, &ksp));
165: PetscCall(KSPGetPC(ksp, &pc));
166: PetscCall(PCSetType(pc, PCSHELL));
167: PetscCall(PCShellSetApply(pc, MatrixFreePreconditioner));
168: }
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Customize nonlinear solver; set runtime options
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: /*
175: Set an optional user-defined monitoring routine
176: */
177: PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 400, 400, &monP.viewer));
178: PetscCall(SNESMonitorSet(snes, Monitor, &monP, 0));
180: /*
181: Set names for some vectors to facilitate monitoring (optional)
182: */
183: PetscCall(PetscObjectSetName((PetscObject)x, "Approximate Solution"));
184: PetscCall(PetscObjectSetName((PetscObject)U, "Exact Solution"));
186: /*
187: Set SNES/KSP/KSP/PC runtime options, e.g.,
188: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
189: */
190: PetscCall(SNESSetFromOptions(snes));
192: /*
193: Set an optional user-defined routine to check the validity of candidate
194: iterates that are determined by line search methods
195: */
196: PetscCall(SNESGetLineSearch(snes, &linesearch));
197: PetscCall(PetscOptionsHasName(NULL, NULL, "-post_check_iterates", &post_check));
199: if (post_check) {
200: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activating post step checking routine\n"));
201: PetscCall(SNESLineSearchSetPostCheck(linesearch, PostCheck, &checkP));
202: PetscCall(VecDuplicate(x, &checkP.last_step));
204: checkP.tolerance = 1.0;
205: checkP.user = &ctx;
207: PetscCall(PetscOptionsGetReal(NULL, NULL, "-check_tol", &checkP.tolerance, NULL));
208: }
210: PetscCall(PetscOptionsHasName(NULL, NULL, "-post_setsubksp", &post_setsubksp));
211: if (post_setsubksp) {
212: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activating post setsubksp\n"));
213: PetscCall(SNESLineSearchSetPostCheck(linesearch, PostSetSubKSP, &checkP1));
214: }
216: PetscCall(PetscOptionsHasName(NULL, NULL, "-pre_check_iterates", &pre_check));
217: if (pre_check) {
218: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Activating pre step checking routine\n"));
219: PetscCall(SNESLineSearchSetPreCheck(linesearch, PreCheck, &checkP));
220: }
222: /*
223: Print parameters used for convergence testing (optional) ... just
224: to demonstrate this routine; this information is also printed with
225: the option -snes_view
226: */
227: PetscCall(SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf));
228: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "atol=%g, rtol=%g, stol=%g, maxit=%" PetscInt_FMT ", maxf=%" PetscInt_FMT "\n", (double)abstol, (double)rtol, (double)stol, maxit, maxf));
230: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
231: Initialize application:
232: Store right-hand side of PDE and exact solution
233: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235: /*
236: Get local grid boundaries (for 1-dimensional DMDA):
237: xs, xm - starting grid index, width of local grid (no ghost points)
238: */
239: PetscCall(DMDAGetCorners(ctx.da, &xs, NULL, NULL, &xm, NULL, NULL));
241: /*
242: Get pointers to vector data
243: */
244: PetscCall(DMDAVecGetArray(ctx.da, F, &FF));
245: PetscCall(DMDAVecGetArray(ctx.da, U, &UU));
247: /*
248: Compute local vector entries
249: */
250: xp = ctx.h * xs;
251: for (i = xs; i < xs + xm; i++) {
252: FF[i] = 6.0 * xp + PetscPowScalar(xp + 1.e-12, 6.0); /* +1.e-12 is to prevent 0^6 */
253: UU[i] = xp * xp * xp;
254: xp += ctx.h;
255: }
257: /*
258: Restore vectors
259: */
260: PetscCall(DMDAVecRestoreArray(ctx.da, F, &FF));
261: PetscCall(DMDAVecRestoreArray(ctx.da, U, &UU));
262: if (viewinitial) {
263: PetscCall(VecView(U, 0));
264: PetscCall(VecView(F, 0));
265: }
267: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268: Evaluate initial guess; then solve nonlinear system
269: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271: /*
272: Note: The user should initialize the vector, x, with the initial guess
273: for the nonlinear solver prior to calling SNESSolve(). In particular,
274: to employ an initial guess of zero, the user should explicitly set
275: this vector to zero by calling VecSet().
276: */
277: PetscCall(FormInitialGuess(x));
278: PetscCall(SNESSolve(snes, NULL, x));
279: PetscCall(SNESGetIterationNumber(snes, &its));
280: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %" PetscInt_FMT "\n", its));
282: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
283: Check solution and clean up
284: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
286: /*
287: Check the error
288: */
289: PetscCall(VecAXPY(x, none, U));
290: PetscCall(VecNorm(x, NORM_2, &norm));
291: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)norm, its));
292: if (ctx.sjerr) {
293: SNESType snestype;
294: PetscCall(SNESGetType(snes, &snestype));
295: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES Type %s\n", snestype));
296: }
298: /*
299: Free work space. All PETSc objects should be destroyed when they
300: are no longer needed.
301: */
302: PetscCall(PetscViewerDestroy(&monP.viewer));
303: if (post_check) PetscCall(VecDestroy(&checkP.last_step));
304: PetscCall(VecDestroy(&x));
305: PetscCall(VecDestroy(&r));
306: PetscCall(VecDestroy(&U));
307: PetscCall(VecDestroy(&F));
308: PetscCall(MatDestroy(&J));
309: PetscCall(SNESDestroy(&snes));
310: PetscCall(DMDestroy(&ctx.da));
311: PetscCall(PetscFinalize());
312: return 0;
313: }
315: /* ------------------------------------------------------------------- */
316: /*
317: FormInitialGuess - Computes initial guess.
319: Input/Output Parameter:
320: . x - the solution vector
321: */
322: PetscErrorCode FormInitialGuess(Vec x)
323: {
324: PetscScalar pfive = .50;
326: PetscFunctionBeginUser;
327: PetscCall(VecSet(x, pfive));
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: /* ------------------------------------------------------------------- */
332: /*
333: FormFunction - Evaluates nonlinear function, F(x).
335: Input Parameters:
336: . snes - the SNES context
337: . x - input vector
338: . ctx - optional user-defined context, as set by SNESSetFunction()
340: Output Parameter:
341: . f - function vector
343: Note:
344: The user-defined context can contain any application-specific
345: data needed for the function evaluation.
346: */
347: PetscErrorCode FormFunction(SNES snes, Vec x, Vec f, void *ctx)
348: {
349: ApplicationCtx *user = (ApplicationCtx *)ctx;
350: DM da = user->da;
351: PetscScalar *ff, d;
352: const PetscScalar *xx, *FF;
353: PetscInt i, M, xs, xm;
354: Vec xlocal;
356: PetscFunctionBeginUser;
357: PetscCall(DMGetLocalVector(da, &xlocal));
358: /*
359: Scatter ghost points to local vector, using the 2-step process
360: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
361: By placing code between these two statements, computations can
362: be done while messages are in transition.
363: */
364: PetscCall(DMGlobalToLocalBegin(da, x, INSERT_VALUES, xlocal));
365: PetscCall(DMGlobalToLocalEnd(da, x, INSERT_VALUES, xlocal));
367: /*
368: Get pointers to vector data.
369: - The vector xlocal includes ghost point; the vectors x and f do
370: NOT include ghost points.
371: - Using DMDAVecGetArray() allows accessing the values using global ordering
372: */
373: PetscCall(DMDAVecGetArrayRead(da, xlocal, (void *)&xx));
374: PetscCall(DMDAVecGetArray(da, f, &ff));
375: PetscCall(DMDAVecGetArrayRead(da, user->F, (void *)&FF));
377: /*
378: Get local grid boundaries (for 1-dimensional DMDA):
379: xs, xm - starting grid index, width of local grid (no ghost points)
380: */
381: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
382: PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
384: /*
385: Set function values for boundary points; define local interior grid point range:
386: xsi - starting interior grid index
387: xei - ending interior grid index
388: */
389: if (xs == 0) { /* left boundary */
390: ff[0] = xx[0];
391: xs++;
392: xm--;
393: }
394: if (xs + xm == M) { /* right boundary */
395: ff[xs + xm - 1] = xx[xs + xm - 1] - 1.0;
396: xm--;
397: }
399: /*
400: Compute function over locally owned part of the grid (interior points only)
401: */
402: d = 1.0 / (user->h * user->h);
403: for (i = xs; i < xs + xm; i++) ff[i] = d * (xx[i - 1] - 2.0 * xx[i] + xx[i + 1]) + xx[i] * xx[i] - FF[i];
405: /*
406: Restore vectors
407: */
408: PetscCall(DMDAVecRestoreArrayRead(da, xlocal, (void *)&xx));
409: PetscCall(DMDAVecRestoreArray(da, f, &ff));
410: PetscCall(DMDAVecRestoreArrayRead(da, user->F, (void *)&FF));
411: PetscCall(DMRestoreLocalVector(da, &xlocal));
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: /* ------------------------------------------------------------------- */
416: /*
417: FormJacobian - Evaluates Jacobian matrix.
419: Input Parameters:
420: . snes - the SNES context
421: . x - input vector
422: . dummy - optional user-defined context (not used here)
424: Output Parameters:
425: . jac - Jacobian matrix
426: . B - optionally different preconditioning matrix
427: . flag - flag indicating matrix structure
428: */
429: PetscErrorCode FormJacobian(SNES snes, Vec x, Mat jac, Mat B, void *ctx)
430: {
431: ApplicationCtx *user = (ApplicationCtx *)ctx;
432: PetscScalar *xx, d, A[3];
433: PetscInt i, j[3], M, xs, xm;
434: DM da = user->da;
436: PetscFunctionBeginUser;
437: if (user->sjerr) {
438: PetscCall(SNESSetJacobianDomainError(snes));
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
441: /*
442: Get pointer to vector data
443: */
444: PetscCall(DMDAVecGetArrayRead(da, x, &xx));
445: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
447: /*
448: Get range of locally owned matrix
449: */
450: PetscCall(DMDAGetInfo(da, NULL, &M, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
452: /*
453: Determine starting and ending local indices for interior grid points.
454: Set Jacobian entries for boundary points.
455: */
457: if (xs == 0) { /* left boundary */
458: i = 0;
459: A[0] = 1.0;
461: PetscCall(MatSetValues(jac, 1, &i, 1, &i, A, INSERT_VALUES));
462: xs++;
463: xm--;
464: }
465: if (xs + xm == M) { /* right boundary */
466: i = M - 1;
467: A[0] = 1.0;
468: PetscCall(MatSetValues(jac, 1, &i, 1, &i, A, INSERT_VALUES));
469: xm--;
470: }
472: /*
473: Interior grid points
474: - Note that in this case we set all elements for a particular
475: row at once.
476: */
477: d = 1.0 / (user->h * user->h);
478: for (i = xs; i < xs + xm; i++) {
479: j[0] = i - 1;
480: j[1] = i;
481: j[2] = i + 1;
482: A[0] = A[2] = d;
483: A[1] = -2.0 * d + 2.0 * xx[i];
484: PetscCall(MatSetValues(jac, 1, &i, 3, j, A, INSERT_VALUES));
485: }
487: /*
488: Assemble matrix, using the 2-step process:
489: MatAssemblyBegin(), MatAssemblyEnd().
490: By placing code between these two statements, computations can be
491: done while messages are in transition.
493: Also, restore vector.
494: */
496: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
497: PetscCall(DMDAVecRestoreArrayRead(da, x, &xx));
498: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: /* ------------------------------------------------------------------- */
503: /*
504: Monitor - Optional user-defined monitoring routine that views the
505: current iterate with an x-window plot. Set by SNESMonitorSet().
507: Input Parameters:
508: snes - the SNES context
509: its - iteration number
510: norm - 2-norm function value (may be estimated)
511: ctx - optional user-defined context for private data for the
512: monitor routine, as set by SNESMonitorSet()
514: Note:
515: See the manpage for PetscViewerDrawOpen() for useful runtime options,
516: such as -nox to deactivate all x-window output.
517: */
518: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
519: {
520: MonitorCtx *monP = (MonitorCtx *)ctx;
521: Vec x;
523: PetscFunctionBeginUser;
524: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ",SNES Function norm %g\n", its, (double)fnorm));
525: PetscCall(SNESGetSolution(snes, &x));
526: PetscCall(VecView(x, monP->viewer));
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /* ------------------------------------------------------------------- */
531: /*
532: PreCheck - Optional user-defined routine that checks the validity of
533: candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
535: Input Parameters:
536: snes - the SNES context
537: xcurrent - current solution
538: y - search direction and length
540: Output Parameters:
541: y - proposed step (search direction and length) (possibly changed)
542: changed_y - tells if the step has changed or not
543: */
544: PetscErrorCode PreCheck(SNESLineSearch linesearch, Vec xcurrent, Vec y, PetscBool *changed_y, void *ctx)
545: {
546: PetscFunctionBeginUser;
547: *changed_y = PETSC_FALSE;
548: PetscFunctionReturn(PETSC_SUCCESS);
549: }
551: /* ------------------------------------------------------------------- */
552: /*
553: PostCheck - Optional user-defined routine that checks the validity of
554: candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
556: Input Parameters:
557: snes - the SNES context
558: ctx - optional user-defined context for private data for the
559: monitor routine, as set by SNESLineSearchSetPostCheck()
560: xcurrent - current solution
561: y - search direction and length
562: x - the new candidate iterate
564: Output Parameters:
565: y - proposed step (search direction and length) (possibly changed)
566: x - current iterate (possibly modified)
568: */
569: PetscErrorCode PostCheck(SNESLineSearch linesearch, Vec xcurrent, Vec y, Vec x, PetscBool *changed_y, PetscBool *changed_x, void *ctx)
570: {
571: PetscInt i, iter, xs, xm;
572: StepCheckCtx *check;
573: ApplicationCtx *user;
574: PetscScalar *xa, *xa_last, tmp;
575: PetscReal rdiff;
576: DM da;
577: SNES snes;
579: PetscFunctionBeginUser;
580: *changed_x = PETSC_FALSE;
581: *changed_y = PETSC_FALSE;
583: PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
584: check = (StepCheckCtx *)ctx;
585: user = check->user;
586: PetscCall(SNESGetIterationNumber(snes, &iter));
588: /* iteration 1 indicates we are working on the second iteration */
589: if (iter > 0) {
590: da = user->da;
591: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Checking candidate step at iteration %" PetscInt_FMT " with tolerance %g\n", iter, (double)check->tolerance));
593: /* Access local array data */
594: PetscCall(DMDAVecGetArray(da, check->last_step, &xa_last));
595: PetscCall(DMDAVecGetArray(da, x, &xa));
596: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
598: /*
599: If we fail the user-defined check for validity of the candidate iterate,
600: then modify the iterate as we like. (Note that the particular modification
601: below is intended simply to demonstrate how to manipulate this data, not
602: as a meaningful or appropriate choice.)
603: */
604: for (i = xs; i < xs + xm; i++) {
605: if (!PetscAbsScalar(xa[i])) rdiff = 2 * check->tolerance;
606: else rdiff = PetscAbsScalar((xa[i] - xa_last[i]) / xa[i]);
607: if (rdiff > check->tolerance) {
608: tmp = xa[i];
609: xa[i] = .5 * (xa[i] + xa_last[i]);
610: *changed_x = PETSC_TRUE;
611: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Altering entry %" PetscInt_FMT ": x=%g, x_last=%g, diff=%g, x_new=%g\n", i, (double)PetscAbsScalar(tmp), (double)PetscAbsScalar(xa_last[i]), (double)rdiff, (double)PetscAbsScalar(xa[i])));
612: }
613: }
614: PetscCall(DMDAVecRestoreArray(da, check->last_step, &xa_last));
615: PetscCall(DMDAVecRestoreArray(da, x, &xa));
616: }
617: PetscCall(VecCopy(x, check->last_step));
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: /* ------------------------------------------------------------------- */
622: /*
623: PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
624: e.g,
625: mpiexec -n 8 ./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -post_setsubksp -sub_ksp_rtol 1.e-16
626: Set by SNESLineSearchSetPostCheck().
628: Input Parameters:
629: linesearch - the LineSearch context
630: xcurrent - current solution
631: y - search direction and length
632: x - the new candidate iterate
634: Output Parameters:
635: y - proposed step (search direction and length) (possibly changed)
636: x - current iterate (possibly modified)
638: */
639: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch, Vec xcurrent, Vec y, Vec x, PetscBool *changed_y, PetscBool *changed_x, void *ctx)
640: {
641: SetSubKSPCtx *check;
642: PetscInt iter, its, sub_its, maxit;
643: KSP ksp, sub_ksp, *sub_ksps;
644: PC pc;
645: PetscReal ksp_ratio;
646: SNES snes;
648: PetscFunctionBeginUser;
649: PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
650: check = (SetSubKSPCtx *)ctx;
651: PetscCall(SNESGetIterationNumber(snes, &iter));
652: PetscCall(SNESGetKSP(snes, &ksp));
653: PetscCall(KSPGetPC(ksp, &pc));
654: PetscCall(PCBJacobiGetSubKSP(pc, NULL, NULL, &sub_ksps));
655: sub_ksp = sub_ksps[0];
656: PetscCall(KSPGetIterationNumber(ksp, &its)); /* outer KSP iteration number */
657: PetscCall(KSPGetIterationNumber(sub_ksp, &sub_its)); /* inner KSP iteration number */
659: if (iter) {
660: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " ...PostCheck snes iteration %" PetscInt_FMT ", ksp_it %" PetscInt_FMT " %" PetscInt_FMT ", subksp_it %" PetscInt_FMT "\n", iter, check->its0, its, sub_its));
661: ksp_ratio = (PetscReal)its / check->its0;
662: maxit = (PetscInt)(ksp_ratio * sub_its + 0.5);
663: if (maxit < 2) maxit = 2;
664: PetscCall(KSPSetTolerances(sub_ksp, PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, maxit));
665: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " ...ksp_ratio %g, new maxit %" PetscInt_FMT "\n\n", (double)ksp_ratio, maxit));
666: }
667: check->its0 = its; /* save current outer KSP iteration number */
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
671: /* ------------------------------------------------------------------- */
672: /*
673: MatrixFreePreconditioner - This routine demonstrates the use of a
674: user-provided preconditioner. This code implements just the null
675: preconditioner, which of course is not recommended for general use.
677: Input Parameters:
678: + pc - preconditioner
679: - x - input vector
681: Output Parameter:
682: . y - preconditioned vector
683: */
684: PetscErrorCode MatrixFreePreconditioner(PC pc, Vec x, Vec y)
685: {
686: PetscFunctionBeginUser;
687: PetscCall(VecCopy(x, y));
688: PetscFunctionReturn(PETSC_SUCCESS);
689: }
691: /*TEST
693: test:
694: args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
696: test:
697: suffix: 2
698: nsize: 3
699: args: -nox -pc_type asm -mat_type mpiaij -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
701: test:
702: suffix: 3
703: nsize: 2
704: args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
706: test:
707: suffix: 4
708: args: -nox -pre_check_iterates -post_check_iterates
710: test:
711: suffix: 5
712: requires: double !complex !single
713: nsize: 2
714: args: -nox -snes_test_jacobian -snes_test_jacobian_view
716: test:
717: suffix: 6
718: requires: double !complex !single
719: nsize: 4
720: args: -test_jacobian_domain_error -snes_converged_reason -snes_check_jacobian_domain_error 1
722: test:
723: suffix: 7
724: requires: double !complex !single
725: nsize: 4
726: args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontr -snes_check_jacobian_domain_error 1
728: test:
729: suffix: 8
730: requires: double !complex !single
731: nsize: 4
732: args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonrsls -snes_check_jacobian_domain_error 1
734: test:
735: suffix: 9
736: requires: double !complex !single
737: nsize: 4
738: args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonssls -snes_check_jacobian_domain_error 1
740: test:
741: suffix: 10
742: requires: double !complex !single
743: nsize: 4
744: args: -test_jacobian_domain_error -snes_converged_reason -snes_type qn -snes_qn_scale_type jacobian -snes_check_jacobian_domain_error 1
746: test:
747: suffix: 11
748: requires: double !complex !single
749: nsize: 4
750: args: -test_jacobian_domain_error -snes_converged_reason -snes_type ms -snes_ms_type m62 -snes_ms_damping 0.9 -snes_check_jacobian_domain_error 1
752: test:
753: suffix: 12
754: args: -view_initial
755: filter: grep -v "type:"
757: test:
758: suffix: 13
759: requires: double !complex !single
760: nsize: 4
761: args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontrdc -snes_check_jacobian_domain_error 1
763: TEST*/