Actual source code: ex3.c

petsc-master 2017-11-20
Report Typos and Errors
  1: static char help[] = "Newton methods to solve u'' + u^{2} = f in parallel.\n\
  2: This example employs a user-defined monitoring routine and optionally a user-defined\n\
  3: routine to check candidate iterates produced by line search routines.\n\
  4: The command line options include:\n\
  5:   -pre_check_iterates : activate checking of iterates\n\
  6:   -post_check_iterates : activate checking of iterates\n\
  7:   -check_tol <tol>: set tolerance for iterate checking\n\
  8:   -user_precond : activate a (trivial) user-defined preconditioner\n\n";

 10: /*T
 11:    Concepts: SNES^basic parallel example
 12:    Concepts: SNES^setting a user-defined monitoring routine
 13:    Processors: n
 14: T*/

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

 31:  #include <petscdm.h>
 32:  #include <petscdmda.h>
 33:  #include <petscsnes.h>

 35: /*
 36:    User-defined routines.
 37: */
 38: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 39: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 40: PetscErrorCode FormInitialGuess(Vec);
 41: PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
 42: PetscErrorCode PreCheck(SNESLineSearch,Vec,Vec,PetscBool*,void*);
 43: PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
 44: PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
 45: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);

 47: /*
 48:    User-defined application context
 49: */
 50: typedef struct {
 51:   DM          da;      /* distributed array */
 52:   Vec         F;       /* right-hand-side of PDE */
 53:   PetscMPIInt rank;    /* rank of processor */
 54:   PetscMPIInt size;    /* size of communicator */
 55:   PetscReal   h;       /* mesh spacing */
 56: } ApplicationCtx;

 58: /*
 59:    User-defined context for monitoring
 60: */
 61: typedef struct {
 62:   PetscViewer viewer;
 63: } MonitorCtx;

 65: /*
 66:    User-defined context for checking candidate iterates that are
 67:    determined by line search methods
 68: */
 69: typedef struct {
 70:   Vec            last_step;  /* previous iterate */
 71:   PetscReal      tolerance;  /* tolerance for changes between successive iterates */
 72:   ApplicationCtx *user;
 73: } StepCheckCtx;

 75: typedef struct {
 76:   PetscInt its0; /* num of prevous outer KSP iterations */
 77: } SetSubKSPCtx;

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

 96:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 97:   MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
 98:   MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
 99:   PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
100:   ctx.h = 1.0/(N-1);

102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103:      Create nonlinear solver context
104:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

106:   SNESCreate(PETSC_COMM_WORLD,&snes);

108:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109:      Create vector data structures; set function evaluation routine
110:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

112:   /*
113:      Create distributed array (DMDA) to manage parallel grid and vectors
114:   */
115:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,1,1,NULL,&ctx.da);
116:   DMSetFromOptions(ctx.da);
117:   DMSetUp(ctx.da);

119:   /*
120:      Extract global and local vectors from DMDA; then duplicate for remaining
121:      vectors that are the same types
122:   */
123:   DMCreateGlobalVector(ctx.da,&x);
124:   VecDuplicate(x,&r);
125:   VecDuplicate(x,&F); ctx.F = F;
126:   VecDuplicate(x,&U);

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

138:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139:      Create matrix data structure; set Jacobian evaluation routine
140:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

142:   MatCreate(PETSC_COMM_WORLD,&J);
143:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
144:   MatSetFromOptions(J);
145:   MatSeqAIJSetPreallocation(J,3,NULL);
146:   MatMPIAIJSetPreallocation(J,3,NULL,3,NULL);

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

158:   /*
159:      Optionally allow user-provided preconditioner
160:    */
161:   PetscOptionsHasName(NULL,NULL,"-user_precond",&flg);
162:   if (flg) {
163:     KSP ksp;
164:     PC  pc;
165:     SNESGetKSP(snes,&ksp);
166:     KSPGetPC(ksp,&pc);
167:     PCSetType(pc,PCSHELL);
168:     PCShellSetApply(pc,MatrixFreePreconditioner);
169:   }

171:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172:      Customize nonlinear solver; set runtime options
173:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

175:   /*
176:      Set an optional user-defined monitoring routine
177:   */
178:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);
179:   SNESMonitorSet(snes,Monitor,&monP,0);

181:   /*
182:      Set names for some vectors to facilitate monitoring (optional)
183:   */
184:   PetscObjectSetName((PetscObject)x,"Approximate Solution");
185:   PetscObjectSetName((PetscObject)U,"Exact Solution");

187:   /*
188:      Set SNES/KSP/KSP/PC runtime options, e.g.,
189:          -snes_view -snes_monitor -ksp_type <ksp> -pc_type <pc>
190:   */
191:   SNESSetFromOptions(snes);

