Actual source code: snesut.c

petsc-3.5.1 2014-08-06
Report Typos and Errors
  2: #include <petsc-private/snesimpl.h>       /*I   "petsc-private/snesimpl.h"   I*/
  3: #include <petscdm.h>
  4: #include <petscblaslapack.h>

  8: /*@C
  9:    SNESMonitorSolution - Monitors progress of the SNES solvers by calling
 10:    VecView() for the approximate solution at each iteration.

 12:    Collective on SNES

 14:    Input Parameters:
 15: +  snes - the SNES context
 16: .  its - iteration number
 17: .  fgnorm - 2-norm of residual
 18: -  dummy - either a viewer or NULL

 20:    Level: intermediate

 22: .keywords: SNES, nonlinear, vector, monitor, view

 24: .seealso: SNESMonitorSet(), SNESMonitorDefault(), VecView()
 25: @*/
 26: PetscErrorCode  SNESMonitorSolution(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
 27: {
 29:   Vec            x;
 30:   PetscViewer    viewer = (PetscViewer) dummy;

 33:   SNESGetSolution(snes,&x);
 34:   if (!viewer) {
 35:     MPI_Comm comm;
 36:     PetscObjectGetComm((PetscObject)snes,&comm);
 37:     viewer = PETSC_VIEWER_DRAW_(comm);
 38:   }
 39:   VecView(x,viewer);
 40:   return(0);
 41: }

 45: /*@C
 46:    SNESMonitorResidual - Monitors progress of the SNES solvers by calling
 47:    VecView() for the residual at each iteration.

 49:    Collective on SNES

 51:    Input Parameters:
 52: +  snes - the SNES context
 53: .  its - iteration number
 54: .  fgnorm - 2-norm of residual
 55: -  dummy - either a viewer or NULL

 57:    Level: intermediate

 59: .keywords: SNES, nonlinear, vector, monitor, view

 61: .seealso: SNESMonitorSet(), SNESMonitorDefault(), VecView()
 62: @*/
 63: PetscErrorCode  SNESMonitorResidual(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
 64: {
 66:   Vec            x;
 67:   PetscViewer    viewer = (PetscViewer) dummy;

 70:   SNESGetFunction(snes,&x,0,0);
 71:   if (!viewer) {
 72:     MPI_Comm comm;
 73:     PetscObjectGetComm((PetscObject)snes,&comm);
 74:     viewer = PETSC_VIEWER_DRAW_(comm);
 75:   }
 76:   VecView(x,viewer);
 77:   return(0);
 78: }

 82: /*@C
 83:    SNESMonitorSolutionUpdate - Monitors progress of the SNES solvers by calling
 84:    VecView() for the UPDATE to the solution at each iteration.

 86:    Collective on SNES

 88:    Input Parameters:
 89: +  snes - the SNES context
 90: .  its - iteration number
 91: .  fgnorm - 2-norm of residual
 92: -  dummy - either a viewer or NULL

 94:    Level: intermediate

 96: .keywords: SNES, nonlinear, vector, monitor, view

 98: .seealso: SNESMonitorSet(), SNESMonitorDefault(), VecView()
 99: @*/
100: PetscErrorCode  SNESMonitorSolutionUpdate(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
101: {
103:   Vec            x;
104:   PetscViewer    viewer = (PetscViewer) dummy;

107:   SNESGetSolutionUpdate(snes,&x);
108:   if (!viewer) {
109:     MPI_Comm comm;
110:     PetscObjectGetComm((PetscObject)snes,&comm);
111:     viewer = PETSC_VIEWER_DRAW_(comm);
112:   }
113:   VecView(x,viewer);
114:   return(0);
115: }

119: /*@C
120:    SNESMonitorDefault - Monitors progress of the SNES solvers (default).

122:    Collective on SNES

124:    Input Parameters:
125: +  snes - the SNES context
126: .  its - iteration number
127: .  fgnorm - 2-norm of residual
128: -  dummy - unused context

130:    Notes:
131:    This routine prints the residual norm at each iteration.

133:    Level: intermediate

135: .keywords: SNES, nonlinear, default, monitor, norm

137: .seealso: SNESMonitorSet(), SNESMonitorSolution()
138: @*/
139: PetscErrorCode  SNESMonitorDefault(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
140: {
142:   PetscViewer    viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));

145:   PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
146:   PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm %14.12e \n",its,(double)fgnorm);
147:   PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
148:   return(0);
149: }

153: PetscErrorCode SNESMonitorJacUpdateSpectrum(SNES snes,PetscInt it,PetscReal fnorm,void *ctx)
154: {
155: #if defined(PETSC_MISSING_LAPACK_GEEV)
156:   SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
157: #elif defined(PETSC_HAVE_ESSL)
158:   SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"GEEV - No support for ESSL Lapack Routines");
159: #else
160:   Vec            X;
161:   Mat            J,dJ,dJdense;
163:   PetscErrorCode (*func)(SNES,Vec,Mat,Mat,void*);
164:   PetscInt       n,i;
165:   PetscBLASInt   nb,lwork;
166:   PetscReal      *eigr,*eigi;
167:   PetscScalar    *work;
168:   PetscScalar    *a;

171:   if (it == 0) return(0);
172:   /* create the difference between the current update and the current jacobian */
173:   SNESGetSolution(snes,&X);
174:   SNESGetJacobian(snes,NULL,&J,&func,NULL);
175:   MatDuplicate(J,MAT_COPY_VALUES,&dJ);
176:   SNESComputeJacobian(snes,X,dJ,dJ);
177:   MatAXPY(dJ,-1.0,J,SAME_NONZERO_PATTERN);

179:   /* compute the spectrum directly */
180:   MatConvert(dJ,MATSEQDENSE,MAT_INITIAL_MATRIX,&dJdense);
181:   MatGetSize(dJ,&n,NULL);
182:   PetscBLASIntCast(n,&nb);
183:   lwork = 3*nb;
184:   PetscMalloc1(n,&eigr);
185:   PetscMalloc1(n,&eigi);
186:   PetscMalloc1(lwork,&work);
187:   MatDenseGetArray(dJdense,&a);
188: #if !defined(PETSC_USE_COMPLEX)
189:   {
190:     PetscBLASInt lierr;
191:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
192:     PetscStackCallBLAS("LAPACKgeev",LAPACKgeev_("N","N",&nb,a,&nb,eigr,eigi,NULL,&nb,NULL,&nb,work,&lwork,&lierr));
193:     if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"geev() error %d",lierr);
194:     PetscFPTrapPop();
195:   }
196: #else
197:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex");
198: #endif
199:   PetscPrintf(PetscObjectComm((PetscObject)snes),"Eigenvalues of J_%d - J_%d:\n",it,it-1);
200:   for (i=0;i<n;i++) {
201:     PetscPrintf(PetscObjectComm((PetscObject)snes),"%5d: %20.5g + %20.5gi\n",i,(double)eigr[i],(double)eigi[i]);
202:   }
203:   MatDenseRestoreArray(dJdense,&a);
204:   MatDestroy(&dJ);
205:   MatDestroy(&dJdense);
206:   PetscFree(eigr);
207:   PetscFree(eigi);
208:   PetscFree(work);
209:   return(0);
210: #endif
211: }

