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