Actual source code: ex3.c

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

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

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

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

 80: int main(int argc,char **argv)
 81: {
 82:   SNES           snes;                 /* SNES context */
 83:   SNESLineSearch linesearch;           /* SNESLineSearch context */
 84:   Mat            J;                    /* Jacobian matrix */
 85:   ApplicationCtx ctx;                  /* user-defined context */
 86:   Vec            x,r,U,F;              /* vectors */
 87:   MonitorCtx     monP;                 /* monitoring context */
 88:   StepCheckCtx   checkP;               /* step-checking context */
 89:   SetSubKSPCtx   checkP1;
 90:   PetscBool      pre_check,post_check,post_setsubksp; /* flag indicating whether we're checking candidate iterates */
 91:   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,viewinitial = PETSC_FALSE;

 96:   PetscInitialize(&argc,&argv,(char*)0,help);
 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);
101:   ctx.sjerr = PETSC_FALSE;
102:   PetscOptionsGetBool(NULL,NULL,"-test_jacobian_domain_error",&ctx.sjerr,NULL);
103:   PetscOptionsGetBool(NULL,NULL,"-view_initial",&viewinitial,NULL);

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

109:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

115:   /*
116:      Create distributed array (DMDA) to manage parallel grid and vectors
117:   */
118:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,N,1,1,NULL,&ctx.da);
119:   DMSetFromOptions(ctx.da);
120:   DMSetUp(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,NULL);
149:   MatMPIAIJSetPreallocation(J,3,NULL,3,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:      Optionally allow user-provided preconditioner
163:    */
164:   PetscOptionsHasName(NULL,NULL,"-user_precond",&flg);
165:   if (flg) {
166:     KSP ksp;
167:     PC  pc;
168:     SNESGetKSP(snes,&ksp);
169:     KSPGetPC(ksp,&pc);
170:     PCSetType(pc,PCSHELL);
171:     PCShellSetApply(pc,MatrixFreePreconditioner);
172:   }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

261:   /*
262:      Restore vectors
263:   */
264:   DMDAVecRestoreArray(ctx.da,F,&FF);
265:   DMDAVecRestoreArray(ctx.da,U,&UU);
266:   if (viewinitial) {
267:     VecView(U,0);
268:     VecView(F,0);
269:   }

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

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

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

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

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

319: /* ------------------------------------------------------------------- */
320: /*
321:    FormInitialGuess - Computes initial guess.

323:    Input/Output Parameter:
324: .  x - the solution vector
325: */
326: PetscErrorCode FormInitialGuess(Vec x)
327: {
328:   PetscScalar    pfive = .50;

331:   VecSet(x,pfive);
332:   return 0;
333: }

335: /* ------------------------------------------------------------------- */
336: /*
337:    FormFunction - Evaluates nonlinear function, F(x).

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

344:    Output Parameter:
345: .  f - function vector

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

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

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

380:   /*
381:      Get local grid boundaries (for 1-dimensional DMDA):
382:        xs, xm  - starting grid index, width of local grid (no ghost points)
383:   */
384:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);
385:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

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

401:   /*
402:      Compute function over locally owned part of the grid (interior points only)
403:   */
404:   d = 1.0/(user->h*user->h);
405:   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];

407:   /*
408:      Restore vectors
409:   */
410:   DMDAVecRestoreArray(da,xlocal,&xx);
411:   DMDAVecRestoreArray(da,f,&ff);
412:   DMDAVecRestoreArray(da,user->F,&FF);
413:   DMRestoreLocalVector(da,&xlocal);
414:   return 0;
415: }

417: /* ------------------------------------------------------------------- */
418: /*
419:    FormJacobian - Evaluates Jacobian matrix.

421:    Input Parameters:
422: .  snes - the SNES context
423: .  x - input vector
424: .  dummy - optional user-defined context (not used here)

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

439:   if (user->sjerr) {
440:     SNESSetJacobianDomainError(snes);
441:     return 0;
442:   }
443:   /*
444:      Get pointer to vector data
445:   */
446:   DMDAVecGetArrayRead(da,x,&xx);
447:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

449:   /*
450:     Get range of locally owned matrix
451:   */
452:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

454:   /*
455:      Determine starting and ending local indices for interior grid points.
456:      Set Jacobian entries for boundary points.
457:   */

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

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

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

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

490:      Also, restore vector.
491:   */

493:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
494:   DMDAVecRestoreArrayRead(da,x,&xx);
495:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

497:   return 0;
498: }

500: /* ------------------------------------------------------------------- */
501: /*
502:    Monitor - Optional user-defined monitoring routine that views the
503:    current iterate with an x-window plot. Set by SNESMonitorSet().

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

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

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

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

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

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

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

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

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

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

578:   *changed_x = PETSC_FALSE;
579:   *changed_y = PETSC_FALSE;

581:   SNESLineSearchGetSNES(linesearch, &snes);
582:   check = (StepCheckCtx*)ctx;
583:   user  = check->user;
584:   SNESGetIterationNumber(snes,&iter);

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

591:     /* Access local array data */
592:     DMDAVecGetArray(da,check->last_step,&xa_last);
593:     DMDAVecGetArray(da,x,&xa);
594:     DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

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

