Actual source code: ex3.c
petsc-3.7.7 2017-09-25
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: -pre_check_iterates : activate checking of iterates\n\
8: -post_check_iterates : activate checking of iterates\n\
9: -check_tol <tol>: set tolerance for iterate checking\n\n";
11: /*T
12: Concepts: SNES^basic parallel example
13: Concepts: SNES^setting a user-defined monitoring routine
14: Concepts: error handling^using the macro __FUNCT__ to define routine names;
15: Processors: n
16: T*/
18: /*
19: Include "petscdraw.h" so that we can use distributed arrays (DMDAs).
20: Include "petscdraw.h" so that we can use PETSc drawing routines.
21: Include "petscsnes.h" so that we can use SNES solvers. Note that this
22: file automatically includes:
23: petscsys.h - base PETSc routines petscvec.h - vectors
24: petscmat.h - matrices
25: petscis.h - index sets petscksp.h - Krylov subspace methods
26: petscviewer.h - viewers petscpc.h - preconditioners
27: petscksp.h - linear solvers
28: */
30: #include <petscdm.h>
31: #include <petscdmda.h>
32: #include <petscsnes.h>
34: /*
35: User-defined routines. Note that immediately before each routine below,
36: we define the macro __FUNCT__ to be a string containing the routine name.
37: If defined, this macro is used in the PETSc error handlers to provide a
38: complete traceback of routine names. All PETSc library routines use this
39: macro, and users can optionally employ it as well in their application
40: codes. Note that users can get a traceback of PETSc errors regardless of
41: whether they define __FUNCT__ in application codes; this macro merely
42: provides the added traceback detail of the application routine names.
43: */
44: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
45: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
46: PetscErrorCode FormInitialGuess(Vec);
47: PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
48: PetscErrorCode PreCheck(SNESLineSearch,Vec,Vec,PetscBool*,void*);
49: PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
50: PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
51: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);
53: /*
54: User-defined application context
55: */
56: typedef struct {
57: DM da; /* distributed array */
58: Vec F; /* right-hand-side of PDE */
59: PetscMPIInt rank; /* rank of processor */
60: PetscMPIInt size; /* size of communicator */
61: PetscReal h; /* mesh spacing */
62: } ApplicationCtx;
64: /*
65: User-defined context for monitoring
66: */
67: typedef struct {
68: PetscViewer viewer;
69: } MonitorCtx;
71: /*
72: User-defined context for checking candidate iterates that are
73: determined by line search methods
74: */
75: typedef struct {
76: Vec last_step; /* previous iterate */
77: PetscReal tolerance; /* tolerance for changes between successive iterates */
78: ApplicationCtx *user;
79: } StepCheckCtx;
81: typedef struct {
82: PetscInt its0; /* num of prevous outer KSP iterations */
83: } SetSubKSPCtx;
87: int main(int argc,char **argv)
88: {
89: SNES snes; /* SNES context */
90: SNESLineSearch linesearch; /* SNESLineSearch context */
91: Mat J; /* Jacobian matrix */
92: ApplicationCtx ctx; /* user-defined context */
93: Vec x,r,U,F; /* vectors */
94: MonitorCtx monP; /* monitoring context */
95: StepCheckCtx checkP; /* step-checking context */
96: SetSubKSPCtx checkP1;
97: PetscBool pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking candidate iterates */
98: PetscScalar xp,*FF,*UU,none = -1.0;
100: PetscInt its,N = 5,i,maxit,maxf,xs,xm;
101: PetscReal abstol,rtol,stol,norm;
102: PetscBool flg;
105: PetscInitialize(&argc,&argv,(char*)0,help);
106: MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
107: MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
108: PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
109: ctx.h = 1.0/(N-1);
111: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112: Create nonlinear solver context
113: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: SNESCreate(PETSC_COMM_WORLD,&snes);
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Create vector data structures; set function evaluation routine
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: /*
122: Create distributed array (DMDA) to manage parallel grid and vectors
123: */
124: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,1,1,NULL,&ctx.da);
126: /*
127: Extract global and local vectors from DMDA; then duplicate for remaining
128: vectors that are the same types
129: */
130: DMCreateGlobalVector(ctx.da,&x);
131: VecDuplicate(x,&r);
132: VecDuplicate(x,&F); ctx.F = F;
133: VecDuplicate(x,&U);
135: /*
136: Set function evaluation routine and vector. Whenever the nonlinear
137: solver needs to compute the nonlinear function, it will call this
138: routine.
139: - Note that the final routine argument is the user-defined
140: context that provides application-specific data for the
141: function evaluation routine.
142: */
143: SNESSetFunction(snes,r,FormFunction,&ctx);
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Create matrix data structure; set Jacobian evaluation routine
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: MatCreate(PETSC_COMM_WORLD,&J);
150: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
151: MatSetFromOptions(J);
152: MatSeqAIJSetPreallocation(J,3,NULL);
153: MatMPIAIJSetPreallocation(J,3,NULL,3,NULL);
155: /*
156: Set Jacobian matrix data structure and default Jacobian evaluation
157: routine. Whenever the nonlinear solver needs to compute the
158: Jacobian matrix, it will call this routine.
159: - Note that the final routine argument is the user-defined
160: context that provides application-specific data for the
161: Jacobian evaluation routine.
162: */
163: SNESSetJacobian(snes,J,J,FormJacobian,&ctx);
165: /*
166: Optional allow user provided preconditioner
167: */
168: PetscOptionsHasName(NULL,NULL,"-user_precond",&flg);
169: if (flg) {
170: KSP ksp;
171: PC pc;
172: SNESGetKSP(snes,&ksp);
173: KSPGetPC(ksp,&pc);
174: PCSetType(pc,PCSHELL);
175: PCShellSetApply(pc,MatrixFreePreconditioner);
176: }
178: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
179: Customize nonlinear solver; set runtime options
180: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182: /*
183: Set an optional user-defined monitoring routine
184: */
185: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
186: SNESMonitorSet(snes,Monitor,&monP,0);
188: /*
189: Set names for some vectors to facilitate monitoring (optional)
190: */
191: PetscObjectSetName((PetscObject)x,"Approximate Solution");
192: PetscObjectSetName((PetscObject)U,"Exact Solution");
194: /*
195: Set SNES/KSP/KSP/PC runtime options, e.g.,
196: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
197: */
198: SNESSetFromOptions(snes);
200: /*
201: Set an optional user-defined routine to check the validity of candidate
202: iterates that are determined by line search methods
203: */
204: SNESGetLineSearch(snes, &linesearch);
205: PetscOptionsHasName(NULL,NULL,"-post_check_iterates",&post_check);
207: if (post_check) {
208: PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
209: SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP);
210: VecDuplicate(x,&(checkP.last_step));
212: checkP.tolerance = 1.0;
213: checkP.user = &ctx;
215: PetscOptionsGetReal(NULL,NULL,"-check_tol",&checkP.tolerance,NULL);
216: }
218: PetscOptionsHasName(NULL,NULL,"-post_setsubksp",&post_setsubksp);
219: if (post_setsubksp) {
220: PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");
221: SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1);
222: }
224: PetscOptionsHasName(NULL,NULL,"-pre_check_iterates",&pre_check);
225: if (pre_check) {
226: PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
227: SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP);
228: }
230: /*
231: Print parameters used for convergence testing (optional) ... just
232: to demonstrate this routine; this information is also printed with
233: the option -snes_view
234: */
235: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
236: PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);
238: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239: Initialize application:
240: Store right-hand-side of PDE and exact solution
241: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
243: /*
244: Get local grid boundaries (for 1-dimensional DMDA):
245: xs, xm - starting grid index, width of local grid (no ghost points)
246: */
247: DMDAGetCorners(ctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
249: /*
250: Get pointers to vector data
251: */
252: DMDAVecGetArray(ctx.da,F,&FF);
253: DMDAVecGetArray(ctx.da,U,&UU);
255: /*
256: Compute local vector entries
257: */
258: xp = ctx.h*xs;
259: for (i=xs; i<xs+xm; i++) {
260: FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
261: UU[i] = xp*xp*xp;
262: xp += ctx.h;
263: }
265: /*
266: Restore vectors
267: */
268: DMDAVecRestoreArray(ctx.da,F,&FF);
269: DMDAVecRestoreArray(ctx.da,U,&UU);
271: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
272: Evaluate initial guess; then solve nonlinear system
273: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
275: /*
276: Note: The user should initialize the vector, x, with the initial guess
277: for the nonlinear solver prior to calling SNESSolve(). In particular,
278: to employ an initial guess of zero, the user should explicitly set
279: this vector to zero by calling VecSet().
280: */
281: FormInitialGuess(x);
282: SNESSolve(snes,NULL,x);
283: SNESGetIterationNumber(snes,&its);
284: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
286: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
287: Check solution and clean up
288: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
290: /*
291: Check the error
292: */
293: VecAXPY(x,none,U);
294: VecNorm(x,NORM_2,&norm);
295: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)norm,its);
297: /*
298: Free work space. All PETSc objects should be destroyed when they
299: are no longer needed.
300: */
301: PetscViewerDestroy(&monP.viewer);
302: if (post_check) {VecDestroy(&checkP.last_step);}
303: VecDestroy(&x);
304: VecDestroy(&r);
305: VecDestroy(&U);
306: VecDestroy(&F);
307: MatDestroy(&J);
308: SNESDestroy(&snes);
309: DMDestroy(&ctx.da);
310: PetscFinalize();
311: return(0);
312: }
313: /* ------------------------------------------------------------------- */
316: /*
317: FormInitialGuess - Computes initial guess.
319: Input/Output Parameter:
320: . x - the solution vector
321: */
322: PetscErrorCode FormInitialGuess(Vec x)
323: {
325: PetscScalar pfive = .50;
328: VecSet(x,pfive);
329: return(0);
330: }
331: /* ------------------------------------------------------------------- */
334: /*
335: FormFunction - Evaluates nonlinear function, F(x).
337: Input Parameters:
338: . snes - the SNES context
339: . x - input vector
340: . ctx - optional user-defined context, as set by SNESSetFunction()
342: Output Parameter:
343: . f - function vector
345: Note:
346: The user-defined context can contain any application-specific
347: data needed for the function evaluation.
348: */
349: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
350: {
351: ApplicationCtx *user = (ApplicationCtx*) ctx;
352: DM da = user->da;
353: PetscScalar *xx,*ff,*FF,d;
355: PetscInt i,M,xs,xm;
356: Vec xlocal;
359: DMGetLocalVector(da,&xlocal);
360: /*
361: Scatter ghost points to local vector, using the 2-step process
362: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
363: By placing code between these two statements, computations can
364: be done while messages are in transition.
365: */
366: DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
367: DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
369: /*
370: Get pointers to vector data.
371: - The vector xlocal includes ghost point; the vectors x and f do
372: NOT include ghost points.
373: - Using DMDAVecGetArray() allows accessing the values using global ordering
374: */
375: DMDAVecGetArray(da,xlocal,&xx);
376: DMDAVecGetArray(da,f,&ff);
377: DMDAVecGetArray(da,user->F,&FF);
379: /*
380: Get local grid boundaries (for 1-dimensional DMDA):
381: xs, xm - starting grid index, width of local grid (no ghost points)
382: */
383: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
384: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,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: DMDAVecGetArrayRead(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,NULL,NULL,NULL,NULL,NULL);
451: /*
452: Determine starting and ending local indices for interior grid points.
453: Set Jacobian entries for boundary points.
454: */
456: if (xs == 0) { /* left boundary */
457: i = 0; A[0] = 1.0;
459: MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
460: xs++;xm--;
461: }
462: if (xs+xm == M) { /* right boundary */
463: i = M-1;
464: A[0] = 1.0;
465: MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
466: xm--;
467: }
469: /*
470: Interior grid points
471: - Note that in this case we set all elements for a particular
472: row at once.
473: */
474: d = 1.0/(user->h*user->h);
475: for (i=xs; i<xs+xm; i++) {
476: j[0] = i - 1; j[1] = i; j[2] = i + 1;
477: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
478: MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
479: }
481: /*
482: Assemble matrix, using the 2-step process:
483: MatAssemblyBegin(), MatAssemblyEnd().
484: By placing code between these two statements, computations can be
485: done while messages are in transition.
487: Also, restore vector.
488: */
490: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
491: DMDAVecRestoreArrayRead(da,x,&xx);
492: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
494: return(0);
495: }
496: /* ------------------------------------------------------------------- */
499: /*
500: Monitor - Optional user-defined monitoring routine that views the
501: current iterate with an x-window plot. Set by SNESMonitorSet().
503: Input Parameters:
504: snes - the SNES context
505: its - iteration number
506: norm - 2-norm function value (may be estimated)
507: ctx - optional user-defined context for private data for the
508: monitor routine, as set by SNESMonitorSet()
510: Note:
511: See the manpage for PetscViewerDrawOpen() for useful runtime options,
512: such as -nox to deactivate all x-window output.
513: */
514: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
515: {
517: MonitorCtx *monP = (MonitorCtx*) ctx;
518: Vec x;
521: PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %g\n",its,(double)fnorm);
522: SNESGetSolution(snes,&x);
523: VecView(x,monP->viewer);
524: return(0);
525: }
527: /* ------------------------------------------------------------------- */
530: /*
531: PreCheck - Optional user-defined routine that checks the validity of
532: candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
534: Input Parameters:
535: snes - the SNES context
536: xcurrent - current solution
537: y - search direction and length
539: Output Parameters:
540: y - proposed step (search direction and length) (possibly changed)
541: changed_y - tells if the step has changed or not
542: */
543: PetscErrorCode PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y, PetscBool *changed_y, void * ctx)
544: {
546: *changed_y = PETSC_FALSE;
547: return(0);
548: }
550: /* ------------------------------------------------------------------- */
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: {
573: PetscInt i,iter,xs,xm;
574: StepCheckCtx *check;
575: ApplicationCtx *user;
576: PetscScalar *xa,*xa_last,tmp;
577: PetscReal rdiff;
578: DM da;
579: SNES snes;
582: *changed_x = PETSC_FALSE;
583: *changed_y = PETSC_FALSE;
585: SNESLineSearchGetSNES(linesearch, &snes);
586: check = (StepCheckCtx*)ctx;
587: user = check->user;
588: SNESGetIterationNumber(snes,&iter);
590: /* iteration 1 indicates we are working on the second iteration */
591: if (iter > 0) {
592: da = user->da;
593: PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %g\n",iter,(double)check->tolerance);
595: /* Access local array data */
596: DMDAVecGetArray(da,check->last_step,&xa_last);
597: DMDAVecGetArray(da,x,&xa);
598: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
600: /*
601: If we fail the user-defined check for validity of the candidate iterate,
602: then modify the iterate as we like. (Note that the particular modification
603: below is intended simply to demonstrate how to manipulate this data, not
604: as a meaningful or appropriate choice.)
605: */
606: for (i=xs; i<xs+xm; i++) {
607: if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance;
608: else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
609: if (rdiff > check->tolerance) {
610: tmp = xa[i];
611: xa[i] = .5*(xa[i] + xa_last[i]);
612: *changed_x = PETSC_TRUE;
613: PetscPrintf(PETSC_COMM_WORLD," Altering entry %D: x=%g, x_last=%g, diff=%g, x_new=%g\n",
614: i,(double)PetscAbsScalar(tmp),(double)PetscAbsScalar(xa_last[i]),(double)rdiff,(double)PetscAbsScalar(xa[i]));
615: }
616: }
617: DMDAVecRestoreArray(da,check->last_step,&xa_last);
618: DMDAVecRestoreArray(da,x,&xa);
619: }
620: VecCopy(x,check->last_step);
621: return(0);
622: }
625: /* ------------------------------------------------------------------- */
628: /*
629: PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
630: e.g,
631: 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
632: Set by SNESLineSearchSetPostCheck().
634: Input Parameters:
635: linesearch - the LineSearch context
636: xcurrent - current solution
637: y - search direction and length
638: x - the new candidate iterate
640: Output Parameters:
641: y - proposed step (search direction and length) (possibly changed)
642: x - current iterate (possibly modified)
644: */
645: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
646: {
648: SetSubKSPCtx *check;
649: PetscInt iter,its,sub_its,maxit;
650: KSP ksp,sub_ksp,*sub_ksps;
651: PC pc;
652: PetscReal ksp_ratio;
653: SNES snes;
656: SNESLineSearchGetSNES(linesearch, &snes);
657: check = (SetSubKSPCtx*)ctx;
658: SNESGetIterationNumber(snes,&iter);
659: SNESGetKSP(snes,&ksp);
660: KSPGetPC(ksp,&pc);
661: PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);
662: sub_ksp = sub_ksps[0];
663: KSPGetIterationNumber(ksp,&its); /* outer KSP iteration number */
664: KSPGetIterationNumber(sub_ksp,&sub_its); /* inner KSP iteration number */
666: if (iter) {
667: PetscPrintf(PETSC_COMM_WORLD," ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);
668: ksp_ratio = ((PetscReal)(its))/check->its0;
669: maxit = (PetscInt)(ksp_ratio*sub_its + 0.5);
670: if (maxit < 2) maxit = 2;
671: KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
672: PetscPrintf(PETSC_COMM_WORLD," ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);
673: }
674: check->its0 = its; /* save current outer KSP iteration number */
675: return(0);
676: }
678: /* ------------------------------------------------------------------- */
679: /*
680: MatrixFreePreconditioner - This routine demonstrates the use of a
681: user-provided preconditioner. This code implements just the null
682: preconditioner, which of course is not recommended for general use.
684: Input Parameters:
685: + pc - preconditioner
686: - x - input vector
688: Output Parameter:
689: . y - preconditioned vector
690: */
691: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
692: {
694: VecCopy(x,y);
695: return 0;
696: }