Actual source code: ex3.c

petsc-master 2016-09-28
Report Typos and Errors
  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);if (ierr) return ierr;
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);
125:   DMSetFromOptions(ctx.da);
126:   DMSetUp(ctx.da);

128:   /*
129:      Extract global and local vectors from DMDA; then duplicate for remaining
130:      vectors that are the same types
131:   */
132:   DMCreateGlobalVector(ctx.da,&x);
133:   VecDuplicate(x,&r);
134:   VecDuplicate(x,&F); ctx.F = F;
135:   VecDuplicate(x,&U);

137:   /*
138:      Set function evaluation routine and vector.  Whenever the nonlinear
139:      solver needs to compute the nonlinear function, it will call this
140:      routine.
141:       - Note that the final routine argument is the user-defined
142:         context that provides application-specific data for the
143:         function evaluation routine.
144:   */
145:   SNESSetFunction(snes,r,FormFunction,&ctx);

147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148:      Create matrix data structure; set Jacobian evaluation routine
149:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

151:   MatCreate(PETSC_COMM_WORLD,&J);
152:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
153:   MatSetFromOptions(J);
154:   MatSeqAIJSetPreallocation(J,3,NULL);
155:   MatMPIAIJSetPreallocation(J,3,NULL,3,NULL);

157:   /*
158:      Set Jacobian matrix data structure and default Jacobian evaluation
159:      routine.  Whenever the nonlinear solver needs to compute the
160:      Jacobian matrix, it will call this routine.
161:       - Note that the final routine argument is the user-defined
162:         context that provides application-specific data for the
163:         Jacobian evaluation routine.
164:   */
165:   SNESSetJacobian(snes,J,J,FormJacobian,&ctx);

167:   /*
168:      Optional allow user provided preconditioner
169:    */
170:   PetscOptionsHasName(NULL,NULL,"-user_precond",&flg);
171:   if (flg) {
172:     KSP ksp;
173:     PC  pc;
174:     SNESGetKSP(snes,&ksp);
175:     KSPGetPC(ksp,&pc);
176:     PCSetType(pc,PCSHELL);
177:     PCShellSetApply(pc,MatrixFreePreconditioner);
178:   }

180:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181:      Customize nonlinear solver; set runtime options
182:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

184:   /*
185:      Set an optional user-defined monitoring routine
186:   */
187:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
188:   SNESMonitorSet(snes,Monitor,&monP,0);

190:   /*
191:      Set names for some vectors to facilitate monitoring (optional)
192:   */
193:   PetscObjectSetName((PetscObject)x,"Approximate Solution");
194:   PetscObjectSetName((PetscObject)U,"Exact Solution");

196:   /*
197:      Set SNES/KSP/KSP/PC runtime options, e.g.,
198:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
199:   */
200:   SNESSetFromOptions(snes);

202:   /*
203:      Set an optional user-defined routine to check the validity of candidate
204:      iterates that are determined by line search methods
205:   */
206:   SNESGetLineSearch(snes, &linesearch);
207:   PetscOptionsHasName(NULL,NULL,"-post_check_iterates",&post_check);

209:   if (post_check) {
210:     PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
211:     SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP);
212:     VecDuplicate(x,&(checkP.last_step));

214:     checkP.tolerance = 1.0;
215:     checkP.user      = &ctx;

217:     PetscOptionsGetReal(NULL,NULL,"-check_tol",&checkP.tolerance,NULL);
218:   }

220:   PetscOptionsHasName(NULL,NULL,"-post_setsubksp",&post_setsubksp);
221:   if (post_setsubksp) {
222:     PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");
223:     SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1);
224:   }

226:   PetscOptionsHasName(NULL,NULL,"-pre_check_iterates",&pre_check);
227:   if (pre_check) {
228:     PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
229:     SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP);
230:   }

232:   /*
233:      Print parameters used for convergence testing (optional) ... just
234:      to demonstrate this routine; this information is also printed with
235:      the option -snes_view
236:   */
237:   SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);
238:   PetscPrintf(PETSC_COMM_WORLD,"atol=%g, rtol=%g, stol=%g, maxit=%D, maxf=%D\n",(double)abstol,(double)rtol,(double)stol,maxit,maxf);

240:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241:      Initialize application:
242:      Store right-hand-side of PDE and exact solution
243:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

