Actual source code: ex3.c
petsc-3.4.5 2014-06-29
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(SNESLineSearch,Vec,Vec,PetscBool*,void*);
47: PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
48: PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
49: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);
51: /*
52: User-defined application context
53: */
54: typedef struct {
55: DM da; /* distributed array */
56: Vec F; /* right-hand-side of PDE */
57: PetscMPIInt rank; /* rank of processor */
58: PetscMPIInt size; /* size of communicator */
59: PetscReal h; /* mesh spacing */
60: } ApplicationCtx;
62: /*
63: User-defined context for monitoring
64: */
65: typedef struct {
66: PetscViewer viewer;
67: } MonitorCtx;
69: /*
70: User-defined context for checking candidate iterates that are
71: determined by line search methods
72: */
73: typedef struct {
74: Vec last_step; /* previous iterate */
75: PetscReal tolerance; /* tolerance for changes between successive iterates */
76: ApplicationCtx *user;
77: } StepCheckCtx;
79: typedef struct {
80: PetscInt its0; /* num of prevous outer KSP iterations */
81: } SetSubKSPCtx;
85: int main(int argc,char **argv)
86: {
87: SNES snes; /* SNES context */
88: SNESLineSearch linesearch; /* SNESLineSearch context */
89: Mat J; /* Jacobian matrix */
90: ApplicationCtx ctx; /* user-defined context */
91: Vec x,r,U,F; /* vectors */
92: MonitorCtx monP; /* monitoring context */
93: StepCheckCtx checkP; /* step-checking context */
94: SetSubKSPCtx checkP1;
95: PetscBool pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking candidate iterates */
96: PetscScalar xp,*FF,*UU,none = -1.0;
98: PetscInt its,N = 5,i,maxit,maxf,xs,xm;
99: PetscReal abstol,rtol,stol,norm;
100: PetscBool flg;
103: PetscInitialize(&argc,&argv,(char*)0,help);
104: MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
105: MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
106: PetscOptionsGetInt(NULL,"-n",&N,NULL);
107: ctx.h = 1.0/(N-1);
109: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110: Create nonlinear solver context
111: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: SNESCreate(PETSC_COMM_WORLD,&snes);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Create vector data structures; set function evaluation routine
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: /*
120: Create distributed array (DMDA) to manage parallel grid and vectors
121: */
122: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,N,1,1,NULL,&ctx.da);
124: /*
125: Extract global and local vectors from DMDA; then duplicate for remaining
126: vectors that are the same types
127: */
128: DMCreateGlobalVector(ctx.da,&x);
129: VecDuplicate(x,&r);
130: VecDuplicate(x,&F); ctx.F = F;
131: VecDuplicate(x,&U);
133: /*
134: Set function evaluation routine and vector. Whenever the nonlinear
135: solver needs to compute the nonlinear function, it will call this
136: routine.
137: - Note that the final routine argument is the user-defined
138: context that provides application-specific data for the
139: function evaluation routine.
140: */
141: SNESSetFunction(snes,r,FormFunction,&ctx);
143: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: Create matrix data structure; set Jacobian evaluation routine
145: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147: MatCreate(PETSC_COMM_WORLD,&J);
148: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
149: MatSetFromOptions(J);
150: MatSeqAIJSetPreallocation(J,3,NULL);
151: MatMPIAIJSetPreallocation(J,3,NULL,3,NULL);
153: /*
154: Set Jacobian matrix data structure and default Jacobian evaluation
155: routine. Whenever the nonlinear solver needs to compute the
156: Jacobian matrix, it will call this routine.
157: - Note that the final routine argument is the user-defined
158: context that provides application-specific data for the
159: Jacobian evaluation routine.
160: */
161: SNESSetJacobian(snes,J,J,FormJacobian,&ctx);
163: /*
164: Optional allow user provided preconditioner
165: */
166: PetscOptionsHasName(NULL,"-user_precond",&flg);
167: if (flg) {
168: KSP ksp;
169: PC pc;
170: SNESGetKSP(snes,&ksp);
171: KSPGetPC(ksp,&pc);
172: PCSetType(pc,PCSHELL);
173: PCShellSetApply(pc,MatrixFreePreconditioner);
174: }
176: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177: Customize nonlinear solver; set runtime options
178: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180: /*
181: Set an optional user-defined monitoring routine
182: */
183: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
184: SNESMonitorSet(snes,Monitor,&monP,0);
186: /*
187: Set names for some vectors to facilitate monitoring (optional)
188: */
189: PetscObjectSetName((PetscObject)x,"Approximate Solution");
190: PetscObjectSetName((PetscObject)U,"Exact Solution");
192: /*
193: Set SNES/KSP/KSP/PC runtime options, e.g.,
194: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
195: */
196: SNESSetFromOptions(snes);
198: /*
199: Set an optional user-defined routine to check the validity of candidate
200: iterates that are determined by line search methods
201: */
202: SNESGetLineSearch(snes, &linesearch);
203: PetscOptionsHasName(NULL,"-post_check_iterates",&post_check);
205: if (post_check) {
206: PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
207: SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP);
208: VecDuplicate(x,&(checkP.last_step));
210: checkP.tolerance = 1.0;
211: checkP.user = &ctx;
213: PetscOptionsGetReal(NULL,"-check_tol",&checkP.tolerance,NULL);
214: }
216: PetscOptionsHasName(NULL,"-post_setsubksp",&post_setsubksp);
217: if (post_setsubksp) {
218: PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");
219: SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1);
220: }
222: PetscOptionsHasName(NULL,"-pre_check_iterates",&pre_check);
223: if (pre_check) {
224: PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
225: SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP);
226: }
228: /*
229: Print parameters used for convergence testing (optional) ... just
230: to demonstrate this routine; this information is also printed with
231: the option -snes_view
232: */
233: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
234: PetscPrintf(PETSC_COMM_WORLD,"atol=%G, rtol=%G, stol=%G, maxit=%D, maxf=%D\n",abstol,rtol,stol,maxit,maxf);
236: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237: Initialize application:
238: Store right-hand-side of PDE and exact solution
239: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
241: /*
242: Get local grid boundaries (for 1-dimensional DMDA):
243: xs, xm - starting grid index, width of local grid (no ghost points)
244: */
245: DMDAGetCorners(ctx.da,&xs,NULL,NULL,&xm,NULL,NULL);
247: /*
248: Get pointers to vector data
249: */
250: DMDAVecGetArray(ctx.da,F,&FF);
251: DMDAVecGetArray(ctx.da,U,&UU);
253: /*
254: Compute local vector entries
255: */
256: xp = ctx.h*xs;
257: for (i=xs; i<xs+xm; i++) {
258: FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
259: UU[i] = xp*xp*xp;
260: xp += ctx.h;
261: }
263: /*
264: Restore vectors
265: */
266: DMDAVecRestoreArray(ctx.da,F,&FF);
267: DMDAVecRestoreArray(ctx.da,U,&UU);
269: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270: Evaluate initial guess; then solve nonlinear system
271: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
273: /*
274: Note: The user should initialize the vector, x, with the initial guess
275: for the nonlinear solver prior to calling SNESSolve(). In particular,
276: to employ an initial guess of zero, the user should explicitly set
277: this vector to zero by calling VecSet().
278: */
279: FormInitialGuess(x);
280: SNESSolve(snes,NULL,x);
281: SNESGetIterationNumber(snes,&its);
282: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n\n",its);
284: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
285: Check solution and clean up
286: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
288: /*
289: Check the error
290: */
291: VecAXPY(x,none,U);
292: VecNorm(x,NORM_2,&norm);
293: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G Iterations %D\n",norm,its);
295: /*
296: Free work space. All PETSc objects should be destroyed when they
297: are no longer needed.
298: */
299: PetscViewerDestroy(&monP.viewer);
300: if (post_check) {VecDestroy(&checkP.last_step);}
301: VecDestroy(&x);
302: VecDestroy(&r);
303: VecDestroy(&U);
304: VecDestroy(&F);
305: MatDestroy(&J);
306: SNESDestroy(&snes);
307: DMDestroy(&ctx.da);
308: PetscFinalize();
309: return(0);
310: }
311: /* ------------------------------------------------------------------- */
314: /*
315: FormInitialGuess - Computes initial guess.
317: Input/Output Parameter:
318: . x - the solution vector
319: */
320: PetscErrorCode FormInitialGuess(Vec x)
321: {
323: PetscScalar pfive = .50;
326: VecSet(x,pfive);
327: return(0);
328: }
329: /* ------------------------------------------------------------------- */
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 *xx,*ff,*FF,d;
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: DMDAVecGetArray(da,xlocal,&xx);
374: DMDAVecGetArray(da,f,&ff);
375: DMDAVecGetArray(da,user->F,&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,
383: NULL,NULL,NULL,NULL,NULL);
385: /*
386: Set function values for boundary points; define local interior grid point range:
387: xsi - starting interior grid index
388: xei - ending interior grid index
389: */
390: if (xs == 0) { /* left boundary */
391: ff[0] = xx[0];
392: xs++;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: DMDAVecRestoreArray(da,xlocal,&xx);
409: DMDAVecRestoreArray(da,f,&ff);
410: DMDAVecRestoreArray(da,user->F,&FF);
411: DMRestoreLocalVector(da,&xlocal);
412: return(0);
413: }
414: /* ------------------------------------------------------------------- */
417: /*
418: FormJacobian - Evaluates Jacobian matrix.
420: Input Parameters:
421: . snes - the SNES context
422: . x - input vector
423: . dummy - optional user-defined context (not used here)
425: Output Parameters:
426: . jac - Jacobian matrix
427: . B - optionally different preconditioning matrix
428: . flag - flag indicating matrix structure
429: */
430: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure *flag,void *ctx)
431: {
432: ApplicationCtx *user = (ApplicationCtx*) ctx;
433: PetscScalar *xx,d,A[3];
435: PetscInt i,j[3],M,xs,xm;
436: DM da = user->da;
439: /*
440: Get pointer to vector data
441: */
442: DMDAVecGetArray(da,x,&xx);
443: DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
445: /*
446: Get range of locally owned matrix
447: */
448: DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,
449: 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: DMDAVecRestoreArray(da,x,&xx);
492: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
494: *flag = SAME_NONZERO_PATTERN;
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,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,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,PetscAbsScalar(tmp),PetscAbsScalar(xa_last[i]),rdiff,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: }