193:   /*
194:      Set an optional user-defined routine to check the validity of candidate
195:      iterates that are determined by line search methods
196:   */
197:   SNESGetLineSearch(snes, &linesearch);
198:   PetscOptionsHasName(NULL,NULL,"-post_check_iterates",&post_check);

200:   if (post_check) {
201:     PetscPrintf(PETSC_COMM_WORLD,"Activating post step checking routine\n");
202:     SNESLineSearchSetPostCheck(linesearch,PostCheck,&checkP);
203:     VecDuplicate(x,&(checkP.last_step));

205:     checkP.tolerance = 1.0;
206:     checkP.user      = &ctx;

208:     PetscOptionsGetReal(NULL,NULL,"-check_tol",&checkP.tolerance,NULL);
209:   }

211:   PetscOptionsHasName(NULL,NULL,"-post_setsubksp",&post_setsubksp);
212:   if (post_setsubksp) {
213:     PetscPrintf(PETSC_COMM_WORLD,"Activating post setsubksp\n");
214:     SNESLineSearchSetPostCheck(linesearch,PostSetSubKSP,&checkP1);
215:   }

217:   PetscOptionsHasName(NULL,NULL,"-pre_check_iterates",&pre_check);
218:   if (pre_check) {
219:     PetscPrintf(PETSC_COMM_WORLD,"Activating pre step checking routine\n");
220:     SNESLineSearchSetPreCheck(linesearch,PreCheck,&checkP);
221:   }

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

231:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232:      Initialize application:
233:      Store right-hand-side of PDE and exact solution
234:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

236:   /*
237:      Get local grid boundaries (for 1-dimensional DMDA):
238:        xs, xm - starting grid index, width of local grid (no ghost points)
239:   */
240:   DMDAGetCorners(ctx.da,&xs,NULL,NULL,&xm,NULL,NULL);

242:   /*
243:      Get pointers to vector data
244:   */
245:   DMDAVecGetArray(ctx.da,F,&FF);
246:   DMDAVecGetArray(ctx.da,U,&UU);

248:   /*
249:      Compute local vector entries
250:   */
251:   xp = ctx.h*xs;
252:   for (i=xs; i<xs+xm; i++) {
253:     FF[i] = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 is to prevent 0^6 */
254:     UU[i] = xp*xp*xp;
255:     xp   += ctx.h;
256:   }

258:   /*
259:      Restore vectors
260:   */
261:   DMDAVecRestoreArray(ctx.da,F,&FF);
262:   DMDAVecRestoreArray(ctx.da,U,&UU);

264:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
265:      Evaluate initial guess; then solve nonlinear system
266:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

279:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
280:      Check solution and clean up
281:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

283:   /*
284:      Check the error
285:   */
286:   VecAXPY(x,none,U);
287:   VecNorm(x,NORM_2,&norm);
288:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)norm,its);

290:   /*
291:      Free work space.  All PETSc objects should be destroyed when they
292:      are no longer needed.
293:   */
294:   PetscViewerDestroy(&monP.viewer);
295:   if (post_check) {VecDestroy(&checkP.last_step);}
296:   VecDestroy(&x);
297:   VecDestroy(&r);
298:   VecDestroy(&U);
299:   VecDestroy(&F);
300:   MatDestroy(&J);
301:   SNESDestroy(&snes);
302:   DMDestroy(&ctx.da);
303:   PetscFinalize();
304:   return ierr;
305: }

307: /* ------------------------------------------------------------------- */
308: /*
309:    FormInitialGuess - Computes initial guess.

311:    Input/Output Parameter:
312: .  x - the solution vector
313: */
314: PetscErrorCode FormInitialGuess(Vec x)
315: {
317:   PetscScalar    pfive = .50;

320:   VecSet(x,pfive);
321:   return(0);
322: }