245:   /*
246:      Get local grid boundaries (for 1-dimensional DMDA):
247:        xs, xm - starting grid index, width of local grid (no ghost points)
248:   */
249:   DMDAGetCorners(ctx.da,&xs,NULL,NULL,&xm,NULL,NULL);

251:   /*
252:      Get pointers to vector data
253:   */
254:   DMDAVecGetArray(ctx.da,F,&FF);
255:   DMDAVecGetArray(ctx.da,U,&UU);

257:   /*
258:      Compute local vector entries
259:   */
260:   xp = ctx.h*xs;
261:   for (i=xs; i<xs+xm; i++) {
262:     FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
263:     UU[i] = xp*xp*xp;
264:     xp   += ctx.h;
265:   }

267:   /*
268:      Restore vectors
269:   */
270:   DMDAVecRestoreArray(ctx.da,F,&FF);
271:   DMDAVecRestoreArray(ctx.da,U,&UU);

273:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274:      Evaluate initial guess; then solve nonlinear system
275:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

277:   /*
278:      Note: The user should initialize the vector, x, with the initial guess
279:      for the nonlinear solver prior to calling SNESSolve().  In particular,
280:      to employ an initial guess of zero, the user should explicitly set
281:      this vector to zero by calling VecSet().
282:   */
283:   FormInitialGuess(x);
284:   SNESSolve(snes,NULL,x);
285:   SNESGetIterationNumber(snes,&its);
286:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);

288:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289:      Check solution and clean up
290:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

292:   /*
293:      Check the error
294:   */
295:   VecAXPY(x,none,U);
296:   VecNorm(x,NORM_2,&norm);
297:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)norm,its);

299:   /*
300:      Free work space.  All PETSc objects should be destroyed when they
301:      are no longer needed.
302:   */
303:   PetscViewerDestroy(&monP.viewer);
304:   if (post_check) {VecDestroy(&checkP.last_step);}
305:   VecDestroy(&x);
306:   VecDestroy(&r);
307:   VecDestroy(&U);
308:   VecDestroy(&F);
309:   MatDestroy(&J);
310:   SNESDestroy(&snes);
311:   DMDestroy(&ctx.da);
312:   PetscFinalize();
313:   return ierr;
314: }
315: /* ------------------------------------------------------------------- */
318: /*
319:    FormInitialGuess - Computes initial guess.

321:    Input/Output Parameter:
322: .  x - the solution vector
323: */
324: PetscErrorCode FormInitialGuess(Vec x)
325: {
327:   PetscScalar    pfive = .50;

330:   VecSet(x,pfive);
331:   return(0);
332: }
333: /* ------------------------------------------------------------------- */
336: /*
337:    FormFunction - Evaluates nonlinear function, F(x).

339:    Input Parameters:
340: .  snes - the SNES context
341: .  x - input vector
342: .  ctx - optional user-defined context, as set by SNESSetFunction()

344:    Output Parameter:
345: .  f - function vector

347:    Note:
348:    The user-defined context can contain any application-specific
349:    data needed for the function evaluation.
350: */
351: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
352: {
353:   ApplicationCtx *user = (ApplicationCtx*) ctx;
354:   DM             da    = user->da;
355:   PetscScalar    *xx,*ff,*FF,d;
357:   PetscInt       i,M,xs,xm;
358:   Vec            xlocal;

361:   DMGetLocalVector(da,&xlocal);
362:   /*
363:      Scatter ghost points to local vector, using the 2-step process
364:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
365:      By placing code between these two statements, computations can
366:      be done while messages are in transition.
367:   */
368:   DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
369:   DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);

371:   /*
372:      Get pointers to vector data.
373:        - The vector xlocal includes ghost point; the vectors x and f do
374:          NOT include ghost points.
375:        - Using DMDAVecGetArray() allows accessing the values using global ordering
376:   */
377:   DMDAVecGetArray(da,xlocal,&xx);
378:   DMDAVecGetArray(da,f,&ff);
379:   DMDAVecGetArray(da,user->F,&FF);

381:   /*
382:      Get local grid boundaries (for 1-dimensional DMDA):
383:        xs, xm  - starting grid index, width of local grid (no ghost points)
384:   */
385:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
386:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

