Actual source code: ex3.c
petsc-3.5.4 2015-05-23
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 <petscdm.h>
30: #include <petscdmda.h>
31: #include <petscsnes.h>
33: /*
34: User-defined routines. Note that immediately before each routine below,
35: we define the macro __FUNCT__ to be a string containing the routine name.
36: If defined, this macro is used in the PETSc error handlers to provide a
37: complete traceback of routine names. All PETSc library routines use this
38: macro, and users can optionally employ it as well in their application
39: codes. Note that users can get a traceback of PETSc errors regardless of
40: whether they define __FUNCT__ in application codes; this macro merely
41: provides the added traceback detail of the application routine names.
42: */
43: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
44: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
45: PetscErrorCode FormInitialGuess(Vec);
46: PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
47: PetscErrorCode PreCheck(SNESLineSearch,Vec,Vec,PetscBool*,void*);
48: PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
49: PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
50: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);
52: /*
53: User-defined application context
54: */
55: typedef struct {
56: DM da; /* distributed array */
57: Vec F; /* right-hand-side of PDE */
58: PetscMPIInt rank; /* rank of processor */
59: PetscMPIInt size; /* size of communicator */
60: PetscReal h; /* mesh spacing */
61: } ApplicationCtx;
63: /*
64: User-defined context for monitoring
65: */
66: typedef struct {
67: PetscViewer viewer;
68: } MonitorCtx;
70: /*
71: User-defined context for checking candidate iterates that are
72: determined by line search methods
73: */
74: typedef struct {
75: Vec last_step; /* previous iterate */
76: PetscReal tolerance; /* tolerance for changes between successive iterates */
77: ApplicationCtx *user;
78: } StepCheckCtx;
80: typedef struct {
81: PetscInt its0; /* num of prevous outer KSP iterations */
82: } SetSubKSPCtx;
86: int main(int argc,char **argv)
87: {
88: SNES snes; /* SNES context */
89: SNESLineSearch linesearch; /* SNESLineSearch context */
90: Mat J; /* Jacobian matrix */
91: ApplicationCtx ctx; /* user-defined context */
92: Vec x,r,U,F; /* vectors */
93: MonitorCtx monP; /* monitoring context */
94: StepCheckCtx checkP; /* step-checking context */
95: SetSubKSPCtx checkP1;
96: PetscBool pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking candidate iterates */
97: PetscScalar xp,*FF,*UU,none = -1.0;
99: PetscInt its,N = 5,i,maxit,maxf,xs,xm;
100: PetscReal abstol,rtol,stol,norm;
101: PetscBool flg;
104: PetscInitialize(&argc,&argv,(char*)0,help);
105: MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
106: MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
107: PetscOptionsGetInt(NULL,"-n",&N,NULL);
108: ctx.h = 1.0/(N-1);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Create nonlinear solver context
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114: SNESCreate(PETSC_COMM_WORLD,&snes);
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Create vector data structures; set function evaluation routine
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: /*
121: Create distributed array (DMDA) to manage parallel grid and vectors
122: */
123: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,1,1,NULL,&ctx.da);
125: /*
126: Extract global and local vectors from DMDA; then duplicate for remaining
127: vectors that are the same types
128: */
129: DMCreateGlobalVector(ctx.da,&x);
130: VecDuplicate(x,&r);
131: VecDuplicate(x,&F); ctx.F = F;
132: VecDuplicate(x,&U);
134: /*
135: Set function evaluation routine and vector. Whenever the nonlinear
136: solver needs to compute the nonlinear function, it will call this
137: routine.
138: - Note that the final routine argument is the user-defined
139: context that provides application-specific data for the
140: function evaluation routine.
141: */
142: SNESSetFunction(snes,r,FormFunction,&ctx);
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Create matrix data structure; set Jacobian evaluation routine
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: MatCreate(PETSC_COMM_WORLD,&J);
149: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
150: MatSetFromOptions(J);
151: MatSeqAIJSetPreallocation(J,3,NULL);
152: MatMPIAIJSetPreallocation(J,3,NULL,3,NULL);
154: /*
155: Set Jacobian matrix data structure and default Jacobian evaluation
156: routine. Whenever the nonlinear solver needs to compute the
157: Jacobian matrix, it will call this routine.
158: - Note that the final routine argument is the user-defined
159: context that provides application-specific data for the
160: Jacobian evaluation routine.
161: */
162: SNESSetJacobian(snes,J,J,FormJacobian,&ctx);
164: /*
165: Optional allow user provided preconditioner
166: */
167: PetscOptionsHasName(NULL,"-user_precond",&flg);
168: if (flg) {
169: KSP ksp;
170: PC pc;
171: SNESGetKSP(snes,&ksp);
172: KSPGetPC(ksp,&pc);
173: PCSetType(pc,PCSHELL);
174: PCShellSetApply(pc,MatrixFreePreconditioner);
175: }
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Customize nonlinear solver; set runtime options
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181: /*
182: Set an optional user-defined monitoring routine
183: */
184: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
185: SNESMonitorSet(snes,Monitor,&monP,0);
187: /*
188: Set names for some vectors to facilitate monitoring (optional)
189: */
190: PetscObjectSetName((PetscObject)x,"Approximate Solution");
191: PetscObjectSetName((PetscObject)U,"Exact Solution");
193: /*
194: Set SNES/KSP/KSP/PC runtime options, e.g.,
195: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
196: */
197: SNESSetFromOptions(snes);
199: /*
200: Set an optional user-defined routine to check the validity of candidate
201: iterates that are determined by line search methods
202: */
203: SNESGetLineSearch(snes, &linesearch);
204: PetscOptionsHasName(NULL,"-post_check_iterates",&post_check);
206: if (post_check) {
207: PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
208: SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP);
209: VecDuplicate(x,&(checkP.last_step));
211: checkP.tolerance = 1.0;
212: checkP.user = &ctx;
214: PetscOptionsGetReal(NULL,"-check_tol",&checkP.tolerance,NULL);
215: }
217: PetscOptionsHasName(NULL,"-post_setsubksp",&post_setsubksp);
218: if (post_setsubksp) {
219: PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");
220: SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1);
221: }
223: PetscOptionsHasName(NULL,"-pre_check_iterates",&pre_check);
224: if (pre_check) {
225: PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
226: SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP);
227: }
229: /*
230: Print parameters used for convergence testing (optional) ... just
231: to demonstrate this routine; this information is also printed with
232: the option -snes_view
233: */
234: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
235: PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);
237: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238: Initialize application:
239: Store right-hand-side of PDE and exact solution
240: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
242: /*
243: Get local grid boundaries (for 1-dimensional DMDA):
244: xs, xm - starting grid index, width of local grid (no ghost points)
245: */
246: DMDAGetCorners(ctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
248: /*
249: Get pointers to vector data
250: */
251: DMDAVecGetArray(ctx.da,F,&FF);
252: DMDAVecGetArray(ctx.da,U,&UU);
254: /*
255: Compute local vector entries
256: */
257: xp = ctx.h*xs;
258: for (i=xs; i<xs+xm; i++) {
259: FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
260: UU[i] = xp*xp*xp;
261: xp += ctx.h;
262: }
264: /*
265: Restore vectors
266: */
267: DMDAVecRestoreArray(ctx.da,F,&FF);
268: DMDAVecRestoreArray(ctx.da,U,&UU);
270: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
271: Evaluate initial guess; then solve nonlinear system
272: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274: /*
275: Note: The user should initialize the vector, x, with the initial guess
276: for the nonlinear solver prior to calling SNESSolve(). In particular,
277: to employ an initial guess of zero, the user should explicitly set
278: this vector to zero by calling VecSet().
279: */
280: FormInitialGuess(x);
281: SNESSolve(snes,NULL,x);
282: SNESGetIterationNumber(snes,&its);
283: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
285: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286: Check solution and clean up
287: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
289: /*
290: Check the error
291: */
292: VecAXPY(x,none,U);
293: VecNorm(x,NORM_2,&norm);
294: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)norm,its);
296: /*
297: Free work space. All PETSc objects should be destroyed when they
298: are no longer needed.
299: */
300: PetscViewerDestroy(&monP.viewer);
301: if (post_check) {VecDestroy(&checkP.last_step);}
302: VecDestroy(&x);
303: VecDestroy(&r);
304: VecDestroy(&U);
305: VecDestroy(&F);
306: MatDestroy(&J);
307: SNESDestroy(&snes);
308: DMDestroy(&ctx.da);
309: PetscFinalize();
310: return(0);
311: }
312: /* ------------------------------------------------------------------- */
315: /*
316: FormInitialGuess - Computes initial guess.
318: Input/Output Parameter:
319: . x - the solution vector
320: */
321: PetscErrorCode FormInitialGuess(Vec x)
322: {
324: PetscScalar pfive = .50;
327: VecSet(x,pfive);
328: return(0);
329: }
330: /* ------------------------------------------------------------------- */
333: /*
334: FormFunction - Evaluates nonlinear function, F(x).
336: Input Parameters:
337: . snes - the SNES context
338: . x - input vector
339: . ctx - optional user-defined context, as set by SNESSetFunction()
341: Output Parameter:
342: . f - function vector
344: Note:
345: The user-defined context can contain any application-specific
346: data needed for the function evaluation.
347: */
348: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
349: {
350: ApplicationCtx *user = (ApplicationCtx*) ctx;
351: DM da = user->da;
352: PetscScalar *xx,*ff,*FF,d;
354: PetscInt i,M,xs,xm;
355: Vec xlocal;
358: DMGetLocalVector(da,&xlocal);
359: /*
360: Scatter ghost points to local vector, using the 2-step process
361: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
362: By placing code between these two statements, computations can
363: be done while messages are in transition.
364: */
365: DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
366: DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
368: /*
369: Get pointers to vector data.
370: - The vector xlocal includes ghost point; the vectors x and f do
371: NOT include ghost points.
372: - Using DMDAVecGetArray() allows accessing the values using global ordering
373: */
374: DMDAVecGetArray(da,xlocal,&xx);
375: DMDAVecGetArray(da,f,&ff);
376: DMDAVecGetArray(da,user->F,&FF);
378: /*
379: Get local grid boundaries (for 1-dimensional DMDA):
380: xs, xm - starting grid index, width of local grid (no ghost points)
381: */
382: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
383: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,
384: NULL,NULL,NULL,NULL,NULL);
386: /*
387: Set function values for boundary points; define local interior grid point range:
388: xsi - starting interior grid index
389: xei - ending interior grid index
390: */
391: if (xs == 0) { /* left boundary */
392: ff[0] = xx[0];
393: xs++;xm--;
394: }
395: if (xs+xm == M) { /* right boundary */
396: ff[xs+xm-1] = xx[xs+xm-1] - 1.0;
397: xm--;
398: }
400: /*
401: Compute function over locally owned part of the grid (interior points only)
402: */
403: d = 1.0/(user->h*user->h);
404: 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];
406: /*
407: Restore vectors
408: */
409: DMDAVecRestoreArray(da,xlocal,&xx);
410: DMDAVecRestoreArray(da,f,&ff);
411: DMDAVecRestoreArray(da,user->F,&FF);
412: DMRestoreLocalVector(da,&xlocal);
413: return(0);
414: }
415: /* ------------------------------------------------------------------- */
418: /*
419: FormJacobian - Evaluates Jacobian matrix.
421: Input Parameters:
422: . snes - the SNES context
423: . x - input vector
424: . dummy - optional user-defined context (not used here)
426: Output Parameters:
427: . jac - Jacobian matrix
428: . B - optionally different preconditioning matrix
429: . flag - flag indicating matrix structure
430: */
431: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
432: {
433: ApplicationCtx *user = (ApplicationCtx*) ctx;
434: PetscScalar *xx,d,A[3];
436: PetscInt i,j[3],M,xs,xm;
437: DM da = user->da;
440: /*
441: Get pointer to vector data
442: */
443: DMDAVecGetArray(da,x,&xx);
444: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
446: /*
447: Get range of locally owned matrix
448: */
449: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,
450: 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; A[0] = 1.0;
460: MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
461: xs++;xm--;
462: }
463: if (xs+xm == M) { /* right boundary */
464: i = M-1;
465: A[0] = 1.0;
466: 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; j[1] = i; j[2] = i + 1;
478: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
479: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
480: }
482: /*
483: Assemble matrix, using the 2-step process:
484: MatAssemblyBegin(), MatAssemblyEnd().
485: By placing code between these two statements, computations can be
486: done while messages are in transition.
488: Also, restore vector.
489: */
491: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
492: DMDAVecRestoreArray(da,x,&xx);
493: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
495: return(0);
496: }
497: /* ------------------------------------------------------------------- */
500: /*
501: Monitor - Optional user-defined monitoring routine that views the
502: current iterate with an x-window plot. Set by SNESMonitorSet().
504: Input Parameters:
505: snes - the SNES context
506: its - iteration number
507: norm - 2-norm function value (may be estimated)
508: ctx - optional user-defined context for private data for the
509: monitor routine, as set by SNESMonitorSet()
511: Note:
512: See the manpage for PetscViewerDrawOpen() for useful runtime options,
513: such as -nox to deactivate all x-window output.
514: */
515: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
516: {
518: MonitorCtx *monP = (MonitorCtx*) ctx;
519: Vec x;
522: PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %g\n",its,(double)fnorm);
523: SNESGetSolution(snes,&x);
524: VecView(x,monP->viewer);
525: return(0);
526: }
528: /* ------------------------------------------------------------------- */
531: /*
532: PreCheck - Optional user-defined routine that checks the validity of
533: candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
535: Input Parameters:
536: snes - the SNES context
537: xcurrent - current solution
538: y - search direction and length
540: Output Parameters:
541: y - proposed step (search direction and length) (possibly changed)
542: changed_y - tells if the step has changed or not
543: */
544: PetscErrorCode PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y, PetscBool *changed_y, void * ctx)
545: {
547: *changed_y = PETSC_FALSE;
548: return(0);
549: }
551: /* ------------------------------------------------------------------- */
554: /*
555: PostCheck - Optional user-defined routine that checks the validity of
556: candidate steps of a line search method. Set by SNESLineSearchSetPostCheck().
558: Input Parameters:
559: snes - the SNES context
560: ctx - optional user-defined context for private data for the
561: monitor routine, as set by SNESLineSearchSetPostCheck()
562: xcurrent - current solution
563: y - search direction and length
564: x - the new candidate iterate
566: Output Parameters:
567: y - proposed step (search direction and length) (possibly changed)
568: x - current iterate (possibly modified)
570: */
571: PetscErrorCode PostCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
572: {
574: PetscInt i,iter,xs,xm;
575: StepCheckCtx *check;
576: ApplicationCtx *user;
577: PetscScalar *xa,*xa_last,tmp;
578: PetscReal rdiff;
579: DM da;
580: SNES snes;
583: *changed_x = PETSC_FALSE;
584: *changed_y = PETSC_FALSE;
586: SNESLineSearchGetSNES(linesearch, &snes);
587: check = (StepCheckCtx*)ctx;
588: user = check->user;
589: SNESGetIterationNumber(snes,&iter);
590: SNESLineSearchGetPreCheck(linesearch, NULL, (void**)&check);
592: /* iteration 1 indicates we are working on the second iteration */
593: if (iter > 0) {
594: da = user->da;
595: PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %g\n",iter,(double)check->tolerance);
597: /* Access local array data */
598: DMDAVecGetArray(da,check->last_step,&xa_last);
599: DMDAVecGetArray(da,x,&xa);
600: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
602: /*
603: If we fail the user-defined check for validity of the candidate iterate,
604: then modify the iterate as we like. (Note that the particular modification
605: below is intended simply to demonstrate how to manipulate this data, not
606: as a meaningful or appropriate choice.)
607: */
608: for (i=xs; i<xs+xm; i++) {
609: if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance;
610: else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
611: if (rdiff > check->tolerance) {
612: tmp = xa[i];
613: xa[i] = .5*(xa[i] + xa_last[i]);
614: *changed_x = PETSC_TRUE;
615: PetscPrintf(PETSC_COMM_WORLD," Altering entry %D: x=%g, x_last=%g, diff=%g, x_new=%g\n",
616: i,(double)PetscAbsScalar(tmp),(double)PetscAbsScalar(xa_last[i]),(double)rdiff,(double)PetscAbsScalar(xa[i]));
617: }
618: }
619: DMDAVecRestoreArray(da,check->last_step,&xa_last);
620: DMDAVecRestoreArray(da,x,&xa);
621: }
622: VecCopy(x,check->last_step);
623: return(0);
624: }
627: /* ------------------------------------------------------------------- */
630: /*
631: PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
632: e.g,
633: 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
634: Set by SNESLineSearchSetPostCheck().
636: Input Parameters:
637: linesearch - the LineSearch context
638: xcurrent - current solution
639: y - search direction and length
640: x - the new candidate iterate
642: Output Parameters:
643: y - proposed step (search direction and length) (possibly changed)
644: x - current iterate (possibly modified)
646: */
647: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
648: {
650: SetSubKSPCtx *check;
651: PetscInt iter,its,sub_its,maxit;
652: KSP ksp,sub_ksp,*sub_ksps;
653: PC pc;
654: PetscReal ksp_ratio;
655: SNES snes;
658: SNESLineSearchGetSNES(linesearch, &snes);
659: check = (SetSubKSPCtx*)ctx;
660: SNESGetIterationNumber(snes,&iter);
661: SNESGetKSP(snes,&ksp);
662: KSPGetPC(ksp,&pc);
663: PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);
664: sub_ksp = sub_ksps[0];
665: KSPGetIterationNumber(ksp,&its); /* outer KSP iteration number */
666: KSPGetIterationNumber(sub_ksp,&sub_its); /* inner KSP iteration number */
668: if (iter) {
669: PetscPrintf(PETSC_COMM_WORLD," ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);
670: ksp_ratio = ((PetscReal)(its))/check->its0;
671: maxit = (PetscInt)(ksp_ratio*sub_its + 0.5);
672: if (maxit < 2) maxit = 2;
673: KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
674: PetscPrintf(PETSC_COMM_WORLD," ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);
675: }
676: check->its0 = its; /* save current outer KSP iteration number */
677: return(0);
678: }
680: /* ------------------------------------------------------------------- */
681: /*
682: MatrixFreePreconditioner - This routine demonstrates the use of a
683: user-provided preconditioner. This code implements just the null
684: preconditioner, which of course is not recommended for general use.
686: Input Parameters:
687: + pc - preconditioner
688: - x - input vector
690: Output Parameter:
691: . y - preconditioned vector
692: */
693: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
694: {
696: VecCopy(x,y);
697: return 0;
698: }