324: /* ------------------------------------------------------------------- */
325: /*
326:    FormFunction - Evaluates nonlinear function, F(x).

328:    Input Parameters:
329: .  snes - the SNES context
330: .  x - input vector
331: .  ctx - optional user-defined context, as set by SNESSetFunction()

333:    Output Parameter:
334: .  f - function vector

336:    Note:
337:    The user-defined context can contain any application-specific
338:    data needed for the function evaluation.
339: */
340: PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
341: {
342:   ApplicationCtx *user = (ApplicationCtx*) ctx;
343:   DM             da    = user->da;
344:   PetscScalar    *xx,*ff,*FF,d;
346:   PetscInt       i,M,xs,xm;
347:   Vec            xlocal;

350:   DMGetLocalVector(da,&xlocal);
351:   /*
352:      Scatter ghost points to local vector, using the 2-step process
353:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
354:      By placing code between these two statements, computations can
355:      be done while messages are in transition.
356:   */
357:   DMGlobalToLocalBegin(da,x,INSERT_VALUES,xlocal);
358:   DMGlobalToLocalEnd(da,x,INSERT_VALUES,xlocal);

360:   /*
361:      Get pointers to vector data.
362:        - The vector xlocal includes ghost point; the vectors x and f do
363:          NOT include ghost points.
364:        - Using DMDAVecGetArray() allows accessing the values using global ordering
365:   */
366:   DMDAVecGetArray(da,xlocal,&xx);
367:   DMDAVecGetArray(da,f,&ff);
368:   DMDAVecGetArray(da,user->F,&FF);

370:   /*
371:      Get local grid boundaries (for 1-dimensional DMDA):
372:        xs, xm  - starting grid index, width of local grid (no ghost points)
373:   */
374:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
375:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

377:   /*
378:      Set function values for boundary points; define local interior grid point range:
379:         xsi - starting interior grid index
380:         xei - ending interior grid index
381:   */
382:   if (xs == 0) { /* left boundary */
383:     ff[0] = xx[0];
384:     xs++;xm--;
385:   }
386:   if (xs+xm == M) {  /* right boundary */
387:     ff[xs+xm-1] = xx[xs+xm-1] - 1.0;
388:     xm--;
389:   }

391:   /*
392:      Compute function over locally owned part of the grid (interior points only)
393:   */
394:   d = 1.0/(user->h*user->h);
395:   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];

397:   /*
398:      Restore vectors
399:   */
400:   DMDAVecRestoreArray(da,xlocal,&xx);
401:   DMDAVecRestoreArray(da,f,&ff);
402:   DMDAVecRestoreArray(da,user->F,&FF);
403:   DMRestoreLocalVector(da,&xlocal);
404:   return(0);
405: }

407: /* ------------------------------------------------------------------- */
408: /*
409:    FormJacobian - Evaluates Jacobian matrix.

411:    Input Parameters:
412: .  snes - the SNES context
413: .  x - input vector
414: .  dummy - optional user-defined context (not used here)

416:    Output Parameters:
417: .  jac - Jacobian matrix
418: .  B - optionally different preconditioning matrix
419: .  flag - flag indicating matrix structure
420: */
421: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
422: {
423:   ApplicationCtx *user = (ApplicationCtx*) ctx;
424:   PetscScalar    *xx,d,A[3];
426:   PetscInt       i,j[3],M,xs,xm;
427:   DM             da = user->da;

430:   /*
431:      Get pointer to vector data
432:   */
433:   DMDAVecGetArrayRead(da,x,&xx);
434:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

436:   /*
437:     Get range of locally owned matrix
438:   */
439:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

441:   /*
442:      Determine starting and ending local indices for interior grid points.
443:      Set Jacobian entries for boundary points.
444:   */

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

449:     MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
450:     xs++;xm--;
451:   }
452:   if (xs+xm == M) { /* right boundary */
453:     i    = M-1;
454:     A[0] = 1.0;
455:     MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
456:     xm--;
457:   }

459:   /*
460:      Interior grid points
461:       - Note that in this case we set all elements for a particular
462:         row at once.
463:   */
464:   d = 1.0/(user->h*user->h);
465:   for (i=xs; i<xs+xm; i++) {
466:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
467:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
468:     MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
469:   }

471:   /*
472:      Assemble matrix, using the 2-step process:
473:        MatAssemblyBegin(), MatAssemblyEnd().
474:      By placing code between these two statements, computations can be
475:      done while messages are in transition.

477:      Also, restore vector.
478:   */

480:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
481:   DMDAVecRestoreArrayRead(da,x,&xx);
482:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

484:   return(0);
485: }

487: /* ------------------------------------------------------------------- */
488: /*
489:    Monitor - Optional user-defined monitoring routine that views the
490:    current iterate with an x-window plot. Set by SNESMonitorSet().

492:    Input Parameters:
493:    snes - the SNES context
494:    its - iteration number
495:    norm - 2-norm function value (may be estimated)
496:    ctx - optional user-defined context for private data for the
497:          monitor routine, as set by SNESMonitorSet()

499:    Note:
500:    See the manpage for PetscViewerDrawOpen() for useful runtime options,
501:    such as -nox to deactivate all x-window output.
502:  */
503: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
504: {
506:   MonitorCtx     *monP = (MonitorCtx*) ctx;
507:   Vec            x;

510:   PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %g\n",its,(double)fnorm);
511:   SNESGetSolution(snes,&x);
512:   VecView(x,monP->viewer);
513:   return(0);
514: }