388:   /*
389:      Set function values for boundary points; define local interior grid point range:
390:         xsi - starting interior grid index
391:         xei - ending interior grid index
392:   */
393:   if (xs == 0) { /* left boundary */
394:     ff[0] = xx[0];
395:     xs++;xm--;
396:   }
397:   if (xs+xm == M) {  /* right boundary */
398:     ff[xs+xm-1] = xx[xs+xm-1] - 1.0;
399:     xm--;
400:   }

402:   /*
403:      Compute function over locally owned part of the grid (interior points only)
404:   */
405:   d = 1.0/(user->h*user->h);
406:   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];

408:   /*
409:      Restore vectors
410:   */
411:   DMDAVecRestoreArray(da,xlocal,&xx);
412:   DMDAVecRestoreArray(da,f,&ff);
413:   DMDAVecRestoreArray(da,user->F,&FF);
414:   DMRestoreLocalVector(da,&xlocal);
415:   return(0);
416: }
417: /* ------------------------------------------------------------------- */
420: /*
421:    FormJacobian - Evaluates Jacobian matrix.

423:    Input Parameters:
424: .  snes - the SNES context
425: .  x - input vector
426: .  dummy - optional user-defined context (not used here)

428:    Output Parameters:
429: .  jac - Jacobian matrix
430: .  B - optionally different preconditioning matrix
431: .  flag - flag indicating matrix structure
432: */
433: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
434: {
435:   ApplicationCtx *user = (ApplicationCtx*) ctx;
436:   PetscScalar    *xx,d,A[3];
438:   PetscInt       i,j[3],M,xs,xm;
439:   DM             da = user->da;

442:   /*
443:      Get pointer to vector data
444:   */
445:   DMDAVecGetArrayRead(da,x,&xx);
446:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

448:   /*
449:     Get range of locally owned matrix
450:   */
451:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

453:   /*
454:      Determine starting and ending local indices for interior grid points.
455:      Set Jacobian entries for boundary points.
456:   */

458:   if (xs == 0) {  /* left boundary */
459:     i = 0; A[0] = 1.0;

461:     MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
462:     xs++;xm--;
463:   }
464:   if (xs+xm == M) { /* right boundary */
465:     i    = M-1;
466:     A[0] = 1.0;
467:     MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
468:     xm--;
469:   }

471:   /*
472:      Interior grid points
473:       - Note that in this case we set all elements for a particular
474:         row at once.
475:   */
476:   d = 1.0/(user->h*user->h);
477:   for (i=xs; i<xs+xm; i++) {
478:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
479:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
480:     MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
481:   }

483:   /*
484:      Assemble matrix, using the 2-step process:
485:        MatAssemblyBegin(), MatAssemblyEnd().
486:      By placing code between these two statements, computations can be
487:      done while messages are in transition.

489:      Also, restore vector.
490:   */

492:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
493:   DMDAVecRestoreArrayRead(da,x,&xx);
494:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

496:   return(0);
497: }
498: /* ------------------------------------------------------------------- */
501: /*
502:    Monitor - Optional user-defined monitoring routine that views the
503:    current iterate with an x-window plot. Set by SNESMonitorSet().

505:    Input Parameters:
506:    snes - the SNES context
507:    its - iteration number
508:    norm - 2-norm function value (may be estimated)
509:    ctx - optional user-defined context for private data for the
510:          monitor routine, as set by SNESMonitorSet()

512:    Note:
513:    See the manpage for PetscViewerDrawOpen() for useful runtime options,
514:    such as -nox to deactivate all x-window output.
515:  */
516: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
517: {
519:   MonitorCtx     *monP = (MonitorCtx*) ctx;
520:   Vec            x;

523:   PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %g\n",its,(double)fnorm);
524:   SNESGetSolution(snes,&x);
525:   VecView(x,monP->viewer);
526:   return(0);
527: }

529: /* ------------------------------------------------------------------- */
532: /*
533:    PreCheck - Optional user-defined routine that checks the validity of
534:    candidate steps of a line search method.  Set by SNESLineSearchSetPreCheck().

536:    Input Parameters:
537:    snes - the SNES context
538:    xcurrent - current solution
539:    y - search direction and length

541:    Output Parameters:
542:    y         - proposed step (search direction and length) (possibly changed)
543:    changed_y - tells if the step has changed or not
544:  */
545: PetscErrorCode PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y, PetscBool *changed_y, void * ctx)
546: {
548:   *changed_y = PETSC_FALSE;
549:   return(0);
550: }

