Actual source code: snesut.c

petsc-3.4.4 2014-03-13
  2: #include <petsc-private/snesimpl.h>       /*I   "petsc-private/snesimpl.h"   I*/
  3: #include <petscblaslapack.h>

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

 11:    Collective on SNES

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

 19:    Level: intermediate

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

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

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

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

 48:    Collective on SNES

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

 56:    Level: intermediate

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

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

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

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

 85:    Collective on SNES

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

 93:    Level: intermediate

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

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

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

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

121:    Collective on SNES

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

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

132:    Level: intermediate

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

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

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

152: PetscErrorCode SNESMonitorJacUpdateSpectrum(SNES snes,PetscInt it,PetscReal fnorm,void *ctx)
153: {
154: #if defined(PETSC_MISSING_LAPACK_GEEV)
155:   SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"GEEV - Lapack routine is unavailable\nNot able to provide eigen values.");
156: #elif defined(PETSC_HAVE_ESSL)
157:   SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_SUP,"GEEV - No support for ESSL Lapack Routines");
158: #else
159:   Vec            X;
160:   Mat            J,dJ,dJdense;
162:   PetscErrorCode (*func)(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
163:   PetscInt       n,i;
164:   PetscBLASInt   nb,lwork;
165:   PetscReal      *eigr,*eigi;
166:   MatStructure   flg = DIFFERENT_NONZERO_PATTERN;
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,&J,NULL,&func,NULL);
175:   MatDuplicate(J,MAT_COPY_VALUES,&dJ);
176:   SNESComputeJacobian(snes,X,&dJ,&dJ,&flg);
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:   PetscMalloc(n*sizeof(PetscReal),&eigr);
185:   PetscMalloc(n*sizeof(PetscReal),&eigi);
186:   PetscMalloc(lwork*sizeof(PetscScalar),&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,eigr[i],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,SNESMonitorRatioContext,&ctx);
373:   SNESGetConvergenceHistory(snes,&history,NULL,NULL);
374:   if (!history) {
375:     PetscMalloc(100*sizeof(PetscReal),&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,fgnorm);
403:   } else if (fgnorm > 1.e-11) {
404:     PetscViewerASCIIPrintf(viewer,"%3D SNES Function norm %5.3e \n",its,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:    SNESSkipConverged - 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  SNESSkipConverged(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: {

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

557:   VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);
558:   PetscLogObjectParents(snes,nw,snes->work);
559:   return(0);
560: }