215: PetscErrorCode  SNESMonitorRange_Private(SNES snes,PetscInt it,PetscReal *per)
216: {
218:   Vec            resid;
219:   PetscReal      rmax,pwork;
220:   PetscInt       i,n,N;
221:   PetscScalar    *r;

224:   SNESGetFunction(snes,&resid,0,0);
225:   VecNorm(resid,NORM_INFINITY,&rmax);
226:   VecGetLocalSize(resid,&n);
227:   VecGetSize(resid,&N);
228:   VecGetArray(resid,&r);
229:   pwork = 0.0;
230:   for (i=0; i<n; i++) {
231:     pwork += (PetscAbsScalar(r[i]) > .20*rmax);
232:   }
233:   MPI_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)snes));
234:   VecRestoreArray(resid,&r);
235:   *per = *per/N;
236:   return(0);
237: }

241: /*@C
242:    SNESMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.

244:    Collective on SNES

246:    Input Parameters:
247: +  snes   - iterative context
248: .  it    - iteration number
249: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
250: -  dummy - unused monitor context

252:    Options Database Key:
253: .  -snes_monitor_range - Activates SNESMonitorRange()

255:    Level: intermediate

257: .keywords: SNES, default, monitor, residual

259: .seealso: SNESMonitorSet(), SNESMonitorDefault(), SNESMonitorLGCreate()
260: @*/
261: PetscErrorCode  SNESMonitorRange(SNES snes,PetscInt it,PetscReal rnorm,void *dummy)
262: {
264:   PetscReal      perc,rel;
265:   PetscViewer    viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
266:   /* should be in a MonitorRangeContext */
267:   static PetscReal prev;

270:   if (!it) prev = rnorm;
271:   SNESMonitorRange_Private(snes,it,&perc);

273:   rel  = (prev - rnorm)/prev;
274:   prev = rnorm;
275:   PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
276:   PetscViewerASCIIPrintf(viewer,"%3D SNES preconditioned resid norm %14.12e Percent values above 20 percent of maximum %5.2f relative decrease %5.2e ratio %5.2e \n",it,(double)rnorm,(double)(100.0*perc),(double)rel,(double)(rel/perc));
277:   PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
278:   return(0);
279: }

