Actual source code: ex3.c
petsc-3.3-p7 2013-05-11
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*);
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: SNESLineSearch linesearch; /* SNESLineSearch context */
88: Mat J; /* Jacobian matrix */
89: ApplicationCtx ctx; /* user-defined context */
90: Vec x,r,U,F; /* vectors */
91: MonitorCtx monP; /* monitoring context */
92: StepCheckCtx checkP; /* step-checking context */
93: SetSubKSPCtx checkP1;
94: PetscBool pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking
95: 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;
101: PetscInitialize(&argc,&argv,(char *)0,help);
102: MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
103: MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
104: PetscOptionsGetInt(PETSC_NULL,"-n",&N,PETSC_NULL);
105: ctx.h = 1.0/(N-1);
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Create nonlinear solver context
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: SNESCreate(PETSC_COMM_WORLD,&snes);
113: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
114: Create vector data structures; set function evaluation routine
115: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
117: /*
118: Create distributed array (DMDA) to manage parallel grid and vectors
119: */
120: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,N,1,1,PETSC_NULL,&ctx.da);
122: /*
123: Extract global and local vectors from DMDA; then duplicate for remaining
124: vectors that are the same types
125: */
126: DMCreateGlobalVector(ctx.da,&x);
127: VecDuplicate(x,&r);
128: VecDuplicate(x,&F); ctx.F = F;
129: VecDuplicate(x,&U);
131: /*
132: Set function evaluation routine and vector. Whenever the nonlinear
133: solver needs to compute the nonlinear function, it will call this
134: routine.
135: - Note that the final routine argument is the user-defined
136: context that provides application-specific data for the
137: function evaluation routine.
138: */
139: SNESSetFunction(snes,r,FormFunction,&ctx);
141: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: Create matrix data structure; set Jacobian evaluation routine
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: MatCreate(PETSC_COMM_WORLD,&J);
146: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
147: MatSetFromOptions(J);
148: MatSeqAIJSetPreallocation(J,3,PETSC_NULL);
149: MatMPIAIJSetPreallocation(J,3,PETSC_NULL,3,PETSC_NULL);
151: /*
152: Set Jacobian matrix data structure and default Jacobian evaluation
153: routine. Whenever the nonlinear solver needs to compute the
154: Jacobian matrix, it will call this routine.
155: - Note that the final routine argument is the user-defined
156: context that provides application-specific data for the
157: Jacobian evaluation routine.
158: */
159: SNESSetJacobian(snes,J,J,FormJacobian,&ctx);
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: Customize nonlinear solver; set runtime options
163: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: /*
166: Set an optional user-defined monitoring routine
167: */
168: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
169: SNESMonitorSet(snes,Monitor,&monP,0);
171: /*
172: Set names for some vectors to facilitate monitoring (optional)
173: */
174: PetscObjectSetName((PetscObject)x,"Approximate Solution");
175: PetscObjectSetName((PetscObject)U,"Exact Solution");
177: /*
178: Set SNES/KSP/KSP/PC runtime options, e.g.,
179: -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
180: */
181: SNESSetFromOptions(snes);
183: /*
184: Set an optional user-defined routine to check the validity of candidate
185: iterates that are determined by line search methods
186: */
187: SNESGetSNESLineSearch(snes, &linesearch);
188: PetscOptionsHasName(PETSC_NULL,"-post_check_iterates",&post_check);
190: if (post_check) {
191: PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
192: SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP);
193: VecDuplicate(x,&(checkP.last_step));
194: checkP.tolerance = 1.0;
195: checkP.user = &ctx;
196: PetscOptionsGetReal(PETSC_NULL,"-check_tol",&checkP.tolerance,PETSC_NULL);
197: }
199: PetscOptionsHasName(PETSC_NULL,"-post_setsubksp",&post_setsubksp);
200: if (post_setsubksp) {
201: PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");
202: SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1);
203: }
205: PetscOptionsHasName(PETSC_NULL,"-pre_check_iterates",&pre_check);
206: if (pre_check) {
207: PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
208: SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP);
209: }
211: /*
212: Print parameters used for convergence testing (optional) ... just
213: to demonstrate this routine; this information is also printed with
214: the option -snes_view
215: */
216: SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
217: PetscPrintf(PETSC_COMM_WORLD,"atol=%G, rtol=%G, stol=%G, maxit=%D, maxf=%D\n",abstol,rtol,stol,maxit,maxf);
219: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
220: Initialize application:
221: Store right-hand-side of PDE and exact solution
222: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
224: /*
225: Get local grid boundaries (for 1-dimensional DMDA):
226: xs, xm - starting grid index, width of local grid (no ghost points)
227: */
228: DMDAGetCorners(ctx.da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
230: /*
231: Get pointers to vector data
232: */
233: DMDAVecGetArray(ctx.da,F,&FF);
234: DMDAVecGetArray(ctx.da,U,&UU);
236: /*
237: Compute local vector entries
238: */
239: xp = ctx.h*xs;
240: for (i=xs; i<xs+xm; i++) {
241: FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
242: UU[i] = xp*xp*xp;
243: xp += ctx.h;
244: }
246: /*
247: Restore vectors
248: */
249: DMDAVecRestoreArray(ctx.da,F,&FF);
250: DMDAVecRestoreArray(ctx.da,U,&UU);
252: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253: Evaluate initial guess; then solve nonlinear system
254: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
256: /*
257: Note: The user should initialize the vector, x, with the initial guess
258: for the nonlinear solver prior to calling SNESSolve(). In particular,
259: to employ an initial guess of zero, the user should explicitly set
260: this vector to zero by calling VecSet().
261: */
262: FormInitialGuess(x);
263: SNESSolve(snes,PETSC_NULL,x);
264: SNESGetIterationNumber(snes,&its);
265: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n\n",its);
267: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268: Check solution and clean up
269: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271: /*
272: Check the error
273: */
274: VecAXPY(x,none,U);
275: VecNorm(x,NORM_2,&norm);
276: PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G Iterations %D\n",norm,its);
278: /*
279: Free work space. All PETSc objects should be destroyed when they
280: are no longer needed.
281: */
282: PetscViewerDestroy(&monP.viewer);
283: if (post_check) {VecDestroy(&checkP.last_step);}
284: VecDestroy(&x);
285: VecDestroy(&r);
286: VecDestroy(&U);
287: VecDestroy(&F);
288: MatDestroy(&J);
289: SNESDestroy(&snes);
290: DMDestroy(&ctx.da);
291: PetscFinalize();
293: return(0);
294: }
295: /* ------------------------------------------------------------------- */
298: /*
299: FormInitialGuess - Computes initial guess.
301: Input/Output Parameter:
302: . x - the solution vector
303: */
304: PetscErrorCode FormInitialGuess(Vec x)
305: {
307: PetscScalar pfive = .50;
310: VecSet(x,pfive);
311: return(0);
312: }
313: /* ------------------------------------------------------------------- */
316: /*
317: FormFunction - Evaluates nonlinear function, F(x).
319: Input Parameters:
320: . snes - the SNES context
321: . x - input vector
322: . ctx - optional user-defined context, as set by SNESSetFunction()
324: Output Parameter:
325: . f - function vector
327: Note:
328: The user-defined context can contain any application-specific
329: data needed for the function evaluation.
330: */
331: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
332: {
333: ApplicationCtx *user = (ApplicationCtx*) ctx;
334: DM da = user->da;
335: PetscScalar *xx,*ff,*FF,d;
337: PetscInt i,M,xs,xm;
338: Vec xlocal;
341: DMGetLocalVector(da,&xlocal);
342: /*
343: Scatter ghost points to local vector, using the 2-step process
344: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
345: By placing code between these two statements, computations can
346: be done while messages are in transition.
347: */
348: DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
349: DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);
351: /*
352: Get pointers to vector data.
353: - The vector xlocal includes ghost point; the vectors x and f do
354: NOT include ghost points.
355: - Using DMDAVecGetArray() allows accessing the values using global ordering
356: */
357: DMDAVecGetArray(da,xlocal,&xx);
358: DMDAVecGetArray(da,f,&ff);
359: DMDAVecGetArray(da,user->F,&FF);
361: /*
362: Get local grid boundaries (for 1-dimensional DMDA):
363: xs, xm - starting grid index, width of local grid (no ghost points)
364: */
365: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
366: DMDAGetInfo(da,PETSC_NULL,&M,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
367: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
369: /*
370: Set function values for boundary points; define local interior grid point range:
371: xsi - starting interior grid index
372: xei - ending interior grid index
373: */
374: if (xs == 0) { /* left boundary */
375: ff[0] = xx[0];
376: xs++;xm--;
377: }
378: if (xs+xm == M) { /* right boundary */
379: ff[xs+xm-1] = xx[xs+xm-1] - 1.0;
380: xm--;
381: }
383: /*
384: Compute function over locally owned part of the grid (interior points only)
385: */
386: d = 1.0/(user->h*user->h);
387: for (i=xs; i<xs+xm; i++) {
388: ff[i] = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - FF[i];
389: }
391: /*
392: Restore vectors
393: */
394: DMDAVecRestoreArray(da,xlocal,&xx);
395: DMDAVecRestoreArray(da,f,&ff);
396: DMDAVecRestoreArray(da,user->F,&FF);
397: DMRestoreLocalVector(da,&xlocal);
398: return(0);
399: }
400: /* ------------------------------------------------------------------- */
403: /*
404: FormJacobian - Evaluates Jacobian matrix.
406: Input Parameters:
407: . snes - the SNES context
408: . x - input vector
409: . dummy - optional user-defined context (not used here)
411: Output Parameters:
412: . jac - Jacobian matrix
413: . B - optionally different preconditioning matrix
414: . flag - flag indicating matrix structure
415: */
416: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *ctx)
417: {
418: ApplicationCtx *user = (ApplicationCtx*) ctx;
419: PetscScalar *xx,d,A[3];
421: PetscInt i,j[3],M,xs,xm;
422: DM da = user->da;
425: /*
426: Get pointer to vector data
427: */
428: DMDAVecGetArray(da,x,&xx);
429: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
431: /*
432: Get range of locally owned matrix
433: */
434: DMDAGetInfo(da,PETSC_NULL,&M,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
435: PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL);
437: /*
438: Determine starting and ending local indices for interior grid points.
439: Set Jacobian entries for boundary points.
440: */
442: if (xs == 0) { /* left boundary */
443: i = 0; A[0] = 1.0;
444: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
445: xs++;xm--;
446: }
447: if (xs+xm == M) { /* right boundary */
448: i = M-1; A[0] = 1.0;
449: MatSetValues(*jac,1,&i,1,&i,A,INSERT_VALUES);
450: xm--;
451: }
453: /*
454: Interior grid points
455: - Note that in this case we set all elements for a particular
456: row at once.
457: */
458: d = 1.0/(user->h*user->h);
459: for (i=xs; i<xs+xm; i++) {
460: j[0] = i - 1; j[1] = i; j[2] = i + 1;
461: A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
462: MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);
463: }
465: /*
466: Assemble matrix, using the 2-step process:
467: MatAssemblyBegin(), MatAssemblyEnd().
468: By placing code between these two statements, computations can be
469: done while messages are in transition.
471: Also, restore vector.
472: */
474: MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);
475: DMDAVecRestoreArray(da,x,&xx);
476: MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);
478: *flag = SAME_NONZERO_PATTERN;
479: return(0);
480: }
481: /* ------------------------------------------------------------------- */
484: /*
485: Monitor - Optional user-defined monitoring routine that views the
486: current iterate with an x-window plot. Set by SNESMonitorSet().
488: Input Parameters:
489: snes - the SNES context
490: its - iteration number
491: norm - 2-norm function value (may be estimated)
492: ctx - optional user-defined context for private data for the
493: monitor routine, as set by SNESMonitorSet()
495: Note:
496: See the manpage for PetscViewerDrawOpen() for useful runtime options,
497: such as -nox to deactivate all x-window output.
498: */
499: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
500: {
502: MonitorCtx *monP = (MonitorCtx*) ctx;
503: Vec x;
506: PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %G\n",its,fnorm);
507: SNESGetSolution(snes,&x);
508: VecView(x,monP->viewer);
509: return(0);
510: }
512: /* ------------------------------------------------------------------- */
515: /*
516: PreCheck - Optional user-defined routine that checks the validity of
517: candidate steps of a line search method. Set by SNESLineSearchSetPreCheck().
519: Input Parameters:
520: snes - the SNES context
521: xcurrent - current solution
522: y - search direction and length
524: Output Parameters:
525: y - proposed step (search direction and length) (possibly changed)
526: changed_y - tells if the step has changed or not
527: */
528: PetscErrorCode PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y, PetscBool *changed_y, void * ctx)
529: {
531: *changed_y = PETSC_FALSE;
532: return(0);
533: }
535: /* ------------------------------------------------------------------- */
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)
553:
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;
569: SNESLineSearchGetSNES(linesearch, &snes);
570: check = (StepCheckCtx *)ctx;
571: user = check->user;
572: SNESGetIterationNumber(snes,&iter);
573: SNESLineSearchGetPreCheck(linesearch, PETSC_NULL, (void **)&check);
574: /* iteration 1 indicates we are working on the second iteration */
575: if (iter > 0) {
576: da = user->da;
577: PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %G\n",iter,check->tolerance);
579: /* Access local array data */
580: DMDAVecGetArray(da,check->last_step,&xa_last);
581: DMDAVecGetArray(da,x,&xa);
582: DMDAGetCorners(da,&xs,PETSC_NULL,PETSC_NULL,&xm,PETSC_NULL,PETSC_NULL);
584: /*
585: If we fail the user-defined check for validity of the candidate iterate,
586: then modify the iterate as we like. (Note that the particular modification
587: below is intended simply to demonstrate how to manipulate this data, not
588: as a meaningful or appropriate choice.)
589: */
590: for (i=xs; i<xs+xm; i++) {
591: if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance; else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
592: if (rdiff > check->tolerance) {
593: tmp = xa[i];
594: xa[i] = .5*(xa[i] + xa_last[i]);
595: *changed_x = PETSC_TRUE;
596: PetscPrintf(PETSC_COMM_WORLD," Altering entry %D: x=%G, x_last=%G, diff=%G, x_new=%G\n",
597: i,PetscAbsScalar(tmp),PetscAbsScalar(xa_last[i]),rdiff,PetscAbsScalar(xa[i]));
598: }
599: }
600: DMDAVecRestoreArray(da,check->last_step,&xa_last);
601: DMDAVecRestoreArray(da,x,&xa);
602: }
603: VecCopy(x,check->last_step);
604: return(0);
605: }
608: /* ------------------------------------------------------------------- */
611: /*
612: PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
613: e.g,
614: 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
615: Set by SNESLineSearchSetPostCheck().
617: Input Parameters:
618: linesearch - the LineSearch context
619: xcurrent - current solution
620: y - search direction and length
621: x - the new candidate iterate
623: Output Parameters:
624: y - proposed step (search direction and length) (possibly changed)
625: x - current iterate (possibly modified)
626:
627: */
628: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool *changed_y,PetscBool *changed_x, void * ctx)
629: {
631: SetSubKSPCtx *check;
632: PetscInt iter,its,sub_its,maxit;
633: KSP ksp,sub_ksp,*sub_ksps;
634: PC pc;
635: PetscReal ksp_ratio;
636: SNES snes;
639: SNESLineSearchGetSNES(linesearch, &snes);
640: check = (SetSubKSPCtx *)ctx;
641: SNESGetIterationNumber(snes,&iter);
642: SNESGetKSP(snes,&ksp);
643: KSPGetPC(ksp,&pc);
644: PCBJacobiGetSubKSP(pc,PETSC_NULL,PETSC_NULL,&sub_ksps);
645: sub_ksp = sub_ksps[0];
646: KSPGetIterationNumber(ksp,&its); /* outer KSP iteration number */
647: KSPGetIterationNumber(sub_ksp,&sub_its); /* inner KSP iteration number */
649: if (iter){
650: PetscPrintf(PETSC_COMM_WORLD," ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);
651: ksp_ratio = ((PetscReal)(its))/check->its0;
652: maxit = (PetscInt)(ksp_ratio*sub_its + 0.5);
653: if (maxit < 2) maxit = 2;
654: KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
655: PetscPrintf(PETSC_COMM_WORLD," ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);
656: }
657: check->its0 = its; /* save current outer KSP iteration number */
658: return(0);
659: }