Actual source code: ex3.c

petsc-3.14.6 2021-03-30
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 "petscdm.h" so that we can use data management objects (DMs)
 20:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 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:   PetscBool   sjerr;   /* if or not to test jacobian domain error */
 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;

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


100:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
101:   MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
102:   MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
103:   PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
104:   ctx.h = 1.0/(N-1);
105:   ctx.sjerr = PETSC_FALSE;
106:   PetscOptionsGetBool(NULL,NULL,"-test_jacobian_domain_error",&ctx.sjerr,NULL);
107:   PetscOptionsGetBool(NULL,NULL,"-view_initial",&viewinitial,NULL);

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

113:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

265:   /*
266:      Restore vectors
267:   */
268:   DMDAVecRestoreArray(ctx.da,F,&FF);
269:   DMDAVecRestoreArray(ctx.da,U,&UU);
270:   if (viewinitial) {
271:     VecView(U,0);
272:     VecView(F,0);
273:   }

275:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276:      Evaluate initial guess; then solve nonlinear system
277:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

290:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
291:      Check solution and clean up
292:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

294:   /*
295:      Check the error
296:   */
297:   VecAXPY(x,none,U);
298:   VecNorm(x,NORM_2,&norm);
299:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)norm,its);
300:   if (ctx.sjerr) {
301:     SNESType       snestype;
302:     SNESGetType(snes,&snestype);
303:     PetscPrintf(PETSC_COMM_WORLD,"SNES Type %s\n",snestype);
304:   }

306:   /*
307:      Free work space.  All PETSc objects should be destroyed when they
308:      are no longer needed.
309:   */
310:   PetscViewerDestroy(&monP.viewer);
311:   if (post_check) {VecDestroy(&checkP.last_step);}
312:   VecDestroy(&x);
313:   VecDestroy(&r);
314:   VecDestroy(&U);
315:   VecDestroy(&F);
316:   MatDestroy(&J);
317:   SNESDestroy(&snes);
318:   DMDestroy(&ctx.da);
319:   PetscFinalize();
320:   return ierr;
321: }

323: /* ------------------------------------------------------------------- */
324: /*
325:    FormInitialGuess - Computes initial guess.

327:    Input/Output Parameter:
328: .  x - the solution vector
329: */
330: PetscErrorCode FormInitialGuess(Vec x)
331: {
333:   PetscScalar    pfive = .50;

336:   VecSet(x,pfive);
337:   return(0);
338: }

340: /* ------------------------------------------------------------------- */
341: /*
342:    FormFunction - Evaluates nonlinear function, F(x).

344:    Input Parameters:
345: .  snes - the SNES context
346: .  x - input vector
347: .  ctx - optional user-defined context, as set by SNESSetFunction()

349:    Output Parameter:
350: .  f - function vector

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

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

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

386:   /*
387:      Get local grid boundaries (for 1-dimensional DMDA):
388:        xs, xm  - starting grid index, width of local grid (no ghost points)
389:   */
390:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
391:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

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

407:   /*
408:      Compute function over locally owned part of the grid (interior points only)
409:   */
410:   d = 1.0/(user->h*user->h);
411:   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];

413:   /*
414:      Restore vectors
415:   */
416:   DMDAVecRestoreArray(da,xlocal,&xx);
417:   DMDAVecRestoreArray(da,f,&ff);
418:   DMDAVecRestoreArray(da,user->F,&FF);
419:   DMRestoreLocalVector(da,&xlocal);
420:   return(0);
421: }

423: /* ------------------------------------------------------------------- */
424: /*
425:    FormJacobian - Evaluates Jacobian matrix.

427:    Input Parameters:
428: .  snes - the SNES context
429: .  x - input vector
430: .  dummy - optional user-defined context (not used here)

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

446:   if (user->sjerr) {
447:     SNESSetJacobianDomainError(snes);
448:     return(0);
449:   }
450:   /*
451:      Get pointer to vector data
452:   */
453:   DMDAVecGetArrayRead(da,x,&xx);
454:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

456:   /*
457:     Get range of locally owned matrix
458:   */
459:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

461:   /*
462:      Determine starting and ending local indices for interior grid points.
463:      Set Jacobian entries for boundary points.
464:   */

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

469:     MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
470:     xs++;xm--;
471:   }
472:   if (xs+xm == M) { /* right boundary */
473:     i    = M-1;
474:     A[0] = 1.0;
475:     MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
476:     xm--;
477:   }

479:   /*
480:      Interior grid points
481:       - Note that in this case we set all elements for a particular
482:         row at once.
483:   */
484:   d = 1.0/(user->h*user->h);
485:   for (i=xs; i<xs+xm; i++) {
486:     j[0] = i - 1; j[1] = i; j[2] = i + 1;
487:     A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
488:     MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
489:   }

491:   /*
492:      Assemble matrix, using the 2-step process:
493:        MatAssemblyBegin(), MatAssemblyEnd().
494:      By placing code between these two statements, computations can be
495:      done while messages are in transition.

497:      Also, restore vector.
498:   */

500:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
501:   DMDAVecRestoreArrayRead(da,x,&xx);
502:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

504:   return(0);
505: }