516: /* ------------------------------------------------------------------- */
517: /*
518:    PreCheck - Optional user-defined routine that checks the validity of
519:    candidate steps of a line search method.  Set by SNESLineSearchSetPreCheck().

521:    Input Parameters:
522:    snes - the SNES context
523:    xcurrent - current solution
524:    y - search direction and length

526:    Output Parameters:
527:    y         - proposed step (search direction and length) (possibly changed)
528:    changed_y - tells if the step has changed or not
529:  */
530: PetscErrorCode PreCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y, PetscBool *changed_y, void * ctx)
531: {
533:   *changed_y = PETSC_FALSE;
534:   return(0);
535: }

537: /* ------------------------------------------------------------------- */
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)

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;

570:   SNESLineSearchGetSNES(linesearch, &snes);
571:   check = (StepCheckCtx*)ctx;
572:   user  = check->user;
573:   SNESGetIterationNumber(snes,&iter);

575:   /* iteration 1 indicates we are working on the second iteration */
576:   if (iter > 0) {
577:     da   = user->da;
578:     PetscPrintf(PETSC_COMM_WORLD,"Checking candidate step at iteration %D with tolerance %g\n",iter,(double)check->tolerance);

580:     /* Access local array data */
581:     DMDAVecGetArray(da,check->last_step,&xa_last);
582:     DMDAVecGetArray(da,x,&xa);
583:     DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

585:     /*
586:        If we fail the user-defined check for validity of the candidate iterate,
587:        then modify the iterate as we like.  (Note that the particular modification
588:        below is intended simply to demonstrate how to manipulate this data, not
589:        as a meaningful or appropriate choice.)
590:     */
591:     for (i=xs; i<xs+xm; i++) {
592:       if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance;
593:       else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
594:       if (rdiff > check->tolerance) {
595:         tmp        = xa[i];
596:         xa[i]      = .5*(xa[i] + xa_last[i]);
597:         *changed_x = PETSC_TRUE;
598:         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]));
599:       }
600:     }
601:     DMDAVecRestoreArray(da,check->last_step,&xa_last);
602:     DMDAVecRestoreArray(da,x,&xa);
603:   }
604:   VecCopy(x,check->last_step);
605:   return(0);
606: }

608: /* ------------------------------------------------------------------- */
609: /*
610:    PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
611:    e.g,
612:      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
613:    Set by SNESLineSearchSetPostCheck().

615:    Input Parameters:
616:    linesearch - the LineSearch context
617:    xcurrent - current solution
618:    y - search direction and length
619:    x    - the new candidate iterate

621:    Output Parameters:
622:    y    - proposed step (search direction and length) (possibly changed)
623:    x    - current iterate (possibly modified)

625:  */
626: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool  *changed_y,PetscBool  *changed_x, void * ctx)
627: {
629:   SetSubKSPCtx   *check;
630:   PetscInt       iter,its,sub_its,maxit;
631:   KSP            ksp,sub_ksp,*sub_ksps;
632:   PC             pc;
633:   PetscReal      ksp_ratio;
634:   SNES           snes;

637:   SNESLineSearchGetSNES(linesearch, &snes);
638:   check   = (SetSubKSPCtx*)ctx;
639:   SNESGetIterationNumber(snes,&iter);
640:   SNESGetKSP(snes,&ksp);
641:   KSPGetPC(ksp,&pc);
642:   PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);
643:   sub_ksp = sub_ksps[0];
644:   KSPGetIterationNumber(ksp,&its);      /* outer KSP iteration number */
645:   KSPGetIterationNumber(sub_ksp,&sub_its); /* inner KSP iteration number */

647:   if (iter) {
648:     PetscPrintf(PETSC_COMM_WORLD,"    ...PostCheck snes iteration %D, ksp_it %d %d, subksp_it %d\n",iter,check->its0,its,sub_its);
649:     ksp_ratio = ((PetscReal)(its))/check->its0;
650:     maxit     = (PetscInt)(ksp_ratio*sub_its + 0.5);
651:     if (maxit < 2) maxit = 2;
652:     KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
653:     PetscPrintf(PETSC_COMM_WORLD,"    ...ksp_ratio %g, new maxit %d\n\n",ksp_ratio,maxit);
654:   }
655:   check->its0 = its; /* save current outer KSP iteration number */
656:   return(0);
657: }

659: /* ------------------------------------------------------------------- */
660: /*
661:    MatrixFreePreconditioner - This routine demonstrates the use of a
662:    user-provided preconditioner.  This code implements just the null
663:    preconditioner, which of course is not recommended for general use.

665:    Input Parameters:
666: +  pc - preconditioner
667: -  x - input vector

669:    Output Parameter:
670: .  y - preconditioned vector
671: */
672: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
673: {
675:   VecCopy(x,y);
676:   return 0;
677: }