619: /* ------------------------------------------------------------------- */
620: /*
621:    PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
622:    e.g,
623:      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
624:    Set by SNESLineSearchSetPostCheck().

626:    Input Parameters:
627:    linesearch - the LineSearch context
628:    xcurrent - current solution
629:    y - search direction and length
630:    x    - the new candidate iterate

632:    Output Parameters:
633:    y    - proposed step (search direction and length) (possibly changed)
634:    x    - current iterate (possibly modified)

636:  */
637: PetscErrorCode PostSetSubKSP(SNESLineSearch linesearch,Vec xcurrent,Vec y,Vec x,PetscBool  *changed_y,PetscBool  *changed_x, void * ctx)
638: {
639:   SetSubKSPCtx   *check;
640:   PetscInt       iter,its,sub_its,maxit;
641:   KSP            ksp,sub_ksp,*sub_ksps;
642:   PC             pc;
643:   PetscReal      ksp_ratio;
644:   SNES           snes;

647:   SNESLineSearchGetSNES(linesearch, &snes);
648:   check   = (SetSubKSPCtx*)ctx;
649:   SNESGetIterationNumber(snes,&iter);
650:   SNESGetKSP(snes,&ksp);
651:   KSPGetPC(ksp,&pc);
652:   PCBJacobiGetSubKSP(pc,NULL,NULL,&sub_ksps);
653:   sub_ksp = sub_ksps[0];
654:   KSPGetIterationNumber(ksp,&its);      /* outer KSP iteration number */
655:   KSPGetIterationNumber(sub_ksp,&sub_its); /* inner KSP iteration number */

657:   if (iter) {
658:     PetscPrintf(PETSC_COMM_WORLD,"    ...PostCheck snes iteration %D, ksp_it %D %D, subksp_it %D\n",iter,check->its0,its,sub_its);
659:     ksp_ratio = ((PetscReal)(its))/check->its0;
660:     maxit     = (PetscInt)(ksp_ratio*sub_its + 0.5);
661:     if (maxit < 2) maxit = 2;
662:     KSPSetTolerances(sub_ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,maxit);
663:     PetscPrintf(PETSC_COMM_WORLD,"    ...ksp_ratio %g, new maxit %D\n\n",(double)ksp_ratio,maxit);
664:   }
665:   check->its0 = its; /* save current outer KSP iteration number */
666:   return 0;
667: }

669: /* ------------------------------------------------------------------- */
670: /*
671:    MatrixFreePreconditioner - This routine demonstrates the use of a
672:    user-provided preconditioner.  This code implements just the null
673:    preconditioner, which of course is not recommended for general use.

675:    Input Parameters:
676: +  pc - preconditioner
677: -  x - input vector

679:    Output Parameter:
680: .  y - preconditioned vector
681: */
682: PetscErrorCode MatrixFreePreconditioner(PC pc,Vec x,Vec y)
683: {
684:   VecCopy(x,y);
685:   return 0;
686: }

688: /*TEST

690:    test:
691:       args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

693:    test:
694:       suffix: 2
695:       nsize: 3
696:       args: -nox -pc_type asm -mat_type mpiaij -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

698:    test:
699:       suffix: 3
700:       nsize: 2
701:       args: -nox -snes_monitor_cancel -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always

703:    test:
704:       suffix: 4
705:       args: -nox -pre_check_iterates -post_check_iterates

707:    test:
708:       suffix: 5
709:       requires: double !complex !single
710:       nsize: 2
711:       args: -nox -snes_test_jacobian  -snes_test_jacobian_view

713:    test:
714:       suffix: 6
715:       requires: double !complex !single
716:       nsize: 4
717:       args: -test_jacobian_domain_error -snes_converged_reason -snes_check_jacobian_domain_error 1

719:    test:
720:       suffix: 7
721:       requires: double !complex !single
722:       nsize: 4
723:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontr -snes_check_jacobian_domain_error 1

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

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

737:    test:
738:       suffix: 10
739:       requires: double !complex !single
740:       nsize: 4
741:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type qn -snes_qn_scale_type jacobian -snes_check_jacobian_domain_error 1

743:    test:
744:       suffix: 11
745:       requires: double !complex !single
746:       nsize: 4
747:       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

749:    test:
750:       suffix: 12
751:       args: -view_initial
752:       filter: grep -v "type:"

754:    test:
755:       suffix: 13
756:       requires: double !complex !single
757:       nsize: 4
758:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontrdc -snes_check_jacobian_domain_error 1

760: TEST*/