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: /*
 11:    Include "petscdm.h" so that we can use data management objects (DMs)
 12:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 13:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 14:    file automatically includes:
 15:      petscsys.h    - base PETSc routines
 16:      petscvec.h    - vectors
 17:      petscmat.h    - matrices
 18:      petscis.h     - index sets
 19:      petscksp.h    - Krylov subspace methods
 20:      petscviewer.h - viewers
 21:      petscpc.h     - preconditioners
 22:      petscksp.h    - linear solvers
 23: */

 25: #include <petscdm.h>
 26: #include <petscdmda.h>
 27: #include <petscsnes.h>

 29: /*
 30:    User-defined routines.
 31: */
 32: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 33: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 34: PetscErrorCode FormInitialGuess(Vec);
 35: PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
 36: PetscErrorCode PreCheck(SNESLineSearch,Vec,Vec,PetscBool*,void*);
 37: PetscErrorCode PostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
 38: PetscErrorCode PostSetSubKSP(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);
 39: PetscErrorCode MatrixFreePreconditioner(PC,Vec,Vec);

 41: /*
 42:    User-defined application context
 43: */
 44: typedef struct {
 45:   DM          da;      /* distributed array */
 46:   Vec         F;       /* right-hand-side of PDE */
 47:   PetscMPIInt rank;    /* rank of processor */
 48:   PetscMPIInt size;    /* size of communicator */
 49:   PetscReal   h;       /* mesh spacing */
 50:   PetscBool   sjerr;   /* if or not to test jacobian domain error */
 51: } ApplicationCtx;

 53: /*
 54:    User-defined context for monitoring
 55: */
 56: typedef struct {
 57:   PetscViewer viewer;
 58: } MonitorCtx;

 60: /*
 61:    User-defined context for checking candidate iterates that are
 62:    determined by line search methods
 63: */
 64: typedef struct {
 65:   Vec            last_step;  /* previous iterate */
 66:   PetscReal      tolerance;  /* tolerance for changes between successive iterates */
 67:   ApplicationCtx *user;
 68: } StepCheckCtx;

 70: typedef struct {
 71:   PetscInt its0; /* num of prevous outer KSP iterations */
 72: } SetSubKSPCtx;

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

 90:   PetscInitialize(&argc,&argv,(char*)0,help);
 91:   MPI_Comm_rank(PETSC_COMM_WORLD,&ctx.rank);
 92:   MPI_Comm_size(PETSC_COMM_WORLD,&ctx.size);
 93:   PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL);
 94:   ctx.h = 1.0/(N-1);
 95:   ctx.sjerr = PETSC_FALSE;
 96:   PetscOptionsGetBool(NULL,NULL,"-test_jacobian_domain_error",&ctx.sjerr,NULL);
 97:   PetscOptionsGetBool(NULL,NULL,"-view_initial",&viewinitial,NULL);

 99:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100:      Create nonlinear solver context
101:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

103:   SNESCreate(PETSC_COMM_WORLD,&snes);

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:      Create vector data structures; set function evaluation routine
107:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:      Create matrix data structure; set Jacobian evaluation routine
137:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

139:   MatCreate(PETSC_COMM_WORLD,&J);
140:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
141:   MatSetFromOptions(J);
142:   MatSeqAIJSetPreallocation(J,3,NULL);
143:   MatMPIAIJSetPreallocation(J,3,NULL,3,NULL);

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

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

168:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169:      Customize nonlinear solver; set runtime options
170:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

202:     checkP.tolerance = 1.0;
203:     checkP.user      = &ctx;

205:     PetscOptionsGetReal(NULL,NULL,"-check_tol",&checkP.tolerance,NULL);
206:   }

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

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

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

228:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229:      Initialize application:
230:      Store right-hand-side of PDE and exact solution
231:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

239:   /*
240:      Get pointers to vector data
241:   */
242:   DMDAVecGetArray(ctx.da,F,&FF);
243:   DMDAVecGetArray(ctx.da,U,&UU);

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

255:   /*
256:      Restore vectors
257:   */
258:   DMDAVecRestoreArray(ctx.da,F,&FF);
259:   DMDAVecRestoreArray(ctx.da,U,&UU);
260:   if (viewinitial) {
261:     VecView(U,0);
262:     VecView(F,0);
263:   }

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

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

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

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

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

313: /* ------------------------------------------------------------------- */
314: /*
315:    FormInitialGuess - Computes initial guess.

317:    Input/Output Parameter:
318: .  x - the solution vector
319: */
320: PetscErrorCode FormInitialGuess(Vec x)
321: {
322:   PetscScalar    pfive = .50;

325:   VecSet(x,pfive);
326:   return 0;
327: }