552: /* ------------------------------------------------------------------- */
555: /*
556:    PostCheck - Optional user-defined routine that checks the validity of
557:    candidate steps of a line search method.  Set by SNESLineSearchSetPostCheck().

559:    Input Parameters:
560:    snes - the SNES context
561:    ctx  - optional user-defined context for private data for the
562:           monitor routine, as set by SNESLineSearchSetPostCheck()
563:    xcurrent - current solution
564:    y - search direction and length
565:    x    - the new candidate iterate

567:    Output Parameters:
568:    y    - proposed step (search direction and length) (possibly changed)
569:    x    - current iterate (possibly modified)

571:  */
572: PetscErrorCode PostCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool  *changed_y,PetscBool  *changed_x, void * ctx)
573: {
575:   PetscInt       i,iter,xs,xm;
576:   StepCheckCtx   *check;
577:   ApplicationCtx *user;
578:   PetscScalar    *xa,*xa_last,tmp;
579:   PetscReal      rdiff;
580:   DM             da;
581:   SNES           snes;

584:   *changed_x = PETSC_FALSE;
585:   *changed_y = PETSC_FALSE;

587:   SNESLineSearchGetSNES(linesearch, &snes);
588:   check = (StepCheckCtx*)ctx;
589:   user  = check->user;
590:   SNESGetIterationNumber(snes,&iter);

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",i,(double)PetscAbsScalar(tmp),(double)PetscAbsScalar(xa_last[i]),(double)rdiff,(double)PetscAbsScalar(xa[i]));
616:       }
617:     }
618:     DMDAVecRestoreArray(da,check->last_step,&xa_last);
619:     DMDAVecRestoreArray(da,x,&xa);
620:   }
621:   VecCopy(x,check->last_step);
622:   return(0);
623: }


626: /* ------------------------------------------------------------------- */
629: /*
630:    PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
631:    e.g,
632:      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
633:    Set by SNESLineSearchSetPostCheck().

635:    Input Parameters:
636:    linesearch - the LineSearch context
637:    xcurrent - current solution
638:    y - search direction and length
639:    x    - the new candidate iterate

641:    Output Parameters:
642:    y    - proposed step (search direction and length) (possibly changed)
643:    x    - current iterate (possibly modified)

645:  */
646: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool  *changed_y,PetscBool  *changed_x, void * ctx)
647: {
649:   SetSubKSPCtx   *check;
650:   PetscInt       iter,its,sub_its,maxit;
651:   KSP            ksp,sub_ksp,*sub_ksps;
652:   PC             pc;
653:   PetscReal      ksp_ratio;
654:   SNES           snes;

657:   SNESLineSearchGetSNES(linesearch, &snes);
658:   check   = (SetSubKSPCtx*)ctx;
659:   SNESGetIterationNumber(snes,&iter);
660:   SNESGetKSP(snes,&ksp);
661:   KSPGetPC(ksp,&pc);
662:   PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);
663:   sub_ksp = sub_ksps[0];
664:   KSPGetIterationNumber(ksp,&its);      /* outer KSP iteration number */
665:   KSPGetIterationNumber(sub_ksp,&sub_its); /* inner KSP iteration number */

667:   if (iter) {
668:     PetscPrintf(PETSC_COMM_WORLD,"    ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);
669:     ksp_ratio = ((PetscReal)(its))/check->its0;
670:     maxit     = (PetscInt)(ksp_ratio*sub_its + 0.5);
671:     if (maxit < 2) maxit = 2;
672:     KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
673:     PetscPrintf(PETSC_COMM_WORLD,"    ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);
674:   }
675:   check->its0 = its; /* save current outer KSP iteration number */
676:   return(0);
677: }

679: /* ------------------------------------------------------------------- */
680: /*
681:    MatrixFreePreconditioner - This routine demonstrates the use of a
682:    user-provided preconditioner.  This code implements just the null
683:    preconditioner, which of course is not recommended for general use.

685:    Input Parameters:
686: +  pc - preconditioner
687: -  x - input vector

689:    Output Parameter:
690: .  y - preconditioned vector
691: */
692: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
693: {
695:   VecCopy(x,y);
696:   return 0;
697: }