507: /* ------------------------------------------------------------------- */
508: /*
509:    Monitor - Optional user-defined monitoring routine that views the
510:    current iterate with an x-window plot. Set by SNESMonitorSet().

512:    Input Parameters:
513:    snes - the SNES context
514:    its - iteration number
515:    norm - 2-norm function value (may be estimated)
516:    ctx - optional user-defined context for private data for the
517:          monitor routine, as set by SNESMonitorSet()

519:    Note:
520:    See the manpage for PetscViewerDrawOpen() for useful runtime options,
521:    such as -nox to deactivate all x-window output.
522:  */
523: PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
524: {
526:   MonitorCtx     *monP = (MonitorCtx*) ctx;
527:   Vec            x;

530:   PetscPrintf(PETSC_COMM_WORLD,"iter = %D,SNES Function norm %g\n",its,(double)fnorm);
531:   SNESGetSolution(snes,&x);
532:   VecView(x,monP->viewer);
533:   return(0);
534: }

536: /* ------------------------------------------------------------------- */
537: /*
538:    PreCheck - Optional user-defined routine that checks the validity of
539:    candidate steps of a line search method.  Set by SNESLineSearchSetPreCheck().

541:    Input Parameters:
542:    snes - the SNES context
543:    xcurrent - current solution
544:    y - search direction and length

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

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

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

570:    Output Parameters:
571:    y    - proposed step (search direction and length) (possibly changed)
572:    x    - current iterate (possibly modified)

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

587:   *changed_x = PETSC_FALSE;
588:   *changed_y = PETSC_FALSE;

590:   SNESLineSearchGetSNES(linesearch, &snes);
591:   check = (StepCheckCtx*)ctx;
592:   user  = check->user;
593:   SNESGetIterationNumber(snes,&iter);

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

600:     /* Access local array data */
601:     DMDAVecGetArray(da,check->last_step,&xa_last);
602:     DMDAVecGetArray(da,x,&xa);
603:     DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

605:     /*
606:        If we fail the user-defined check for validity of the candidate iterate,
607:        then modify the iterate as we like.  (Note that the particular modification
608:        below is intended simply to demonstrate how to manipulate this data, not
609:        as a meaningful or appropriate choice.)
610:     */
611:     for (i=xs; i<xs+xm; i++) {
612:       if (!PetscAbsScalar(xa[i])) rdiff = 2*check->tolerance;
613:       else rdiff = PetscAbsScalar((xa[i] - xa_last[i])/xa[i]);
614:       if (rdiff > check->tolerance) {
615:         tmp        = xa[i];
616:         xa[i]      = .5*(xa[i] + xa_last[i]);
617:         *changed_x = PETSC_TRUE;
618:         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]));
619:       }
620:     }
621:     DMDAVecRestoreArray(da,check->last_step,&xa_last);
622:     DMDAVecRestoreArray(da,x,&xa);
623:   }
624:   VecCopy(x,check->last_step);
625:   return(0);
626: }

628: /* ------------------------------------------------------------------- */
629: /*
630:    PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
631:    e.g,
632:      mpiexec -n 8 ./ex3 -nox -n 10000 -ksp_type fgmres -pc_type bjacobi -pc_bjacobi_blocks 4 -sub_ksp_type gmres -sub_ksp_max_it 3 -post_setsubksp -sub_ksp_rtol 1.e-16
633:    Set by SNESLineSearchSetPostCheck().

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

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

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

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

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

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

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

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


700: /*TEST

702:    test:
703:       args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

705:    test:
706:       suffix: 2
707:       nsize: 3
708:       args: -nox -pc_type asm -mat_type mpiaij -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

710:    test:
711:       suffix: 3
712:       nsize: 2
713:       args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

715:    test:
716:       suffix: 4
717:       args: -nox -pre_check_iterates -post_check_iterates

719:    test:
720:       suffix: 5
721:       requires: double !complex !single
722:       nsize: 2
723:       args: -nox -snes_test_jacobian  -snes_test_jacobian_view

725:    test:
726:       suffix: 6
727:       requires: double !complex !single
728:       nsize: 4
729:       args: -test_jacobian_domain_error -snes_converged_reason -snes_check_jacobian_domain_error 1

731:    test:
732:       suffix: 7
733:       requires: double !complex !single
734:       nsize: 4
735:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontr -snes_check_jacobian_domain_error 1

737:    test:
738:       suffix: 8
739:       requires: double !complex !single
740:       nsize: 4
741:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonrsls -snes_check_jacobian_domain_error 1

743:    test:
744:       suffix: 9
745:       requires: double !complex !single
746:       nsize: 4
747:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type vinewtonssls -snes_check_jacobian_domain_error 1

749:    test:
750:       suffix: 10
751:       requires: double !complex !single
752:       nsize: 4
753:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type qn -snes_qn_scale_type jacobian -snes_check_jacobian_domain_error 1

755:    test:
756:       suffix: 11
757:       requires: double !complex !single
758:       nsize: 4
759:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type ms -snes_ms_type m62 -snes_ms_damping 0.9 -snes_check_jacobian_domain_error 1

761:    test:
762:       suffix: 12
763:       args: -view_initial
764:       filter: grep -v "type:"

766: TEST*/