329: /* ------------------------------------------------------------------- */
330: /*
331:    FormFunction - Evaluates nonlinear function, F(x).

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

338:    Output Parameter:
339: .  f - function vector

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

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

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

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

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

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

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

411: /* ------------------------------------------------------------------- */
412: /*
413:    FormJacobian - Evaluates Jacobian matrix.

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

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

433:   if (user->sjerr) {
434:     SNESSetJacobianDomainError(snes);
435:     return 0;
436:   }
437:   /*
438:      Get pointer to vector data
439:   */
440:   DMDAVecGetArrayRead(da,x,&xx);
441:   DMDAGetCorners(da,&xs,NULL,NULL,&xm,NULL,NULL);

443:   /*
444:     Get range of locally owned matrix
445:   */
446:   DMDAGetInfo(da,NULL,&M,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);

448:   /*
449:      Determine starting and ending local indices for interior grid points.
450:      Set Jacobian entries for boundary points.
451:   */

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

456:     MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
457:     xs++;xm--;
458:   }
459:   if (xs+xm == M) { /* right boundary */
460:     i    = M-1;
461:     A[0] = 1.0;
462:     MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
463:     xm--;
464:   }

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

478:   /*
479:      Assemble matrix, using the 2-step process:
480:        MatAssemblyBegin(), MatAssemblyEnd().
481:      By placing code between these two statements, computations can be
482:      done while messages are in transition.

484:      Also, restore vector.
485:   */

487:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
488:   DMDAVecRestoreArrayRead(da,x,&xx);
489:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

491:   return 0;
492: }

494: /* ------------------------------------------------------------------- */
495: /*
496:    Monitor - Optional user-defined monitoring routine that views the
497:    current iterate with an x-window plot. Set by SNESMonitorSet().

499:    Input Parameters:
500:    snes - the SNES context
501:    its - iteration number
502:    norm - 2-norm function value (may be estimated)
503:    ctx - optional user-defined context for private data for the
504:          monitor routine, as set by SNESMonitorSet()

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

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

522: /* ------------------------------------------------------------------- */
523: /*
524:    PreCheck - Optional user-defined routine that checks the validity of
525:    candidate steps of a line search method.  Set by SNESLineSearchSetPreCheck().

527:    Input Parameters:
528:    snes - the SNES context
529:    xcurrent - current solution
530:    y - search direction and length

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

543: /* ------------------------------------------------------------------- */
544: /*
545:    PostCheck - Optional user-defined routine that checks the validity of
546:    candidate steps of a line search method.  Set by SNESLineSearchSetPostCheck().

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

556:    Output Parameters:
557:    y    - proposed step (search direction and length) (possibly changed)
558:    x    - current iterate (possibly modified)

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

572:   *changed_x = PETSC_FALSE;
573:   *changed_y = PETSC_FALSE;

575:   SNESLineSearchGetSNES(linesearch, &snes);
576:   check = (StepCheckCtx*)ctx;
577:   user  = check->user;
578:   SNESGetIterationNumber(snes,&iter);

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

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

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

613: /* ------------------------------------------------------------------- */
614: /*
615:    PostSetSubKSP - Optional user-defined routine that reset SubKSP options when hierarchical bjacobi PC is used
616:    e.g,
617:      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
618:    Set by SNESLineSearchSetPostCheck().

620:    Input Parameters:
621:    linesearch - the LineSearch context
622:    xcurrent - current solution
623:    y - search direction and length
624:    x    - the new candidate iterate

626:    Output Parameters:
627:    y    - proposed step (search direction and length) (possibly changed)
628:    x    - current iterate (possibly modified)

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

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

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

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

669:    Input Parameters:
670: +  pc - preconditioner
671: -  x - input vector

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

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:
702:       suffix: 5
703:       requires: double !complex !single
704:       nsize: 2
705:       args: -nox -snes_test_jacobian  -snes_test_jacobian_view

707:    test:
708:       suffix: 6
709:       requires: double !complex !single
710:       nsize: 4
711:       args: -test_jacobian_domain_error -snes_converged_reason -snes_check_jacobian_domain_error 1

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

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

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

731:    test:
732:       suffix: 10
733:       requires: double !complex !single
734:       nsize: 4
735:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type qn -snes_qn_scale_type jacobian -snes_check_jacobian_domain_error 1

737:    test:
738:       suffix: 11
739:       requires: double !complex !single
740:       nsize: 4
741:       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

743:    test:
744:       suffix: 12
745:       args: -view_initial
746:       filter: grep -v "type:"

748:    test:
749:       suffix: 13
750:       requires: double !complex !single
751:       nsize: 4
752:       args: -test_jacobian_domain_error -snes_converged_reason -snes_type newtontrdc -snes_check_jacobian_domain_error 1

754: TEST*/