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;
91: PetscInitialize(&argc, &argv, (char *)0, help);
92: MPI_Comm_rank(PETSC_COMM_WORLD, &ctx.rank);
93: MPI_Comm_size(PETSC_COMM_WORLD, &ctx.size);
94: PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL);
95: ctx.h = 1.0 / (N - 1);
96: ctx.sjerr = PETSC_FALSE;
97: PetscOptionsGetBool(NULL, NULL, "-test_jacobian_domain_error", &ctx.sjerr, NULL);
98: PetscOptionsGetBool(NULL, NULL, "-view_initial", &viewinitial, NULL);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Create nonlinear solver context
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: 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: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, N, 1, 1, NULL, &ctx.da);
114: DMSetFromOptions(ctx.da);
115: 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: DMCreateGlobalVector(ctx.da, &x);
122: VecDuplicate(x, &r);
123: VecDuplicate(x, &F);
124: ctx.F = F;
125: 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: SNESSetFunction(snes, r, FormFunction, &ctx);
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Create matrix data structure; set Jacobian evaluation routine
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: MatCreate(PETSC_COMM_WORLD, &J);
142: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, N, N);
143: MatSetFromOptions(J);
144: MatSeqAIJSetPreallocation(J, 3, NULL);
145: 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: SNESSetJacobian(snes, J, J, FormJacobian, &ctx);
157: /*
158: Optionally allow user-provided preconditioner
159: */
160: PetscOptionsHasName(NULL, NULL, "-user_precond", &flg);
161: if (flg) {
162: KSP ksp;
163: PC pc;
164: SNESGetKSP(snes, &ksp);
165: KSPGetPC(ksp, &pc);
166: PCSetType(pc, PCSHELL);
167: PCShellSetApply(pc, MatrixFreePreconditioner);
168: }
170: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171: Customize nonlinear solver; set runtime options
172: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174: /*
175: Set an optional user-defined monitoring routine
176: */
177: PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 400, 400, &monP.viewer);
178: SNESMonitorSet(snes, Monitor, &monP, 0);
180: /*
181: Set names for some vectors to facilitate monitoring (optional)
182: */
183: PetscObjectSetName((PetscObject)x, "Approximate Solution");
184: 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: 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: SNESGetLineSearch(snes, &linesearch);
197: PetscOptionsHasName(NULL, NULL, "-post_check_iterates", &post_check);
199: if (post_check) {
200: PetscPrintf(PETSC_COMM_WORLD, "Activating post step checking routine\n");
201: SNESLineSearchSetPostCheck(linesearch, PostCheck, &checkP);
202: VecDuplicate(x, &(checkP.last_step));
204: checkP.tolerance = 1.0;
205: checkP.user = &ctx;
207: PetscOptionsGetReal(NULL, NULL, "-check_tol", &checkP.tolerance, NULL);
208: }
210: PetscOptionsHasName(NULL, NULL, "-post_setsubksp", &post_setsubksp);
211: if (post_setsubksp) {
212: PetscPrintf(PETSC_COMM_WORLD, "Activating post setsubksp\n");
213: SNESLineSearchSetPostCheck(linesearch, PostSetSubKSP, &checkP1);
214: }
216: PetscOptionsHasName(NULL, NULL, "-pre_check_iterates", &pre_check);
217: if (pre_check) {
218: PetscPrintf(PETSC_COMM_WORLD, "Activating pre step checking routine\n");
219: 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: SNESGetTolerances(snes, &abstol, &rtol, &stol, &maxit, &maxf);
228: 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: DMDAGetCorners(ctx.da, &xs, NULL, NULL, &xm, NULL, NULL);
241: /*
242: Get pointers to vector data
243: */
244: DMDAVecGetArray(ctx.da, F, &FF);
245: 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: DMDAVecRestoreArray(ctx.da, F, &FF);
261: DMDAVecRestoreArray(ctx.da, U, &UU);
262: if (viewinitial) {
263: VecView(U, 0);
264: 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: FormInitialGuess(x);
278: SNESSolve(snes, NULL, x);
279: SNESGetIterationNumber(snes, &its);
280: 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: VecAXPY(x, none, U);
290: VecNorm(x, NORM_2, &norm);
291: PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g Iterations %" PetscInt_FMT "\n", (double)norm, its);
292: if (ctx.sjerr) {
293: SNESType snestype;
294: SNESGetType(snes, &snestype);
295: 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: PetscViewerDestroy(&monP.viewer);
303: if (post_check) VecDestroy(&checkP.last_step);
304: VecDestroy(&x);
305: VecDestroy(&r);
306: VecDestroy(&U);
307: VecDestroy(&F);
308: MatDestroy(&J);
309: SNESDestroy(&snes);
310: DMDestroy(&ctx.da);
311: 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;
327: VecSet(x, pfive);
328: return 0;
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;
357: 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: DMGlobalToLocalBegin(da, x, INSERT_VALUES, xlocal);
365: 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: DMDAVecGetArrayRead(da, xlocal, (void *)&xx);
374: DMDAVecGetArray(da, f, &ff);
375: 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: DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);
382: 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: DMDAVecRestoreArrayRead(da, xlocal, (void *)&xx);
409: DMDAVecRestoreArray(da, f, &ff);
410: DMDAVecRestoreArrayRead(da, user->F, (void *)&FF);
411: DMRestoreLocalVector(da, &xlocal);
412: return 0;
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;
437: if (user->sjerr) {
438: SNESSetJacobianDomainError(snes);
439: return 0;
440: }
441: /*
442: Get pointer to vector data
443: */
444: DMDAVecGetArrayRead(da, x, &xx);
445: DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL);
447: /*
448: Get range of locally owned matrix
449: */
450: 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: 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: 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: 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: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
497: DMDAVecRestoreArrayRead(da, x, &xx);
498: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
500: return 0;
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;
525: PetscPrintf(PETSC_COMM_WORLD, "iter = %" PetscInt_FMT ",SNES Function norm %g\n", its, (double)fnorm);
526: SNESGetSolution(snes, &x);
527: VecView(x, monP->viewer);
528: return 0;
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: {
548: *changed_y = PETSC_FALSE;
549: return 0;
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;
581: *changed_x = PETSC_FALSE;
582: *changed_y = PETSC_FALSE;
584: SNESLineSearchGetSNES(linesearch, &snes);
585: check = (StepCheckCtx *)ctx;
586: user = check->user;
587: SNESGetIterationNumber(snes, &iter);
589: /* iteration 1 indicates we are working on the second iteration */
590: if (iter > 0) {
591: da = user->da;
592: 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: DMDAVecGetArray(da, check->last_step, &xa_last);
596: DMDAVecGetArray(da, x, &xa);
597: 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: 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: DMDAVecRestoreArray(da, check->last_step, &xa_last);
616: DMDAVecRestoreArray(da, x, &xa);
617: }
618: VecCopy(x, check->last_step);
619: return 0;
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;
650: SNESLineSearchGetSNES(linesearch, &snes);
651: check = (SetSubKSPCtx *)ctx;
652: SNESGetIterationNumber(snes, &iter);
653: SNESGetKSP(snes, &ksp);
654: KSPGetPC(ksp, &pc);
655: PCBJacobiGetSubKSP(pc, NULL, NULL, &sub_ksps);
656: sub_ksp = sub_ksps[0];
657: KSPGetIterationNumber(ksp, &its); /* outer KSP iteration number */
658: KSPGetIterationNumber(sub_ksp, &sub_its); /* inner KSP iteration number */
660: if (iter) {
661: 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: KSPSetTolerances(sub_ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, maxit);
666: 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: return 0;
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: VecCopy(x, y);
688: return 0;
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*/