Actual source code: snesut.c
petsc-3.5.4 2015-05-23
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: }