Actual source code: snesut.c
petsc-3.4.5 2014-06-29
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: }