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, (char *)0, 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));
500: PetscFunctionReturn(PETSC_SUCCESS);
501: }
503: /* ------------------------------------------------------------------- */
504: /*
505: Monitor - Optional user-defined monitoring routine that views the
506: current iterate with an x-window plot. Set by SNESMonitorSet().
508: Input Parameters:
509: snes - the SNES context
510: its - iteration number
511: norm - 2-norm function value (may be estimated)
512: ctx - optional user-defined context for private data for the
513: monitor routine, as set by SNESMonitorSet()
515: Note:
516: See the manpage for PetscViewerDrawOpen() for useful runtime options,
517: such as -nox to deactivate all x-window output.
518: */
519: PetscErrorCode Monitor(SNES snes, PetscInt its, PetscReal fnorm, void *ctx)
520: {
521: MonitorCtx *monP = (MonitorCtx *)ctx;
522: Vec x;
524: PetscFunctionBeginUser;
525: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ",SNES Function norm %g\n", its, (double)fnorm));
526: PetscCall(SNESGetSolution(snes, &x));
527: PetscCall(VecView(x, monP->viewer));
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: /* ------------------------------------------------------------------- */
532: /*
533: PreCheck - Optional user-defined routine that checks the validity of
534: candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
536: Input Parameters:
537: snes - the SNES context
538: xcurrent - current solution
539: y - search direction and length
541: Output Parameters:
542: y - proposed step (search direction and length) (possibly changed)
543: changed_y - tells if the step has changed or not
544: */
545: PetscErrorCode PreCheck(SNESLineSearch linesearch, Vec xcurrent, Vec y, PetscBool *changed_y, void *ctx)
546: {
547: PetscFunctionBeginUser;
548: *changed_y = PETSC_FALSE;
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: /* ------------------------------------------------------------------- */
553: /*
554: PostCheck - Optional user-defined routine that checks the validity of
555: candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
557: Input Parameters:
558: snes - the SNES context
559: ctx - optional user-defined context for private data for the
560: monitor routine, as set by SNESLineSearchSetPostCheck()
561: xcurrent - current solution
562: y - search direction and length
563: x - the new candidate iterate
565: Output Parameters:
566: y - proposed step (search direction and length) (possibly changed)
567: x - current iterate (possibly modified)
569: */
570: PetscErrorCode PostCheck(SNESLineSearch linesearch, Vec xcurrent, Vec y, Vec x, PetscBool *changed_y, PetscBool *changed_x, void *ctx)
571: {
572: PetscInt i, iter, xs, xm;
573: StepCheckCtx *check;
574: ApplicationCtx *user;
575: PetscScalar *xa, *xa_last, tmp;
576: PetscReal rdiff;
577: DM da;
578: SNES snes;
580: PetscFunctionBeginUser;
581: *changed_x = PETSC_FALSE;
582: *changed_y = PETSC_FALSE;
584: PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
585: check = (StepCheckCtx *)ctx;
586: user = check->user;
587: PetscCall(SNESGetIterationNumber(snes, &iter));
589: /* iteration 1 indicates we are working on the second iteration */
590: if (iter > 0) {
591: da = user->da;
592: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Checking candidate step at iteration %" PetscInt_FMT " with tolerance %g\n", iter, (double)check->tolerance));
594: /* Access local array data */
595: PetscCall(DMDAVecGetArray(da, check->last_step, &xa_last));
596: PetscCall(DMDAVecGetArray(da, x, &xa));
597: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
599: /*
600: If we fail the user-defined check for validity of the candidate iterate,
601: then modify the iterate as we like. (Note that the particular modification
602: below is intended simply to demonstrate how to manipulate this data, not
603: as a meaningful or appropriate choice.)
604: */
605: for (i = xs; i < xs + xm; i++) {
606: if (!PetscAbsScalar(xa[i])) rdiff = 2 * check->tolerance;
607: else rdiff = PetscAbsScalar((xa[i] - xa_last[i]) / xa[i]);
608: if (rdiff > check->tolerance) {
609: tmp = xa[i];
610: xa[i] = .5 * (xa[i] + xa_last[i]);
611: *changed_x = PETSC_TRUE;
612: 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])));
613: }
614: }
615: PetscCall(DMDAVecRestoreArray(da, check->last_step, &xa_last));
616: PetscCall(DMDAVecRestoreArray(da, x, &xa));
617: }
618: PetscCall(VecCopy(x, check->last_step));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: /* ------------------------------------------------------------------- */
623: /*
624: PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
625: e.g,
626: 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
627: Set by SNESLineSearchSetPostCheck().
629: Input Parameters:
630: linesearch - the LineSearch context
631: xcurrent - current solution
632: y - search direction and length
633: x - the new candidate iterate
635: Output Parameters:
636: y - proposed step (search direction and length) (possibly changed)
637: x - current iterate (possibly modified)
639: */
640: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch, Vec xcurrent, Vec y, Vec x, PetscBool *changed_y, PetscBool *changed_x, void *ctx)
641: {
642: SetSubKSPCtx *check;
643: PetscInt iter, its, sub_its, maxit;
644: KSP ksp, sub_ksp, *sub_ksps;
645: PC pc;
646: PetscReal ksp_ratio;
647: SNES snes;
649: PetscFunctionBeginUser;
650: PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
651: check = (SetSubKSPCtx *)ctx;
652: PetscCall(SNESGetIterationNumber(snes, &iter));
653: PetscCall(SNESGetKSP(snes, &ksp));
654: PetscCall(KSPGetPC(ksp, &pc));
655: PetscCall(PCBJacobiGetSubKSP(pc, NULL, NULL, &sub_ksps));
656: sub_ksp = sub_ksps[0];
657: PetscCall(KSPGetIterationNumber(ksp, &its)); /* outer KSP iteration number */
658: PetscCall(KSPGetIterationNumber(sub_ksp, &sub_its)); /* inner KSP iteration number */
660: if (iter) {
661: 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));
662: ksp_ratio = ((PetscReal)(its)) / check->its0;
663: maxit = (PetscInt)(ksp_ratio * sub_its + 0.5);
664: if (maxit < 2) maxit = 2;
665: PetscCall(KSPSetTolerances(sub_ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, maxit));
666: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " ...ksp_ratio %g, new maxit %" PetscInt_FMT "\n\n", (double)ksp_ratio, maxit));
667: }
668: check->its0 = its; /* save current outer KSP iteration number */
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: /* ------------------------------------------------------------------- */
673: /*
674: MatrixFreePreconditioner - This routine demonstrates the use of a
675: user-provided preconditioner. This code implements just the null
676: preconditioner, which of course is not recommended for general use.
678: Input Parameters:
679: + pc - preconditioner
680: - x - input vector
682: Output Parameter:
683: . y - preconditioned vector
684: */
685: PetscErrorCode MatrixFreePreconditioner(PC pc, Vec x, Vec y)
686: {
687: PetscFunctionBeginUser;
688: PetscCall(VecCopy(x, y));
689: PetscFunctionReturn(PETSC_SUCCESS);
690: }
692: /*TEST
694: test:
695: args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
697: test:
698: suffix: 2
699: nsize: 3
700: args: -nox -pc_type asm -mat_type mpiaij -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
702: test:
703: suffix: 3
704: nsize: 2
705: args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always
707: test:
708: suffix: 4
709: args: -nox -pre_check_iterates -post_check_iterates
711: test:
712: suffix: 5
713: requires: double !complex !single
714: nsize: 2
715: args: -nox -snes_test_jacobian -snes_test_jacobian_view
717: test:
718: suffix: 6
719: requires: double !complex !single
720: nsize: 4
721: args: -test_jacobian_domain_error -snes_converged_reason -snes_check_jacobian_domain_error 1
723: test:
724: suffix: 7
725: requires: double !complex !single
726: nsize: 4
727: args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontr -snes_check_jacobian_domain_error 1
729: test:
730: suffix: 8
731: requires: double !complex !single
732: nsize: 4
733: args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonrsls -snes_check_jacobian_domain_error 1
735: test:
736: suffix: 9
737: requires: double !complex !single
738: nsize: 4
739: args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonssls -snes_check_jacobian_domain_error 1
741: test:
742: suffix: 10
743: requires: double !complex !single
744: nsize: 4
745: args: -test_jacobian_domain_error -snes_converged_reason -snes_type qn -snes_qn_scale_type jacobian -snes_check_jacobian_domain_error 1
747: test:
748: suffix: 11
749: requires: double !complex !single
750: nsize: 4
751: 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
753: test:
754: suffix: 12
755: args: -view_initial
756: filter: grep -v "type:"
758: test:
759: suffix: 13
760: requires: double !complex !single
761: nsize: 4
762: args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontrdc -snes_check_jacobian_domain_error 1
764: TEST*/