Actual source code: ex3.c
2: static char help[] = "Newton methods to solve u'' + u^{2} = f in parallel.\n\
3: This example employs a user-defined monitoring routine and optionally a user-defined\n\
4: routine to check candidate iterates produced by line search routines. This code also\n\
5: demonstrates use of the macro __FUNCT__ to define routine names for use in error handling.\n\
6: The command line options include:\n\
7: -check_iterates : activate checking of iterates\n\
8: -check_tol <tol>: set tolerance for iterate checking\n\n";
10: /*T
11: Concepts: SNES^basic parallel example
12: Concepts: SNES^setting a user-defined monitoring routine
13: Concepts: error handling^using the macro __FUNCT__ to define routine names;
14: Processors: n
15: T*/
17: /*
18: Include "petscdraw.h" so that we can use distributed arrays (DMDAs).
19: Include "petscdraw.h" so that we can use PETSc drawing routines.
20: Include "petscsnes.h" so that we can use SNES solvers. Note that this
21: file automatically includes:
22: petscsys.h - base PETSc routines petscvec.h - vectors
23: petscmat.h - matrices
24: petscis.h - index sets petscksp.h - Krylov subspace methods
25: petscviewer.h - viewers petscpc.h - preconditioners
26: petscksp.h - linear solvers
27: */
29: #include <petscdmda.h>
30: #include <petscsnes.h>
32: /*
33: User-defined routines. Note that immediately before each routine below,
34: we define the macro __FUNCT__ to be a string containing the routine name.
35: If defined, this macro is used in the PETSc error handlers to provide a
36: complete traceback of routine names. All PETSc library routines use this
37: macro, and users can optionally employ it as well in their application
38: codes. Note that users can get a traceback of PETSc errors regardless of
39: whether they define __FUNCT__ in application codes; this macro merely
40: provides the added traceback detail of the application routine names.
41: */
42: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
43: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
44: PetscErrorCode FormInitialGuess(Vec);
45: PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
46: PetscErrorCode PreCheck(SNES,Vec,Vec,void*,PetscBool*);
47: PetscErrorCode PostCheck(SNES,Vec,Vec,Vec,void*,PetscBool*,PetscBool*);
48: PetscErrorCode PostSetSubKSP(SNES,Vec,Vec,Vec,void*,PetscBool*,PetscBool*);
50: /*
51: User-defined application context
52: */
53: typedef struct {
54: DM da; /* distributed array */
55: Vec F; /* right-hand-side of PDE */
56: PetscMPIInt rank; /* rank of processor */
57: PetscMPIInt size; /* size of communicator */
58: PetscReal h; /* mesh spacing */
59: } ApplicationCtx;
61: /*
62: User-defined context for monitoring
63: */
64: typedef struct {
65: PetscViewer viewer;
66: } MonitorCtx;
68: /*
69: User-defined context for checking candidate iterates that are
70: determined by line search methods
71: */
72: typedef struct {
73: Vec last_step; /* previous iterate */
74: PetscReal tolerance; /* tolerance for changes between successive iterates */
75: ApplicationCtx *user;
76: } StepCheckCtx;
78: typedef struct {
79: PetscInt its0; /* num of prevous outer KSP iterations */
80: } SetSubKSPCtx;
84: int main(int argc,char **argv)
85: {
86: SNES snes; /* SNES context */
87: Mat J; /* Jacobian matrix */
88: ApplicationCtx ctx; /* user-defined context */
89: Vec x,r,U,F; /* vectors */
90: MonitorCtx monP; /* monitoring context */
91: StepCheckCtx checkP; /* step-checking context */
92: SetSubKSPCtx checkP1;
93: PetscBool pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking
94: candidate iterates */
95: PetscScalar xp,*FF,*UU,none = -1.0;
97: PetscInt its,N = 5,i,maxit,maxf,xs,xm;
98: PetscReal abstol,rtol,stol,norm;
100: PetscInitialize(&argc,&argv,(char *)0,help);
101: MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
102: MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
103: PetscOptionsGetInt(PETSC_NULL,"-n",&N,PETSC_NULL);
104: ctx.h = 1.0/(N-1);
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Create nonlinear solver context
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: SNESCreate(PETSC_COMM_WORLD,&snes);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Create vector data structures; set function evaluation routine
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: /*
117: Create distributed array (DMDA) to manage parallel grid and vectors
118: */
119: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,N,1,1,PETSC_NULL,&ctx.da);
121: /*
122: Extract global and local vectors from DMDA; then duplicate for remaining
123: vectors that are the same types
124: */
125: DMCreateGlobalVector(ctx.da,&x);
126: VecDuplicate(x,&r);
127: VecDuplicate(x,&F); ctx.F = F;
128: VecDuplicate(x,&U);
130: /*
131: Set function evaluation routine and vector. Whenever the nonlinear
132: solver needs to compute the nonlinear function, it will call this
133: routine.
134: - Note that the final routine argument is the user-defined
135: context that provides application-specific data for the
136: function evaluation routine.
137: */
138: SNESSetFunction(snes,r,FormFunction,&ctx);
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Create matrix data structure; set Jacobian evaluation routine
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: MatCreate(PETSC_COMM_WORLD,&J);
145: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
146: MatSetFromOptions(J);
147: MatSeqAIJSetPreallocation(J,3,PETSC_NULL);
148: MatMPIAIJSetPreallocation(J,3,PETSC_NULL,3,PETSC_NULL);
150: /*
151: Set Jacobian matrix data structure and default Jacobian evaluation
152: routine. Whenever the nonlinear solver needs to compute the
153: Jacobian matrix, it will call this routine.
154: - Note that the final routine argument is the user-defined
155: context that provides application-specific data for the
156: Jacobian evaluation routine.
157: */
158: SNESSetJacobian(snes,J,J,FormJacobian,&ctx);
160: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161: Customize nonlinear solver; set runtime options
162: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164: /*
165: Set an optional user-defined monitoring routine
166: */
167: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
168: SNESMonitorSet(snes,Monitor,&monP,0);
170: /*
171: Set names for some vectors to facilitate monitoring (optional)
172: */
173: PetscObjectSetName((PetscObject)x,"Approximate Solution");
174: PetscObjectSetName((PetscObject)U,"Exact Solution");
176: /*
177: Set SNES/KSP/KSP/PC runtime options, e.g.,
178: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
179: */
180: SNESSetFromOptions(snes);
182: /*
183: Set an optional user-defined routine to check the validity of candidate
184: iterates that are determined by line search methods
185: */
186: PetscOptionsHasName(PETSC_NULL,"-post_check_iterates",&post_check);
187: if (post_check) {
188: PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
189: SNESLineSearchSetPostCheck(snes,PostCheck,&checkP);
190: VecDuplicate(x,&(checkP.last_step));
191: checkP.tolerance = 1.0;
192: checkP.user = &ctx;
193: PetscOptionsGetReal(PETSC_NULL,"-check_tol",&checkP.tolerance,PETSC_NULL);
194: }
196: PetscOptionsHasName(PETSC_NULL,"-post_setsubksp",&post_setsubksp);
197: if (post_setsubksp) {
198: PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");
199: SNESLineSearchSetPostCheck(snes,PostSetSubKSP,&checkP1);
200: }
202: PetscOptionsHasName(PETSC_NULL,"-pre_check_iterates",&pre_check);
203: if (pre_check) {
204: PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
205: SNESLineSearchSetPreCheck(snes,PreCheck,&checkP);
206: }
208: /*
209: Print parameters used for convergence testing (optional) ... just
210: to demonstrate this routine; this information is also printed with
211: the option -snes_view
212: */
213: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
214: PetscPrintf(PETSC_COMM_WORLD,"atol=%G, rtol=%G, stol=%G, maxit=%D, maxf=%D\n",abstol,rtol,stol,maxit,maxf);
216: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: Initialize application:
218: Store right-hand-side of PDE and exact solution
219: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
221: /*
222: Get local grid boundaries (for 1-dimensional DMDA):
223: xs, xm - starting grid index, width of local grid (no ghost points)
224: */
225: DMDAGetCorners(ctx.da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
227: /*
228: Get pointers to vector data
229: */
230: DMDAVecGetArray(ctx.da,F,&FF);
231: DMDAVecGetArray(ctx.da,U,&UU);
233: /*
234: Compute local vector entries
235: */
236: xp = ctx.h*xs;
237: for (i=xs; i<xs+xm; i++) {
238: FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
239: UU[i] = xp*xp*xp;
240: xp += ctx.h;
241: }
243: /*
244: Restore vectors
245: */
246: DMDAVecRestoreArray(ctx.da,F,&FF);
247: DMDAVecRestoreArray(ctx.da,U,&UU);
249: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
250: Evaluate initial guess; then solve nonlinear system
251: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
253: /*
254: Note: The user should initialize the vector, x, with the initial guess
255: for the nonlinear solver prior to calling SNESSolve(). In particular,
256: to employ an initial guess of zero, the user should explicitly set
257: this vector to zero by calling VecSet().
258: */
259: FormInitialGuess(x);
260: SNESSolve(snes,PETSC_NULL,x);
261: SNESGetIterationNumber(snes,&its);
262: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n\n",its);
264: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
265: Check solution and clean up
266: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
268: /*
269: Check the error
270: */
271: VecAXPY(x,none,U);
272: VecNorm(x,NORM_2,&norm);
273: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G Iterations %D\n",norm,its);
275: /*
276: Free work space. All PETSc objects should be destroyed when they
277: are no longer needed.
278: */
279: PetscViewerDestroy(&monP.viewer);
280: if (post_check) {VecDestroy(&checkP.last_step);}
281: VecDestroy(&x);
282: VecDestroy(&r);
283: VecDestroy(&U);
284: VecDestroy(&F);
285: MatDestroy(&J);
286: SNESDestroy(&snes);
287: DMDestroy(&ctx.da);
288: PetscFinalize();
290: return(0);
291: }
292: /* ------------------------------------------------------------------- */
295: /*
296: FormInitialGuess - Computes initial guess.
298: Input/Output Parameter:
299: . x - the solution vector
300: */
301: PetscErrorCode FormInitialGuess(Vec x)
302: {
304: PetscScalar pfive = .50;
307: VecSet(x,pfive);
308: return(0);
309: }
310: /* ------------------------------------------------------------------- */
313: /*
314: FormFunction - Evaluates nonlinear function, F(x).
316: Input Parameters:
317: . snes - the SNES context
318: . x - input vector
319: . ctx - optional user-defined context, as set by SNESSetFunction()
321: Output Parameter:
322: . f - function vector
324: Note:
325: The user-defined context can contain any application-specific
326: data needed for the function evaluation.
327: */
328: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
329: {
330: ApplicationCtx *user = (ApplicationCtx*) ctx;
331: DM da = user->da;
332: PetscScalar *xx,*ff,*FF,d;
334: PetscInt i,M,xs,xm;
335: Vec xlocal;
338: DMGetLocalVector(da,&xlocal);
339: /*
340: Scatter ghost points to local vector, using the 2-step process
341: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
342: By placing code between these two statements, computations can
343: be done while messages are in transition.
344: */
345: DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
346: DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
348: /*
349: Get pointers to vector data.
350: - The vector xlocal includes ghost point; the vectors x and f do
351: NOT include ghost points.
352: - Using DMDAVecGetArray() allows accessing the values using global ordering
353: */
354: DMDAVecGetArray(da,xlocal,&xx);
355: DMDAVecGetArray(da,f,&ff);
356: DMDAVecGetArray(da,user->F,&FF);
358: /*
359: Get local grid boundaries (for 1-dimensional DMDA):
360: xs, xm - starting grid index, width of local grid (no ghost points)
361: */
362: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
363: DMDAGetInfo(da,PETSC_NULL,&M,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
364: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
366: /*
367: Set function values for boundary points; define local interior grid point range:
368: xsi - starting interior grid index
369: xei - ending interior grid index
370: */
371: if (xs == 0) { /* left boundary */
372: ff[0] = xx[0];
373: xs++;xm--;
374: }
375: if (xs+xm == M) { /* right boundary */
376: ff[xs+xm-1] = xx[xs+xm-1] - 1.0;
377: xm--;
378: }
380: /*
381: Compute function over locally owned part of the grid (interior points only)
382: */
383: d = 1.0/(user->h*user->h);
384: for (i=xs; i<xs+xm; i++) {
385: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
386: }
388: /*
389: Restore vectors
390: */
391: DMDAVecRestoreArray(da,xlocal,&xx);
392: DMDAVecRestoreArray(da,f,&ff);
393: DMDAVecRestoreArray(da,user->F,&FF);
394: DMRestoreLocalVector(da,&xlocal);
395: return(0);
396: }
397: /* ------------------------------------------------------------------- */
400: /*
401: FormJacobian - Evaluates Jacobian matrix.
403: Input Parameters:
404: . snes - the SNES context
405: . x - input vector
406: . dummy - optional user-defined context (not used here)
408: Output Parameters:
409: . jac - Jacobian matrix
410: . B - optionally different preconditioning matrix
411: . flag - flag indicating matrix structure
412: */
413: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *ctx)
414: {
415: ApplicationCtx *user = (ApplicationCtx*) ctx;
416: PetscScalar *xx,d,A[3];
418: PetscInt i,j[3],M,xs,xm;
419: DM da = user->da;
422: /*
423: Get pointer to vector data
424: */
425: DMDAVecGetArray(da,x,&xx);
426: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
428: /*
429: Get range of locally owned matrix
430: */
431: DMDAGetInfo(da,PETSC_NULL,&M,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
432: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
434: /*
435: Determine starting and ending local indices for interior grid points.
436: Set Jacobian entries for boundary points.
437: */
439: if (xs == 0) { /* left boundary */
440: i = 0; A[0] = 1.0;
441: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
442: xs++;xm--;
443: }
444: if (xs+xm == M) { /* right boundary */
445: i = M-1; A[0] = 1.0;
446: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
447: xm--;
448: }
450: /*
451: Interior grid points
452: - Note that in this case we set all elements for a particular
453: row at once.
454: */
455: d = 1.0/(user->h*user->h);
456: for (i=xs; i<xs+xm; i++) {
457: j[0] = i - 1; j[1] = i; j[2] = i + 1;
458: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
459: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
460: }
462: /*
463: Assemble matrix, using the 2-step process:
464: MatAssemblyBegin(), MatAssemblyEnd().
465: By placing code between these two statements, computations can be
466: done while messages are in transition.
468: Also, restore vector.
469: */
471: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
472: DMDAVecRestoreArray(da,x,&xx);
473: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
475: *flag = SAME_NONZERO_PATTERN;
476: return(0);
477: }
478: /* ------------------------------------------------------------------- */
481: /*
482: Monitor - Optional user-defined monitoring routine that views the
483: current iterate with an x-window plot. Set by SNESMonitorSet().
485: Input Parameters:
486: snes - the SNES context
487: its - iteration number
488: norm - 2-norm function value (may be estimated)
489: ctx - optional user-defined context for private data for the
490: monitor routine, as set by SNESMonitorSet()
492: Note:
493: See the manpage for PetscViewerDrawOpen() for useful runtime options,
494: such as -nox to deactivate all x-window output.
495: */
496: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
497: {
499: MonitorCtx *monP = (MonitorCtx*) ctx;
500: Vec x;
503: PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %G\n",its,fnorm);
504: SNESGetSolution(snes,&x);
505: VecView(x,monP->viewer);
506: return(0);
507: }
509: /* ------------------------------------------------------------------- */
512: /*
513: PreCheck - Optional user-defined routine that checks the validity of
514: candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
516: Input Parameters:
517: snes - the SNES context
518: ctx - optional user-defined context for private data for the
519: monitor routine, as set by SNESLineSearchSetPostCheck()
520: xcurrent - current solution
521: y - search direction and length
523: Output Parameters:
524: y - proposed step (search direction and length) (possibly changed)
525:
526: */
527: PetscErrorCode PreCheck(SNES snes,Vec xcurrent,Vec y,void *ctx,PetscBool *changed_y)
528: {
530: *changed_y = PETSC_FALSE;
531: return(0);
532: }
534: /* ------------------------------------------------------------------- */
537: /*
538: PostCheck - Optional user-defined routine that checks the validity of
539: candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
541: Input Parameters:
542: snes - the SNES context
543: ctx - optional user-defined context for private data for the
544: monitor routine, as set by SNESLineSearchSetPostCheck()
545: xcurrent - current solution
546: y - search direction and length
547: x - the new candidate iterate
549: Output Parameters:
550: y - proposed step (search direction and length) (possibly changed)
551: x - current iterate (possibly modified)
552:
553: */
554: PetscErrorCode PostCheck(SNES snes,Vec xcurrent,Vec y,Vec x,void *ctx,PetscBool *changed_y,PetscBool *changed_x)
555: {
557: PetscInt i,iter,xs,xm;
558: StepCheckCtx *check = (StepCheckCtx*) ctx;
559: ApplicationCtx *user = check->user;
560: PetscScalar *xa,*xa_last,tmp;
561: PetscReal rdiff;
562: DM da;
565: *changed_x = PETSC_FALSE;
566: *changed_y = PETSC_FALSE;
567: SNESGetIterationNumber(snes,&iter);
569: /* iteration 1 indicates we are working on the second iteration */
570: if (iter > 0) {
571: da = user->da;
572: PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %G\n",iter,check->tolerance);
574: /* Access local array data */
575: DMDAVecGetArray(da,check->last_step,&xa_last);
576: DMDAVecGetArray(da,x,&xa);
577: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
579: /*
580: If we fail the user-defined check for validity of the candidate iterate,
581: then modify the iterate as we like. (Note that the particular modification
582: below is intended simply to demonstrate how to manipulate this data, not
583: as a meaningful or appropriate choice.)
584: */
585: for (i=xs; i<xs+xm; i++) {
586: if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance; else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
587: if (rdiff > check->tolerance) {
588: tmp = xa[i];
589: xa[i] = .5*(xa[i] + xa_last[i]);
590: *changed_x = PETSC_TRUE;
591: PetscPrintf(PETSC_COMM_WORLD," Altering entry %D: x=%G, x_last=%G, diff=%G, x_new=%G\n",
592: i,PetscAbsScalar(tmp),PetscAbsScalar(xa_last[i]),rdiff,PetscAbsScalar(xa[i]));
593: }
594: }
595: DMDAVecRestoreArray(da,check->last_step,&xa_last);
596: DMDAVecRestoreArray(da,x,&xa);
597: }
598: VecCopy(x,check->last_step);
599: return(0);
600: }
603: /* ------------------------------------------------------------------- */
606: /*
607: PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
608: e.g,
609: 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
610: Set by SNESLineSearchSetPostCheck().
612: Input Parameters:
613: snes - the SNES context
614: ctx - optional user-defined context for private data for the
615: monitor routine, as set by SNESLineSearchSetPostCheck()
616: xcurrent - current solution
617: y - search direction and length
618: x - the new candidate iterate
620: Output Parameters:
621: y - proposed step (search direction and length) (possibly changed)
622: x - current iterate (possibly modified)
623:
624: */
625: PetscErrorCode PostSetSubKSP(SNES snes,Vec xcurrent,Vec y,Vec x,void *ctx,PetscBool *changed_y,PetscBool *changed_x)
626: {
628: SetSubKSPCtx *check = (SetSubKSPCtx*)ctx;
629: PetscInt iter,its,sub_its,maxit;
630: KSP ksp,sub_ksp,*sub_ksps;
631: PC pc;
632: PetscReal ksp_ratio;
633:
635: SNESGetIterationNumber(snes,&iter);
636: SNESGetKSP(snes,&ksp);
637: KSPGetPC(ksp,&pc);
638: PCBJacobiGetSubKSP(pc,PETSC_NULL,PETSC_NULL,&sub_ksps);
639: sub_ksp = sub_ksps[0];
640: KSPGetIterationNumber(ksp,&its); /* outer KSP iteration number */
641: KSPGetIterationNumber(sub_ksp,&sub_its); /* inner KSP iteration number */
643: if (iter){
644: PetscPrintf(PETSC_COMM_WORLD," ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);
645: ksp_ratio = ((PetscReal)(its))/check->its0;
646: maxit = (PetscInt)(ksp_ratio*sub_its + 0.5);
647: if (maxit < 2) maxit = 2;
648: KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
649: PetscPrintf(PETSC_COMM_WORLD," ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);
650: }
651: check->its0 = its; /* save current outer KSP iteration number */
652: return(0);
653: }