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