281: typedef struct {
282:   PetscViewer viewer;
283:   PetscReal   *history;
284: } SNESMonitorRatioContext;

288: /*@C
289:    SNESMonitorRatio - Monitors progress of the SNES solvers by printing the ratio
290:    of residual norm at each iteration to the previous.

292:    Collective on SNES

294:    Input Parameters:
295: +  snes - the SNES context
296: .  its - iteration number
297: .  fgnorm - 2-norm of residual (or gradient)
298: -  dummy -  context of monitor

300:    Level: intermediate

302: .keywords: SNES, nonlinear, monitor, norm

304: .seealso: SNESMonitorSet(), SNESMonitorSolution()
305: @*/
306: PetscErrorCode  SNESMonitorRatio(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
307: {
308:   PetscErrorCode          ierr;
309:   PetscInt                len;
310:   PetscReal               *history;
311:   SNESMonitorRatioContext *ctx = (SNESMonitorRatioContext*)dummy;

314:   SNESGetConvergenceHistory(snes,&history,NULL,&len);
315:   PetscViewerASCIIAddTab(ctx->viewer,((PetscObject)snes)->tablevel);
316:   if (!its || !history || its > len) {
317:     PetscViewerASCIIPrintf(ctx->viewer,"%3D SNES Function norm %14.12e \n",its,(double)fgnorm);
318:   } else {
319:     PetscReal ratio = fgnorm/history[its-1];
320:     PetscViewerASCIIPrintf(ctx->viewer,"%3D SNES Function norm %14.12e %14.12e \n",its,(double)fgnorm,(double)ratio);
321:   }
322:   PetscViewerASCIISubtractTab(ctx->viewer,((PetscObject)snes)->tablevel);
323:   return(0);
324: }

326: /*
327:    If the we set the history monitor space then we must destroy it
328: */
331: PetscErrorCode SNESMonitorRatioDestroy(void **ct)
332: {
333:   PetscErrorCode          ierr;
334:   SNESMonitorRatioContext *ctx = *(SNESMonitorRatioContext**)ct;

337:   PetscFree(ctx->history);
338:   PetscViewerDestroy(&ctx->viewer);
339:   PetscFree(ctx);
340:   return(0);
341: }

345: /*@C
346:    SNESMonitorSetRatio - Sets SNES to use a monitor that prints the
347:    ratio of the function norm at each iteration.

349:    Collective on SNES

351:    Input Parameters:
352: +   snes - the SNES context
353: -   viewer - ASCII viewer to print output

355:    Level: intermediate

357: .keywords: SNES, nonlinear, monitor, norm

359: .seealso: SNESMonitorSet(), SNESMonitorSolution(), SNESMonitorDefault()
360: @*/
361: PetscErrorCode  SNESMonitorSetRatio(SNES snes,PetscViewer viewer)
362: {
363:   PetscErrorCode          ierr;
364:   SNESMonitorRatioContext *ctx;
365:   PetscReal               *history;

368:   if (!viewer) {
369:     PetscViewerASCIIOpen(PetscObjectComm((PetscObject)snes),"stdout",&viewer);
370:     PetscObjectReference((PetscObject)viewer);
371:   }
372:   PetscNewLog(snes,&ctx);
373:   SNESGetConvergenceHistory(snes,&history,NULL,NULL);
374:   if (!history) {
375:     PetscMalloc1(100,&ctx->history);
376:     SNESSetConvergenceHistory(snes,ctx->history,0,100,PETSC_TRUE);
377:   }
378:   ctx->viewer = viewer;

380:   SNESMonitorSet(snes,SNESMonitorRatio,ctx,SNESMonitorRatioDestroy);
381:   return(0);
382: }

384: /* ---------------------------------------------------------------- */
387: /*
388:      Default (short) SNES Monitor, same as SNESMonitorDefault() except
389:   it prints fewer digits of the residual as the residual gets smaller.
390:   This is because the later digits are meaningless and are often
391:   different on different machines; by using this routine different
392:   machines will usually generate the same output.
393: */
394: PetscErrorCode  SNESMonitorDefaultShort(SNES snes,PetscInt its,PetscReal fgnorm,void *dummy)
395: {
397:   PetscViewer    viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));

