Actual source code: ex3.c
petsc-3.8.4 2018-03-24
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: /*T
11: Concepts: SNES^basic parallel example
12: Concepts: SNES^setting a user-defined monitoring routine
13: Processors: n
14: T*/
16: /*
17: Include "petscdraw.h" so that we can use distributed arrays (DMDAs).
18: Include "petscdraw.h" so that we can use PETSc drawing routines.
19: Include "petscsnes.h" so that we can use SNES solvers. Note that this
20: file automatically includes:
21: petscsys.h - base PETSc routines
22: petscvec.h - vectors
23: petscmat.h - matrices
24: petscis.h - index sets
25: petscksp.h - Krylov subspace methods
26: petscviewer.h - viewers
27: petscpc.h - preconditioners
28: petscksp.h - linear solvers
29: */
31: #include <petscdm.h>
32: #include <petscdmda.h>
33: #include <petscsnes.h>
35: /*
36: User-defined routines.
37: */
38: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
39: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
40: PetscErrorCode FormInitialGuess(Vec);
41: PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
42: PetscErrorCode PreCheck(SNESLineSearch,Vec,Vec,PetscBool*,void*);
43: PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
44: PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
45: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);
47: /*
48: User-defined application context
49: */
50: typedef struct {
51: DM da; /* distributed array */
52: Vec F; /* right-hand-side of PDE */
53: PetscMPIInt rank; /* rank of processor */
54: PetscMPIInt size; /* size of communicator */
55: PetscReal h; /* mesh spacing */
56: } ApplicationCtx;
58: /*
59: User-defined context for monitoring
60: */
61: typedef struct {
62: PetscViewer viewer;
63: } MonitorCtx;
65: /*
66: User-defined context for checking candidate iterates that are
67: determined by line search methods
68: */
69: typedef struct {
70: Vec last_step; /* previous iterate */
71: PetscReal tolerance; /* tolerance for changes between successive iterates */
72: ApplicationCtx *user;
73: } StepCheckCtx;
75: typedef struct {
76: PetscInt its0; /* num of prevous outer KSP iterations */
77: } SetSubKSPCtx;
79: int main(int argc,char **argv)
80: {
81: SNES snes; /* SNES context */
82: SNESLineSearch linesearch; /* SNESLineSearch context */
83: Mat J; /* Jacobian matrix */
84: ApplicationCtx ctx; /* user-defined context */
85: Vec x,r,U,F; /* vectors */
86: MonitorCtx monP; /* monitoring context */
87: StepCheckCtx checkP; /* step-checking context */
88: SetSubKSPCtx checkP1;
89: PetscBool pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking candidate iterates */
90: PetscScalar xp,*FF,*UU,none = -1.0;
92: PetscInt its,N = 5,i,maxit,maxf,xs,xm;
93: PetscReal abstol,rtol,stol,norm;
94: PetscBool flg;
96: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
97: MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
98: MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
99: PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
100: ctx.h = 1.0/(N-1);
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Create nonlinear solver context
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: SNESCreate(PETSC_COMM_WORLD,&snes);
108: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109: Create vector data structures; set function evaluation routine
110: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112: /*
113: Create distributed array (DMDA) to manage parallel grid and vectors
114: */
115: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,1,1,NULL,&ctx.da);
116: DMSetFromOptions(ctx.da);
117: DMSetUp(ctx.da);
119: /*
120: Extract global and local vectors from DMDA; then duplicate for remaining
121: vectors that are the same types
122: */
123: DMCreateGlobalVector(ctx.da,&x);
124: VecDuplicate(x,&r);
125: VecDuplicate(x,&F); ctx.F = F;
126: VecDuplicate(x,&U);
128: /*
129: Set function evaluation routine and vector. Whenever the nonlinear
130: solver needs to compute the nonlinear function, it will call this
131: routine.
132: - Note that the final routine argument is the user-defined
133: context that provides application-specific data for the
134: function evaluation routine.
135: */
136: SNESSetFunction(snes,r,FormFunction,&ctx);
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Create matrix data structure; set Jacobian evaluation routine
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: MatCreate(PETSC_COMM_WORLD,&J);
143: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
144: MatSetFromOptions(J);
145: MatSeqAIJSetPreallocation(J,3,NULL);
146: MatMPIAIJSetPreallocation(J,3,NULL,3,NULL);
148: /*
149: Set Jacobian matrix data structure and default Jacobian evaluation
150: routine. Whenever the nonlinear solver needs to compute the
151: Jacobian matrix, it will call this routine.
152: - Note that the final routine argument is the user-defined
153: context that provides application-specific data for the
154: Jacobian evaluation routine.
155: */
156: SNESSetJacobian(snes,J,J,FormJacobian,&ctx);
158: /*
159: Optionally allow user-provided preconditioner
160: */
161: PetscOptionsHasName(NULL,NULL,"-user_precond",&flg);
162: if (flg) {
163: KSP ksp;
164: PC pc;
165: SNESGetKSP(snes,&ksp);
166: KSPGetPC(ksp,&pc);
167: PCSetType(pc,PCSHELL);
168: PCShellSetApply(pc,MatrixFreePreconditioner);
169: }
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Customize nonlinear solver; set runtime options
173: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: /*
176: Set an optional user-defined monitoring routine
177: */
178: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
179: SNESMonitorSet(snes,Monitor,&monP,0);
181: /*
182: Set names for some vectors to facilitate monitoring (optional)
183: */
184: PetscObjectSetName((PetscObject)x,"Approximate Solution");
185: PetscObjectSetName((PetscObject)U,"Exact Solution");
187: /*
188: Set SNES/KSP/KSP/PC runtime options, e.g.,
189: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
190: */
191: SNESSetFromOptions(snes);
193: /*
194: Set an optional user-defined routine to check the validity of candidate
195: iterates that are determined by line search methods
196: */
197: SNESGetLineSearch(snes, &linesearch);
198: PetscOptionsHasName(NULL,NULL,"-post_check_iterates",&post_check);
200: if (post_check) {
201: PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
202: SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP);
203: VecDuplicate(x,&(checkP.last_step));
205: checkP.tolerance = 1.0;
206: checkP.user = &ctx;
208: PetscOptionsGetReal(NULL,NULL,"-check_tol",&checkP.tolerance,NULL);
209: }
211: PetscOptionsHasName(NULL,NULL,"-post_setsubksp",&post_setsubksp);
212: if (post_setsubksp) {
213: PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");
214: SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1);
215: }
217: PetscOptionsHasName(NULL,NULL,"-pre_check_iterates",&pre_check);
218: if (pre_check) {
219: PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
220: SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP);
221: }
223: /*
224: Print parameters used for convergence testing (optional) ... just
225: to demonstrate this routine; this information is also printed with
226: the option -snes_view
227: */
228: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
229: PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);
231: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232: Initialize application:
233: Store right-hand-side of PDE and exact solution
234: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
236: /*
237: Get local grid boundaries (for 1-dimensional DMDA):
238: xs, xm - starting grid index, width of local grid (no ghost points)
239: */
240: DMDAGetCorners(ctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
242: /*
243: Get pointers to vector data
244: */
245: DMDAVecGetArray(ctx.da,F,&FF);
246: DMDAVecGetArray(ctx.da,U,&UU);
248: /*
249: Compute local vector entries
250: */
251: xp = ctx.h*xs;
252: for (i=xs; i<xs+xm; i++) {
253: FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
254: UU[i] = xp*xp*xp;
255: xp += ctx.h;
256: }
258: /*
259: Restore vectors
260: */
261: DMDAVecRestoreArray(ctx.da,F,&FF);
262: DMDAVecRestoreArray(ctx.da,U,&UU);
264: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
265: Evaluate initial guess; then solve nonlinear system
266: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
268: /*
269: Note: The user should initialize the vector, x, with the initial guess
270: for the nonlinear solver prior to calling SNESSolve(). In particular,
271: to employ an initial guess of zero, the user should explicitly set
272: this vector to zero by calling VecSet().
273: */
274: FormInitialGuess(x);
275: SNESSolve(snes,NULL,x);
276: SNESGetIterationNumber(snes,&its);
277: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
279: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
280: Check solution and clean up
281: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
283: /*
284: Check the error
285: */
286: VecAXPY(x,none,U);
287: VecNorm(x,NORM_2,&norm);
288: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)norm,its);
290: /*
291: Free work space. All PETSc objects should be destroyed when they
292: are no longer needed.
293: */
294: PetscViewerDestroy(&monP.viewer);
295: if (post_check) {VecDestroy(&checkP.last_step);}
296: VecDestroy(&x);
297: VecDestroy(&r);
298: VecDestroy(&U);
299: VecDestroy(&F);
300: MatDestroy(&J);
301: SNESDestroy(&snes);
302: DMDestroy(&ctx.da);
303: PetscFinalize();
304: return ierr;
305: }
307: /* ------------------------------------------------------------------- */
308: /*
309: FormInitialGuess - Computes initial guess.
311: Input/Output Parameter:
312: . x - the solution vector
313: */
314: PetscErrorCode FormInitialGuess(Vec x)
315: {
317: PetscScalar pfive = .50;
320: VecSet(x,pfive);
321: return(0);
322: }
324: /* ------------------------------------------------------------------- */
325: /*
326: FormFunction - Evaluates nonlinear function, F(x).
328: Input Parameters:
329: . snes - the SNES context
330: . x - input vector
331: . ctx - optional user-defined context, as set by SNESSetFunction()
333: Output Parameter:
334: . f - function vector
336: Note:
337: The user-defined context can contain any application-specific
338: data needed for the function evaluation.
339: */
340: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
341: {
342: ApplicationCtx *user = (ApplicationCtx*) ctx;
343: DM da = user->da;
344: PetscScalar *xx,*ff,*FF,d;
346: PetscInt i,M,xs,xm;
347: Vec xlocal;
350: DMGetLocalVector(da,&xlocal);
351: /*
352: Scatter ghost points to local vector, using the 2-step process
353: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
354: By placing code between these two statements, computations can
355: be done while messages are in transition.
356: */
357: DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
358: DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
360: /*
361: Get pointers to vector data.
362: - The vector xlocal includes ghost point; the vectors x and f do
363: NOT include ghost points.
364: - Using DMDAVecGetArray() allows accessing the values using global ordering
365: */
366: DMDAVecGetArray(da,xlocal,&xx);
367: DMDAVecGetArray(da,f,&ff);
368: DMDAVecGetArray(da,user->F,&FF);
370: /*
371: Get local grid boundaries (for 1-dimensional DMDA):
372: xs, xm - starting grid index, width of local grid (no ghost points)
373: */
374: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
375: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
377: /*
378: Set function values for boundary points; define local interior grid point range:
379: xsi - starting interior grid index
380: xei - ending interior grid index
381: */
382: if (xs == 0) { /* left boundary */
383: ff[0] = xx[0];
384: xs++;xm--;
385: }
386: if (xs+xm == M) { /* right boundary */
387: ff[xs+xm-1] = xx[xs+xm-1] - 1.0;
388: xm--;
389: }
391: /*
392: Compute function over locally owned part of the grid (interior points only)
393: */
394: d = 1.0/(user->h*user->h);
395: 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];
397: /*
398: Restore vectors
399: */
400: DMDAVecRestoreArray(da,xlocal,&xx);
401: DMDAVecRestoreArray(da,f,&ff);
402: DMDAVecRestoreArray(da,user->F,&FF);
403: DMRestoreLocalVector(da,&xlocal);
404: return(0);
405: }
407: /* ------------------------------------------------------------------- */
408: /*
409: FormJacobian - Evaluates Jacobian matrix.
411: Input Parameters:
412: . snes - the SNES context
413: . x - input vector
414: . dummy - optional user-defined context (not used here)
416: Output Parameters:
417: . jac - Jacobian matrix
418: . B - optionally different preconditioning matrix
419: . flag - flag indicating matrix structure
420: */
421: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
422: {
423: ApplicationCtx *user = (ApplicationCtx*) ctx;
424: PetscScalar *xx,d,A[3];
426: PetscInt i,j[3],M,xs,xm;
427: DM da = user->da;
430: /*
431: Get pointer to vector data
432: */
433: DMDAVecGetArrayRead(da,x,&xx);
434: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
436: /*
437: Get range of locally owned matrix
438: */
439: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
441: /*
442: Determine starting and ending local indices for interior grid points.
443: Set Jacobian entries for boundary points.
444: */
446: if (xs == 0) { /* left boundary */
447: i = 0; A[0] = 1.0;
449: MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
450: xs++;xm--;
451: }
452: if (xs+xm == M) { /* right boundary */
453: i = M-1;
454: A[0] = 1.0;
455: MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
456: xm--;
457: }
459: /*
460: Interior grid points
461: - Note that in this case we set all elements for a particular
462: row at once.
463: */
464: d = 1.0/(user->h*user->h);
465: for (i=xs; i<xs+xm; i++) {
466: j[0] = i - 1; j[1] = i; j[2] = i + 1;
467: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
468: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
469: }
471: /*
472: Assemble matrix, using the 2-step process:
473: MatAssemblyBegin(), MatAssemblyEnd().
474: By placing code between these two statements, computations can be
475: done while messages are in transition.
477: Also, restore vector.
478: */
480: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
481: DMDAVecRestoreArrayRead(da,x,&xx);
482: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
484: return(0);
485: }
487: /* ------------------------------------------------------------------- */
488: /*
489: Monitor - Optional user-defined monitoring routine that views the
490: current iterate with an x-window plot. Set by SNESMonitorSet().
492: Input Parameters:
493: snes - the SNES context
494: its - iteration number
495: norm - 2-norm function value (may be estimated)
496: ctx - optional user-defined context for private data for the
497: monitor routine, as set by SNESMonitorSet()
499: Note:
500: See the manpage for PetscViewerDrawOpen() for useful runtime options,
501: such as -nox to deactivate all x-window output.
502: */
503: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
504: {
506: MonitorCtx *monP = (MonitorCtx*) ctx;
507: Vec x;
510: PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %g\n",its,(double)fnorm);
511: SNESGetSolution(snes,&x);
512: VecView(x,monP->viewer);
513: return(0);
514: }
516: /* ------------------------------------------------------------------- */
517: /*
518: PreCheck - Optional user-defined routine that checks the validity of
519: candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
521: Input Parameters:
522: snes - the SNES context
523: xcurrent - current solution
524: y - search direction and length
526: Output Parameters:
527: y - proposed step (search direction and length) (possibly changed)
528: changed_y - tells if the step has changed or not
529: */
530: PetscErrorCode PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y, PetscBool *changed_y, void * ctx)
531: {
533: *changed_y = PETSC_FALSE;
534: return(0);
535: }
537: /* ------------------------------------------------------------------- */
538: /*
539: PostCheck - Optional user-defined routine that checks the validity of
540: candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
542: Input Parameters:
543: snes - the SNES context
544: ctx - optional user-defined context for private data for the
545: monitor routine, as set by SNESLineSearchSetPostCheck()
546: xcurrent - current solution
547: y - search direction and length
548: x - the new candidate iterate
550: Output Parameters:
551: y - proposed step (search direction and length) (possibly changed)
552: x - current iterate (possibly modified)
554: */
555: PetscErrorCode PostCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
556: {
558: PetscInt i,iter,xs,xm;
559: StepCheckCtx *check;
560: ApplicationCtx *user;
561: PetscScalar *xa,*xa_last,tmp;
562: PetscReal rdiff;
563: DM da;
564: SNES snes;
567: *changed_x = PETSC_FALSE;
568: *changed_y = PETSC_FALSE;
570: SNESLineSearchGetSNES(linesearch, &snes);
571: check = (StepCheckCtx*)ctx;
572: user = check->user;
573: SNESGetIterationNumber(snes,&iter);
575: /* iteration 1 indicates we are working on the second iteration */
576: if (iter > 0) {
577: da = user->da;
578: PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %g\n",iter,(double)check->tolerance);
580: /* Access local array data */
581: DMDAVecGetArray(da,check->last_step,&xa_last);
582: DMDAVecGetArray(da,x,&xa);
583: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
585: /*
586: If we fail the user-defined check for validity of the candidate iterate,
587: then modify the iterate as we like. (Note that the particular modification
588: below is intended simply to demonstrate how to manipulate this data, not
589: as a meaningful or appropriate choice.)
590: */
591: for (i=xs; i<xs+xm; i++) {
592: if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance;
593: else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
594: if (rdiff > check->tolerance) {
595: tmp = xa[i];
596: xa[i] = .5*(xa[i] + xa_last[i]);
597: *changed_x = PETSC_TRUE;
598: PetscPrintf(PETSC_COMM_WORLD," Altering entry %D: 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]));
599: }
600: }
601: DMDAVecRestoreArray(da,check->last_step,&xa_last);
602: DMDAVecRestoreArray(da,x,&xa);
603: }
604: VecCopy(x,check->last_step);
605: return(0);
606: }
608: /* ------------------------------------------------------------------- */
609: /*
610: PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
611: e.g,
612: 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
613: Set by SNESLineSearchSetPostCheck().
615: Input Parameters:
616: linesearch - the LineSearch context
617: xcurrent - current solution
618: y - search direction and length
619: x - the new candidate iterate
621: Output Parameters:
622: y - proposed step (search direction and length) (possibly changed)
623: x - current iterate (possibly modified)
625: */
626: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
627: {
629: SetSubKSPCtx *check;
630: PetscInt iter,its,sub_its,maxit;
631: KSP ksp,sub_ksp,*sub_ksps;
632: PC pc;
633: PetscReal ksp_ratio;
634: SNES snes;
637: SNESLineSearchGetSNES(linesearch, &snes);
638: check = (SetSubKSPCtx*)ctx;
639: SNESGetIterationNumber(snes,&iter);
640: SNESGetKSP(snes,&ksp);
641: KSPGetPC(ksp,&pc);
642: PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);
643: sub_ksp = sub_ksps[0];
644: KSPGetIterationNumber(ksp,&its); /* outer KSP iteration number */
645: KSPGetIterationNumber(sub_ksp,&sub_its); /* inner KSP iteration number */
647: if (iter) {
648: PetscPrintf(PETSC_COMM_WORLD," ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);
649: ksp_ratio = ((PetscReal)(its))/check->its0;
650: maxit = (PetscInt)(ksp_ratio*sub_its + 0.5);
651: if (maxit < 2) maxit = 2;
652: KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
653: PetscPrintf(PETSC_COMM_WORLD," ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);
654: }
655: check->its0 = its; /* save current outer KSP iteration number */
656: return(0);
657: }
659: /* ------------------------------------------------------------------- */
660: /*
661: MatrixFreePreconditioner - This routine demonstrates the use of a
662: user-provided preconditioner. This code implements just the null
663: preconditioner, which of course is not recommended for general use.
665: Input Parameters:
666: + pc - preconditioner
667: - x - input vector
669: Output Parameter:
670: . y - preconditioned vector
671: */
672: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
673: {
675: VecCopy(x,y);
676: return 0;
677: }