Actual source code: ex3.c

petsc-3.6.1 2015-08-06
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:   -check_iterates : activate checking of iterates\n\
  8:   -check_tol <tol>: set tolerance for iterate checking\n\n";

 10: /*T
 11:    Concepts: SNES^basic parallel example
 12:    Concepts: SNES^setting a user-defined monitoring routine
 13:    Concepts: error handling^using the macro __FUNCT__ to define routine names;
 14:    Processors: n
 15: T*/

 17: /*
 18:    Include "petscdraw.h" so that we can use distributed arrays (DMDAs).
 19:    Include "petscdraw.h" so that we can use PETSc drawing routines.
 20:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 21:    file automatically includes:
 22:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 23:      petscmat.h - matrices
 24:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 25:      petscviewer.h - viewers               petscpc.h  - preconditioners
 26:      petscksp.h   - linear solvers
 27: */

 29: #include <petscdm.h>
 30: #include <petscdmda.h>
 31: #include <petscsnes.h>

 33: /*
 34:    User-defined routines.  Note that immediately before each routine below,
 35:    we define the macro __FUNCT__ to be a string containing the routine name.
 36:    If defined, this macro is used in the PETSc error handlers to provide a
 37:    complete traceback of routine names.  All PETSc library routines use this
 38:    macro, and users can optionally employ it as well in their application
 39:    codes.  Note that users can get a traceback of PETSc errors regardless of
 40:    whether they define __FUNCT__ in application codes; this macro merely
 41:    provides the added traceback detail of the application routine names.
 42: */
 43: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 44: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 45: PetscErrorCode FormInitialGuess(Vec);
 46: PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
 47: PetscErrorCode PreCheck(SNESLineSearch,Vec,Vec,PetscBool*,void*);
 48: PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
 49: PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
 50: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);

 52: /*
 53:    User-defined application context
 54: */
 55: typedef struct {
 56:   DM          da;      /* distributed array */
 57:   Vec         F;       /* right-hand-side of PDE */
 58:   PetscMPIInt rank;    /* rank of processor */
 59:   PetscMPIInt size;    /* size of communicator */
 60:   PetscReal   h;       /* mesh spacing */
 61: } ApplicationCtx;

 63: /*
 64:    User-defined context for monitoring
 65: */
 66: typedef struct {
 67:   PetscViewer viewer;
 68: } MonitorCtx;

 70: /*
 71:    User-defined context for checking candidate iterates that are
 72:    determined by line search methods
 73: */
 74: typedef struct {
 75:   Vec            last_step;  /* previous iterate */
 76:   PetscReal      tolerance;  /* tolerance for changes between successive iterates */
 77:   ApplicationCtx *user;
 78: } StepCheckCtx;

 80: typedef struct {
 81:   PetscInt its0; /* num of prevous outer KSP iterations */
 82: } SetSubKSPCtx;

 86: int main(int argc,char **argv)
 87: {
 88:   SNES           snes;                 /* SNES context */
 89:   SNESLineSearch linesearch;          /* SNESLineSearch context */
 90:   Mat            J;                    /* Jacobian matrix */
 91:   ApplicationCtx ctx;                  /* user-defined context */
 92:   Vec            x,r,U,F;              /* vectors */
 93:   MonitorCtx     monP;                 /* monitoring context */
 94:   StepCheckCtx   checkP;               /* step-checking context */
 95:   SetSubKSPCtx   checkP1;
 96:   PetscBool      pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking candidate iterates */
 97:   PetscScalar    xp,*FF,*UU,none = -1.0;
 99:   PetscInt       its,N = 5,i,maxit,maxf,xs,xm;
100:   PetscReal      abstol,rtol,stol,norm;
101:   PetscBool      flg;


104:   PetscInitialize(&argc,&argv,(char*)0,help);
105:   MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
106:   MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
107:   PetscOptionsGetInt(NULL,"-n",&N,NULL);
108:   ctx.h = 1.0/(N-1);

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:      Create nonlinear solver context
112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

114:   SNESCreate(PETSC_COMM_WORLD,&snes);

116:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117:      Create vector data structures; set function evaluation routine
118:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

120:   /*
121:      Create distributed array (DMDA) to manage parallel grid and vectors
122:   */
123:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,1,1,NULL,&ctx.da);

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

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

144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145:      Create matrix data structure; set Jacobian evaluation routine
146:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

148:   MatCreate(PETSC_COMM_WORLD,&J);
149:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
150:   MatSetFromOptions(J);
151:   MatSeqAIJSetPreallocation(J,3,NULL);
152:   MatMPIAIJSetPreallocation(J,3,NULL,3,NULL);

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

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

177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:      Customize nonlinear solver; set runtime options
179:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

211:     checkP.tolerance = 1.0;
212:     checkP.user      = &ctx;

214:     PetscOptionsGetReal(NULL,"-check_tol",&checkP.tolerance,NULL);
215:   }

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

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

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

237:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238:      Initialize application:
239:      Store right-hand-side of PDE and exact solution
240:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

248:   /*
249:      Get pointers to vector data
250:   */
251:   DMDAVecGetArray(ctx.da,F,&FF);
252:   DMDAVecGetArray(ctx.da,U,&UU);

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

264:   /*
265:      Restore vectors
266:   */
267:   DMDAVecRestoreArray(ctx.da,F,&FF);
268:   DMDAVecRestoreArray(ctx.da,U,&UU);

270:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
271:      Evaluate initial guess; then solve nonlinear system
272:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

285:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286:      Check solution and clean up
287:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

327:   VecSet(x,pfive);
328:   return(0);
329: }
330: /* ------------------------------------------------------------------- */
333: /*
334:    FormFunction - Evaluates nonlinear function, F(x).

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

341:    Output Parameter:
342: .  f - function vector

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

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

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

378:   /*
379:      Get local grid boundaries (for 1-dimensional DMDA):
380:        xs, xm  - starting grid index, width of local grid (no ghost points)
381:   */
382:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
383:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,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,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:   DMDAVecGetArrayRead(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,NULL,NULL,NULL,NULL,NULL);

450:   /*
451:      Determine starting and ending local indices for interior grid points.
452:      Set Jacobian entries for boundary points.
453:   */

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

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

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

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

486:      Also, restore vector.
487:   */

489:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
490:   DMDAVecRestoreArrayRead(da,x,&xx);
491:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

493:   return(0);
494: }
495: /* ------------------------------------------------------------------- */
498: /*
499:    Monitor - Optional user-defined monitoring routine that views the
500:    current iterate with an x-window plot. Set by SNESMonitorSet().

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

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

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

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

533:    Input Parameters:
534:    snes - the SNES context
535:    xcurrent - current solution
536:    y - search direction and length

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

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

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

564:    Output Parameters:
565:    y    - proposed step (search direction and length) (possibly changed)
566:    x    - current iterate (possibly modified)

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

581:   *changed_x = PETSC_FALSE;
582:   *changed_y = PETSC_FALSE;

584:   SNESLineSearchGetSNES(linesearch, &snes);
585:   check = (StepCheckCtx*)ctx;
586:   user  = check->user;
587:   SNESGetIterationNumber(snes,&iter);
588:   SNESLineSearchGetPreCheck(linesearch, NULL, (void**)&check);

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: }