400:   PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
401:   if (fgnorm > 1.e-9) {
402:     PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm %g \n",its,(double)fgnorm);
403:   } else if (fgnorm > 1.e-11) {
404:     PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm %5.3e \n",its,(double)fgnorm);
405:   } else {
406:     PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm < 1.e-11\n",its);
407:   }
408:   PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);
409:   return(0);
410: }
411: /* ---------------------------------------------------------------- */
414: /*@C
415:    SNESConvergedDefault - Convergence test of the solvers for
416:    systems of nonlinear equations (default).

418:    Collective on SNES

420:    Input Parameters:
421: +  snes - the SNES context
422: .  it - the iteration (0 indicates before any Newton steps)
423: .  xnorm - 2-norm of current iterate
424: .  snorm - 2-norm of current step
425: .  fnorm - 2-norm of function at current iterate
426: -  dummy - unused context

428:    Output Parameter:
429: .   reason  - one of
430: $  SNES_CONVERGED_FNORM_ABS       - (fnorm < abstol),
431: $  SNES_CONVERGED_SNORM_RELATIVE  - (snorm < stol*xnorm),
432: $  SNES_CONVERGED_FNORM_RELATIVE  - (fnorm < rtol*fnorm0),
433: $  SNES_DIVERGED_FUNCTION_COUNT   - (nfct > maxf),
434: $  SNES_DIVERGED_FNORM_NAN        - (fnorm == NaN),
435: $  SNES_CONVERGED_ITERATING       - (otherwise),

437:    where
438: +    maxf - maximum number of function evaluations,
439:             set with SNESSetTolerances()
440: .    nfct - number of function evaluations,
441: .    abstol - absolute function norm tolerance,
442:             set with SNESSetTolerances()
443: -    rtol - relative function norm tolerance, set with SNESSetTolerances()

445:    Level: intermediate

447: .keywords: SNES, nonlinear, default, converged, convergence

449: .seealso: SNESSetConvergenceTest()
450: @*/
451: PetscErrorCode  SNESConvergedDefault(SNES snes,PetscInt it,PetscReal xnorm,PetscReal snorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
452: {


459:   *reason = SNES_CONVERGED_ITERATING;

461:   if (!it) {
462:     /* set parameter for default relative tolerance convergence test */
463:     snes->ttol = fnorm*snes->rtol;
464:   }
465:   if (PetscIsInfOrNanReal(fnorm)) {
466:     PetscInfo(snes,"Failed to converged, function norm is NaN\n");
467:     *reason = SNES_DIVERGED_FNORM_NAN;
468:   } else if (fnorm < snes->abstol) {
469:     PetscInfo2(snes,"Converged due to function norm %14.12e < %14.12e\n",(double)fnorm,(double)snes->abstol);
470:     *reason = SNES_CONVERGED_FNORM_ABS;
471:   } else if (snes->nfuncs >= snes->max_funcs) {
472:     PetscInfo2(snes,"Exceeded maximum number of function evaluations: %D > %D\n",snes->nfuncs,snes->max_funcs);
473:     *reason = SNES_DIVERGED_FUNCTION_COUNT;
474:   }

476:   if (it && !*reason) {
477:     if (fnorm <= snes->ttol) {
478:       PetscInfo2(snes,"Converged due to function norm %14.12e < %14.12e (relative tolerance)\n",(double)fnorm,(double)snes->ttol);
479:       *reason = SNES_CONVERGED_FNORM_RELATIVE;
480:     } else if (snorm < snes->stol*xnorm) {
481:       PetscInfo3(snes,"Converged due to small update length: %14.12e < %14.12e * %14.12e\n",(double)snorm,(double)snes->stol,(double)xnorm);
482:       *reason = SNES_CONVERGED_SNORM_RELATIVE;
483:     }
484:   }
485:   return(0);
486: }

490: /*@C
491:    SNESConvergedSkip - Convergence test for SNES that NEVER returns as
492:    converged, UNLESS the maximum number of iteration have been reached.

494:    Logically Collective on SNES

496:    Input Parameters:
497: +  snes - the SNES context
498: .  it - the iteration (0 indicates before any Newton steps)
499: .  xnorm - 2-norm of current iterate
500: .  snorm - 2-norm of current step
501: .  fnorm - 2-norm of function at current iterate
502: -  dummy - unused context

504:    Output Parameter:
505: .   reason  - SNES_CONVERGED_ITERATING, SNES_CONVERGED_ITS, or SNES_DIVERGED_FNORM_NAN

507:    Notes:
508:    Convergence is then declared after a fixed number of iterations have been used.

510:    Level: advanced

512: .keywords: SNES, nonlinear, skip, converged, convergence

514: .seealso: SNESSetConvergenceTest()
515: @*/
516: PetscErrorCode  SNESConvergedSkip(SNES snes,PetscInt it,PetscReal xnorm,PetscReal snorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
517: {


524:   *reason = SNES_CONVERGED_ITERATING;

526:   if (fnorm != fnorm) {
527:     PetscInfo(snes,"Failed to converged, function norm is NaN\n");
528:     *reason = SNES_DIVERGED_FNORM_NAN;
529:   } else if (it == snes->max_its) {
530:     *reason = SNES_CONVERGED_ITS;
531:   }
532:   return(0);
533: }

537: /*@C
538:   SNESSetWorkVecs - Gets a number of work vectors.

540:   Input Parameters:
541: . snes  - the SNES context
542: . nw - number of work vectors to allocate

544:    Level: developer

546:    Developers Note: This is PETSC_EXTERN because it may be used by user written plugin SNES implementations

548: @*/
549: PetscErrorCode SNESSetWorkVecs(SNES snes,PetscInt nw)
550: {
551:   DM             dm;
552:   Vec            v;

556:   if (snes->work) {VecDestroyVecs(snes->nwork,&snes->work);}
557:   snes->nwork = nw;

559:   SNESGetDM(snes, &dm);
560:   DMGetGlobalVector(dm, &v);
561:   VecDuplicateVecs(v,snes->nwork,&snes->work);
562:   DMRestoreGlobalVector(dm, &v);
563:   PetscLogObjectParents(snes,nw,snes->work);
564:   return(0);
565: }