Actual source code: ex3.c

petsc-master 2018-05-22
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*/



 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
 24:      petscvec.h    - vectors
 25:      petscmat.h    - matrices
 26:      petscis.h     - index sets
 27:      petscksp.h    - Krylov subspace methods
 28:      petscviewer.h - viewers
 29:      petscpc.h     - preconditioners
 30:      petscksp.h    - linear solvers
 31: */

 33:  #include <petscdm.h>
 34:  #include <petscdmda.h>
 35:  #include <petscsnes.h>

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

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

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

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

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

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

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

104:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105:      Create nonlinear solver context
106:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

108:   SNESCreate(PETSC_COMM_WORLD,&snes);

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:      Create vector data structures; set function evaluation routine
112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141:      Create matrix data structure; set Jacobian evaluation routine
142:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

173:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
174:      Customize nonlinear solver; set runtime options
175:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

207:     checkP.tolerance = 1.0;
208:     checkP.user      = &ctx;

210:     PetscOptionsGetReal(NULL,NULL,"-check_tol",&checkP.tolerance,NULL);
211:   }

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

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

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

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

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

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

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

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

266:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267:      Evaluate initial guess; then solve nonlinear system
268:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

281:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
282:      Check solution and clean up
283:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

309: /* ------------------------------------------------------------------- */
310: /*
311:    FormInitialGuess - Computes initial guess.

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

322:   VecSet(x,pfive);
323:   return(0);
324: }

326: /* ------------------------------------------------------------------- */
327: /*
328:    FormFunction - Evaluates nonlinear function, F(x).

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

335:    Output Parameter:
336: .  f - function vector

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

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

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

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

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

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

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

409: /* ------------------------------------------------------------------- */
410: /*
411:    FormJacobian - Evaluates Jacobian matrix.

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

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

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

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

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

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

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

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

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

479:      Also, restore vector.
480:   */

482:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
483:   DMDAVecRestoreArrayRead(da,x,&xx);
484:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

486:   return(0);
487: }

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

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

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

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

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

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

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

539: /* ------------------------------------------------------------------- */
540: /*
541:    PostCheck - Optional user-defined routine that checks the validity of
542:    candidate steps of a line search method.  Set by SNESLineSearchSetPostCheck().

544:    Input Parameters:
545:    snes - the SNES context
546:    ctx  - optional user-defined context for private data for the
547:           monitor routine, as set by SNESLineSearchSetPostCheck()
548:    xcurrent - current solution
549:    y - search direction and length
550:    x    - the new candidate iterate

552:    Output Parameters:
553:    y    - proposed step (search direction and length) (possibly changed)
554:    x    - current iterate (possibly modified)

556:  */
557: PetscErrorCode PostCheck(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool  *changed_y,PetscBool  *changed_x, void * ctx)
558: {
560:   PetscInt       i,iter,xs,xm;
561:   StepCheckCtx   *check;
562:   ApplicationCtx *user;
563:   PetscScalar    *xa,*xa_last,tmp;
564:   PetscReal      rdiff;
565:   DM             da;
566:   SNES           snes;

569:   *changed_x = PETSC_FALSE;
570:   *changed_y = PETSC_FALSE;

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

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

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

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

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

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

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

667:    Input Parameters:
668: +  pc - preconditioner
669: -  x - input vector

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


682: /*TEST

684:    test:
685:       args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

687:    test:
688:       suffix: 2
689:       nsize: 3
690:       args: -nox -pc_type asm -mat_type mpiaij -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

692:    test:
693:       suffix: 3
694:       nsize: 2
695:       args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

697:    test:
698:       suffix: 4
699:       args: -nox -pre_check_iterates -post_check_iterates

701: TEST*/