Actual source code: iterativ.c
petsc-3.10.5 2019-03-28
1: /*
2: This file contains some simple default routines.
3: These routines should be SHORT, since they will be included in every
4: executable image that uses the iterative routines (note that, through
5: the registry system, we provide a way to load only the truely necessary
6: files)
7: */
8: #include <petsc/private/kspimpl.h>
9: #include <petscdmshell.h>
11: /*@
12: KSPGetResidualNorm - Gets the last (approximate preconditioned)
13: residual norm that has been computed.
15: Not Collective
17: Input Parameters:
18: . ksp - the iterative context
20: Output Parameters:
21: . rnorm - residual norm
23: Level: intermediate
25: .keywords: KSP, get, residual norm
27: .seealso: KSPBuildResidual()
28: @*/
29: PetscErrorCode KSPGetResidualNorm(KSP ksp,PetscReal *rnorm)
30: {
34: *rnorm = ksp->rnorm;
35: return(0);
36: }
38: /*@
39: KSPGetIterationNumber - Gets the current iteration number; if the
40: KSPSolve() is complete, returns the number of iterations
41: used.
43: Not Collective
45: Input Parameters:
46: . ksp - the iterative context
48: Output Parameters:
49: . its - number of iterations
51: Level: intermediate
53: Notes:
54: During the ith iteration this returns i-1
55: .keywords: KSP, get, residual norm
57: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetTotalIterations()
58: @*/
59: PetscErrorCode KSPGetIterationNumber(KSP ksp,PetscInt *its)
60: {
64: *its = ksp->its;
65: return(0);
66: }
68: /*@
69: KSPGetTotalIterations - Gets the total number of iterations this KSP object has performed since was created, counted over all linear solves
71: Not Collective
73: Input Parameters:
74: . ksp - the iterative context
76: Output Parameters:
77: . its - total number of iterations
79: Level: intermediate
81: Notes:
82: Use KSPGetIterationNumber() to get the count for the most recent solve only
83: If this is called within a linear solve (such as in a KSPMonitor routine) then it does not include iterations within that current solve
85: .keywords: KSP, get, residual norm
87: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetIterationNumber()
88: @*/
89: PetscErrorCode KSPGetTotalIterations(KSP ksp,PetscInt *its)
90: {
94: *its = ksp->totalits;
95: return(0);
96: }
98: /*@C
99: KSPMonitorSingularValue - Prints the two norm of the true residual and
100: estimation of the extreme singular values of the preconditioned problem
101: at each iteration.
103: Logically Collective on KSP
105: Input Parameters:
106: + ksp - the iterative context
107: . n - the iteration
108: - rnorm - the two norm of the residual
110: Options Database Key:
111: . -ksp_monitor_singular_value - Activates KSPMonitorSingularValue()
113: Notes:
114: The CG solver uses the Lanczos technique for eigenvalue computation,
115: while GMRES uses the Arnoldi technique; other iterative methods do
116: not currently compute singular values.
118: Level: intermediate
120: .keywords: KSP, CG, default, monitor, extreme, singular values, Lanczos, Arnoldi
122: .seealso: KSPComputeExtremeSingularValues()
123: @*/
124: PetscErrorCode KSPMonitorSingularValue(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
125: {
126: PetscReal emin,emax,c;
128: PetscViewer viewer = dummy->viewer;
133: PetscViewerPushFormat(viewer,dummy->format);
134: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
135: if (!ksp->calc_sings) {
136: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
137: } else {
138: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
139: c = emax/emin;
140: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e %% max %14.12e min %14.12e max/min %14.12e\n",n,(double)rnorm,(double)emax,(double)emin,(double)c);
141: }
142: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
143: PetscViewerPopFormat(viewer);
144: return(0);
145: }
147: /*@C
148: KSPMonitorSolution - Monitors progress of the KSP solvers by calling
149: VecView() for the approximate solution at each iteration.
151: Collective on KSP
153: Input Parameters:
154: + ksp - the KSP context
155: . its - iteration number
156: . fgnorm - 2-norm of residual (or gradient)
157: - dummy - a viewer
159: Level: intermediate
161: Notes:
162: For some Krylov methods such as GMRES constructing the solution at
163: each iteration is expensive, hence using this will slow the code.
165: .keywords: KSP, nonlinear, vector, monitor, view
167: .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView()
168: @*/
169: PetscErrorCode KSPMonitorSolution(KSP ksp,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *dummy)
170: {
172: Vec x;
173: PetscViewer viewer = dummy->viewer;
177: KSPBuildSolution(ksp,NULL,&x);
178: PetscViewerPushFormat(viewer,dummy->format);
179: VecView(x,viewer);
180: PetscViewerPopFormat(viewer);
181: return(0);
182: }
184: /*@C
185: KSPMonitorDefault - Print the residual norm at each iteration of an
186: iterative solver.
188: Collective on KSP
190: Input Parameters:
191: + ksp - iterative context
192: . n - iteration number
193: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
194: - dummy - an ASCII PetscViewer
196: Level: intermediate
198: .keywords: KSP, default, monitor, residual
200: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
201: @*/
202: PetscErrorCode KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
203: {
205: PetscViewer viewer = dummy->viewer;
209: PetscViewerPushFormat(viewer,dummy->format);
210: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
211: if (n == 0 && ((PetscObject)ksp)->prefix) {
212: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
213: }
214: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
215: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
216: PetscViewerPopFormat(viewer);
217: return(0);
218: }
220: /*@C
221: KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
222: residual norm at each iteration of an iterative solver.
224: Collective on KSP
226: Input Parameters:
227: + ksp - iterative context
228: . n - iteration number
229: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
230: - dummy - an ASCII PetscViewer
232: Options Database Key:
233: . -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()
235: Notes:
236: When using right preconditioning, these values are equivalent.
238: Level: intermediate
240: .keywords: KSP, default, monitor, residual
242: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
243: @*/
244: PetscErrorCode KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
245: {
247: Vec resid;
248: PetscReal truenorm,bnorm;
249: PetscViewer viewer = dummy->viewer;
250: char normtype[256];
254: PetscViewerPushFormat(viewer,dummy->format);
255: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
256: if (n == 0 && ((PetscObject)ksp)->prefix) {
257: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
258: }
259: KSPBuildResidual(ksp,NULL,NULL,&resid);
260: VecNorm(resid,NORM_2,&truenorm);
261: VecDestroy(&resid);
262: VecNorm(ksp->vec_rhs,NORM_2,&bnorm);
263: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
264: PetscStrtolower(normtype);
265: PetscViewerASCIIPrintf(viewer,"%3D KSP %s resid norm %14.12e true resid norm %14.12e ||r(i)||/||b|| %14.12e\n",n,normtype,(double)rnorm,(double)truenorm,(double)(truenorm/bnorm));
266: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
267: PetscViewerPopFormat(viewer);
268: return(0);
269: }
271: /*@C
272: KSPMonitorTrueResidualMaxNorm - Prints the true residual max norm each iteration of an iterative solver.
274: Collective on KSP
276: Input Parameters:
277: + ksp - iterative context
278: . n - iteration number
279: . rnorm - norm (preconditioned) residual value (may be estimated).
280: - dummy - an ASCII viewer
282: Options Database Key:
283: . -ksp_monitor_max - Activates KSPMonitorTrueResidualMaxNorm()
285: Notes:
286: This could be implemented (better) with a flag in ksp.
288: Level: intermediate
290: .keywords: KSP, default, monitor, residual
292: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualNorm()
293: @*/
294: PetscErrorCode KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
295: {
297: Vec resid;
298: PetscReal truenorm,bnorm;
299: PetscViewer viewer = dummy->viewer;
300: char normtype[256];
304: PetscViewerPushFormat(viewer,dummy->format);
305: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
306: if (n == 0 && ((PetscObject)ksp)->prefix) {
307: PetscViewerASCIIPrintf(viewer," Residual norms (max) for %s solve.\n",((PetscObject)ksp)->prefix);
308: }
309: KSPBuildResidual(ksp,NULL,NULL,&resid);
310: VecNorm(resid,NORM_INFINITY,&truenorm);
311: VecDestroy(&resid);
312: VecNorm(ksp->vec_rhs,NORM_INFINITY,&bnorm);
313: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
314: PetscStrtolower(normtype);
315: PetscViewerASCIIPrintf(viewer,"%3D KSP true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n",n,(double)truenorm,(double)(truenorm/bnorm));
316: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
317: PetscViewerPopFormat(viewer);
318: return(0);
319: }
321: PetscErrorCode KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
322: {
324: Vec resid;
325: PetscReal rmax,pwork;
326: PetscInt i,n,N;
327: const PetscScalar *r;
330: KSPBuildResidual(ksp,NULL,NULL,&resid);
331: VecNorm(resid,NORM_INFINITY,&rmax);
332: VecGetLocalSize(resid,&n);
333: VecGetSize(resid,&N);
334: VecGetArrayRead(resid,&r);
335: pwork = 0.0;
336: for (i=0; i<n; i++) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
337: VecRestoreArrayRead(resid,&r);
338: VecDestroy(&resid);
339: MPIU_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
340: *per = *per/N;
341: return(0);
342: }
344: /*@C
345: KSPMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.
347: Collective on KSP
349: Input Parameters:
350: + ksp - iterative context
351: . it - iteration number
352: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
353: - dummy - an ASCII viewer
355: Options Database Key:
356: . -ksp_monitor_range - Activates KSPMonitorRange()
358: Level: intermediate
360: .keywords: KSP, default, monitor, residual
362: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
363: @*/
364: PetscErrorCode KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,PetscViewerAndFormat *dummy)
365: {
366: PetscErrorCode ierr;
367: PetscReal perc,rel;
368: PetscViewer viewer = dummy->viewer;
369: /* should be in a MonitorRangeContext */
370: static PetscReal prev;
374: PetscViewerPushFormat(viewer,dummy->format);
375: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
376: if (!it) prev = rnorm;
377: if (it == 0 && ((PetscObject)ksp)->prefix) {
378: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
379: }
380: KSPMonitorRange_Private(ksp,it,&perc);
382: rel = (prev - rnorm)/prev;
383: prev = rnorm;
384: PetscViewerASCIIPrintf(viewer,"%3D KSP 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));
385: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
386: PetscViewerPopFormat(viewer);
387: return(0);
388: }
390: /*@C
391: KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
392: outer iteration in an adaptive way.
394: Collective on KSP
396: Input Parameters:
397: + ksp - iterative context
398: . n - iteration number (not used)
399: . fnorm - the current residual norm
400: . dummy - some context as a C struct. fields:
401: coef: a scaling coefficient. default 1.0. can be passed through
402: -sub_ksp_dynamic_tolerance_param
403: bnrm: norm of the right-hand side. store it to avoid repeated calculation
405: Notes:
406: This may be useful for a flexibly preconditioner Krylov method to
407: control the accuracy of the inner solves needed to gaurantee the
408: convergence of the outer iterations.
410: Level: advanced
412: .keywords: KSP, inner tolerance
414: .seealso: KSPMonitorDynamicToleranceDestroy()
415: @*/
416: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
417: {
419: PC pc;
420: PetscReal outer_rtol, outer_abstol, outer_dtol, inner_rtol;
421: PetscInt outer_maxits,nksp,first,i;
422: KSPDynTolCtx *scale = (KSPDynTolCtx*)dummy;
423: KSP *subksp = NULL;
424: PetscBool isksp;
427: KSPGetPC(ksp, &pc);
429: /* compute inner_rtol */
430: if (scale->bnrm < 0.0) {
431: Vec b;
432: KSPGetRhs(ksp, &b);
433: VecNorm(b, NORM_2, &(scale->bnrm));
434: }
435: KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
436: inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
437: /*PetscPrintf(PETSC_COMM_WORLD, " Inner rtol = %g\n", (double)inner_rtol);*/
439: /* if pc is ksp */
440: PetscObjectTypeCompare((PetscObject)pc,PCKSP,&isksp);
441: if (isksp) {
442: KSP kspinner;
444: PCKSPGetKSP(pc, &kspinner);
445: KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
446: return(0);
447: }
449: /* if pc is bjacobi */
450: PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
451: if (subksp) {
452: for (i=0; i<nksp; i++) {
453: KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
454: }
455: return(0);
456: }
458: /* todo: dynamic tolerance may apply to other types of pc too */
459: return(0);
460: }
462: /*
463: Destroy the dummy context used in KSPMonitorDynamicTolerance()
464: */
465: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
466: {
470: PetscFree(*dummy);
471: return(0);
472: }
474: /*
475: Default (short) KSP Monitor, same as KSPMonitorDefault() except
476: it prints fewer digits of the residual as the residual gets smaller.
477: This is because the later digits are meaningless and are often
478: different on different machines; by using this routine different
479: machines will usually generate the same output.
481: Deprecated: Intentionally has no manual page
482: */
483: PetscErrorCode KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *dummy)
484: {
486: PetscViewer viewer = dummy->viewer;
490: PetscViewerPushFormat(viewer,dummy->format);
491: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
492: if (its == 0 && ((PetscObject)ksp)->prefix) {
493: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
494: }
496: if (fnorm > 1.e-9) {
497: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %g \n",its,(double)fnorm);
498: } else if (fnorm > 1.e-11) {
499: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,(double)fnorm);
500: } else {
501: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);
502: }
503: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
504: PetscViewerPopFormat(viewer);
505: return(0);
506: }
508: /*@C
509: KSPConvergedSkip - Convergence test that do not return as converged
510: until the maximum number of iterations is reached.
512: Collective on KSP
514: Input Parameters:
515: + ksp - iterative context
516: . n - iteration number
517: . rnorm - 2-norm residual value (may be estimated)
518: - dummy - unused convergence context
520: Returns:
521: . reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS
523: Notes:
524: This should be used as the convergence test with the option
525: KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
526: not computed. Convergence is then declared after the maximum number
527: of iterations have been reached. Useful when one is using CG or
528: BiCGStab as a smoother.
530: Level: advanced
532: .keywords: KSP, default, convergence, residual
534: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
535: @*/
536: PetscErrorCode KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
537: {
541: *reason = KSP_CONVERGED_ITERATING;
542: if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
543: return(0);
544: }
547: /*@C
548: KSPConvergedDefaultCreate - Creates and initializes the space used by the KSPConvergedDefault() function context
550: Collective on KSP
552: Output Parameter:
553: . ctx - convergence context
555: Level: intermediate
557: .keywords: KSP, default, convergence, residual
559: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
560: KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
561: @*/
562: PetscErrorCode KSPConvergedDefaultCreate(void **ctx)
563: {
564: PetscErrorCode ierr;
565: KSPConvergedDefaultCtx *cctx;
568: PetscNew(&cctx);
569: *ctx = cctx;
570: return(0);
571: }
573: /*@
574: KSPConvergedDefaultSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
575: instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
576: is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.
578: Collective on KSP
580: Input Parameters:
581: . ksp - iterative context
583: Options Database:
584: . -ksp_converged_use_initial_residual_norm
586: Notes:
587: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
589: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
590: are defined in petscksp.h.
592: If the convergence test is not KSPConvergedDefault() then this is ignored.
594: If right preconditioning is being used then B does not appear in the above formula.
597: Level: intermediate
599: .keywords: KSP, default, convergence, residual
601: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm()
602: @*/
603: PetscErrorCode KSPConvergedDefaultSetUIRNorm(KSP ksp)
604: {
605: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
609: if (ksp->converged != KSPConvergedDefault) return(0);
610: if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
611: ctx->initialrtol = PETSC_TRUE;
612: return(0);
613: }
615: /*@
616: KSPConvergedDefaultSetUMIRNorm - makes the default convergence test use min(|| B*(b - A*(initial guess))||,|| B*b ||)
617: In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
618: is used there is no B in the above formula. UMIRNorm is short for Use Minimum Initial Residual Norm.
620: Collective on KSP
622: Input Parameters:
623: . ksp - iterative context
625: Options Database:
626: . -ksp_converged_use_min_initial_residual_norm
628: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
630: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
631: are defined in petscksp.h.
633: Level: intermediate
635: .keywords: KSP, default, convergence, residual
637: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm()
638: @*/
639: PetscErrorCode KSPConvergedDefaultSetUMIRNorm(KSP ksp)
640: {
641: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
645: if (ksp->converged != KSPConvergedDefault) return(0);
646: if (ctx->initialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
647: ctx->mininitialrtol = PETSC_TRUE;
648: return(0);
649: }
651: /*@C
652: KSPConvergedDefault - Determines convergence of the linear iterative solvers by default
654: Collective on KSP
656: Input Parameters:
657: + ksp - iterative context
658: . n - iteration number
659: . rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
660: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
662: Output Parameter:
663: + positive - if the iteration has converged;
664: . negative - if residual norm exceeds divergence threshold;
665: - 0 - otherwise.
667: Notes:
668: KSPConvergedDefault() reaches convergence when rnorm < MAX (rtol * rnorm_0, abstol);
669: Divergence is detected if rnorm > dtol * rnorm_0,
671: where:
672: + rtol = relative tolerance,
673: . abstol = absolute tolerance.
674: . dtol = divergence tolerance,
675: - rnorm_0 is the two norm of the right hand side (or the preconditioned norm, depending on what was set with
676: KSPSetNormType(). When initial guess is non-zero you
677: can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
678: as the starting point for relative norm convergence testing, that is as rnorm_0
680: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
682: Use KSPSetNormType() (or -ksp_norm_type <none,preconditioned,unpreconditioned,natural>) to change the norm used for computing rnorm
684: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
686: This routine is used by KSP by default so the user generally never needs call it directly.
688: Use KSPSetConvergenceTest() to provide your own test instead of using this one.
690: Level: intermediate
692: .keywords: KSP, default, convergence, residual
694: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
695: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
696: @*/
697: PetscErrorCode KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
698: {
699: PetscErrorCode ierr;
700: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
701: KSPNormType normtype;
706: *reason = KSP_CONVERGED_ITERATING;
708: KSPGetNormType(ksp,&normtype);
709: if (normtype == KSP_NORM_NONE) return(0);
711: if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
712: if (!n) {
713: /* if user gives initial guess need to compute norm of b */
714: if (!ksp->guess_zero && !cctx->initialrtol) {
715: PetscReal snorm = 0.0;
716: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
717: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
718: VecNorm(ksp->vec_rhs,NORM_2,&snorm); /* <- b'*b */
719: } else {
720: Vec z;
721: /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
722: VecDuplicate(ksp->vec_rhs,&z);
723: KSP_PCApply(ksp,ksp->vec_rhs,z);
724: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
725: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
726: VecNorm(z,NORM_2,&snorm); /* dp <- b'*B'*B*b */
727: } else if (ksp->normtype == KSP_NORM_NATURAL) {
728: PetscScalar norm;
729: PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
730: VecDot(ksp->vec_rhs,z,&norm);
731: snorm = PetscSqrtReal(PetscAbsScalar(norm)); /* dp <- b'*B*b */
732: }
733: VecDestroy(&z);
734: }
735: /* handle special case of zero RHS and nonzero guess */
736: if (!snorm) {
737: PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
738: snorm = rnorm;
739: }
740: if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
741: else ksp->rnorm0 = snorm;
742: } else {
743: ksp->rnorm0 = rnorm;
744: }
745: ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
746: }
748: if (n <= ksp->chknorm) return(0);
750: if (PetscIsInfOrNanReal(rnorm)) {
751: PCFailedReason pcreason;
752: PetscInt sendbuf,pcreason_max;
753: PCGetSetUpFailedReason(ksp->pc,&pcreason);
754: sendbuf = (PetscInt)pcreason;
755: MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp));
756: if (pcreason_max) {
757: *reason = KSP_DIVERGED_PCSETUP_FAILED;
758: VecSetInf(ksp->vec_sol);
759: PetscInfo(ksp,"Linear solver pcsetup fails, declaring divergence \n");
760: } else {
761: *reason = KSP_DIVERGED_NANORINF;
762: PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
763: }
764: } else if (rnorm <= ksp->ttol) {
765: if (rnorm < ksp->abstol) {
766: PetscInfo3(ksp,"Linear solver has converged. Residual norm %14.12e is less than absolute tolerance %14.12e at iteration %D\n",(double)rnorm,(double)ksp->abstol,n);
767: *reason = KSP_CONVERGED_ATOL;
768: } else {
769: if (cctx->initialrtol) {
770: PetscInfo4(ksp,"Linear solver has converged. Residual norm %14.12e is less than relative tolerance %14.12e times initial residual norm %14.12e at iteration %D\n",(double)rnorm,(double)ksp->rtol,(double)ksp->rnorm0,n);
771: } else {
772: PetscInfo4(ksp,"Linear solver has converged. Residual norm %14.12e is less than relative tolerance %14.12e times initial right hand side norm %14.12e at iteration %D\n",(double)rnorm,(double)ksp->rtol,(double)ksp->rnorm0,n);
773: }
774: *reason = KSP_CONVERGED_RTOL;
775: }
776: } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
777: PetscInfo3(ksp,"Linear solver is diverging. Initial right hand size norm %14.12e, current residual norm %14.12e at iteration %D\n",(double)ksp->rnorm0,(double)rnorm,n);
778: *reason = KSP_DIVERGED_DTOL;
779: }
780: return(0);
781: }
783: /*@C
784: KSPConvergedDefaultDestroy - Frees the space used by the KSPConvergedDefault() function context
786: Collective on KSP
788: Input Parameters:
789: . ctx - convergence context
791: Level: intermediate
793: .keywords: KSP, default, convergence, residual
795: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
796: KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
797: @*/
798: PetscErrorCode KSPConvergedDefaultDestroy(void *ctx)
799: {
800: PetscErrorCode ierr;
801: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
804: VecDestroy(&cctx->work);
805: PetscFree(ctx);
806: return(0);
807: }
809: /*
810: KSPBuildSolutionDefault - Default code to create/move the solution.
812: Input Parameters:
813: + ksp - iterative context
814: - v - pointer to the user's vector
816: Output Parameter:
817: . V - pointer to a vector containing the solution
819: Level: advanced
821: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
823: .keywords: KSP, build, solution, default
825: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
826: */
827: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
828: {
832: if (ksp->pc_side == PC_RIGHT) {
833: if (ksp->pc) {
834: if (v) {
835: KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
836: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
837: } else {
838: if (v) {
839: VecCopy(ksp->vec_sol,v); *V = v;
840: } else *V = ksp->vec_sol;
841: }
842: } else if (ksp->pc_side == PC_SYMMETRIC) {
843: if (ksp->pc) {
844: if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
845: if (v) {
846: PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
847: *V = v;
848: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
849: } else {
850: if (v) {
851: VecCopy(ksp->vec_sol,v); *V = v;
852: } else *V = ksp->vec_sol;
853: }
854: } else {
855: if (v) {
856: VecCopy(ksp->vec_sol,v); *V = v;
857: } else *V = ksp->vec_sol;
858: }
859: return(0);
860: }
862: /*
863: KSPBuildResidualDefault - Default code to compute the residual.
865: Input Parameters:
866: . ksp - iterative context
867: . t - pointer to temporary vector
868: . v - pointer to user vector
870: Output Parameter:
871: . V - pointer to a vector containing the residual
873: Level: advanced
875: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
877: .keywords: KSP, build, residual, default
879: .seealso: KSPBuildSolutionDefault()
880: */
881: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
882: {
884: Mat Amat,Pmat;
887: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
888: PCGetOperators(ksp->pc,&Amat,&Pmat);
889: KSPBuildSolution(ksp,t,NULL);
890: KSP_MatMult(ksp,Amat,t,v);
891: VecAYPX(v,-1.0,ksp->vec_rhs);
892: *V = v;
893: return(0);
894: }
896: /*@C
897: KSPCreateVecs - Gets a number of work vectors.
899: Input Parameters:
900: + ksp - iterative context
901: . rightn - number of right work vectors
902: - leftn - number of left work vectors to allocate
904: Output Parameter:
905: + right - the array of vectors created
906: - left - the array of left vectors
908: Note: The right vector has as many elements as the matrix has columns. The left
909: vector has as many elements as the matrix has rows.
911: The vectors are new vectors that are not owned by the KSP, they should be destroyed with calls to VecDestroyVecs() when no longer needed.
913: Developers Note: First tries to duplicate the rhs and solution vectors of the KSP, if they do not exist tries to get them from the matrix, if
914: that does not exist tries to get them from the DM (if it is provided).
916: Level: advanced
918: .seealso: MatCreateVecs(), VecDestroyVecs()
920: @*/
921: PetscErrorCode KSPCreateVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
922: {
924: Vec vecr = NULL,vecl = NULL;
925: PetscBool matset,pmatset;
926: Mat mat = NULL;
929: if (rightn) {
930: if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
931: if (ksp->vec_sol) vecr = ksp->vec_sol;
932: else {
933: if (ksp->pc) {
934: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
935: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
936: if (matset) {
937: PCGetOperators(ksp->pc,&mat,NULL);
938: MatCreateVecs(mat,&vecr,NULL);
939: } else if (pmatset) {
940: PCGetOperators(ksp->pc,NULL,&mat);
941: MatCreateVecs(mat,&vecr,NULL);
942: }
943: }
944: if (!vecr) {
945: if (ksp->dm) {
946: DMGetGlobalVector(ksp->dm,&vecr);
947: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
948: }
949: }
950: VecDuplicateVecs(vecr,rightn,right);
951: if (!ksp->vec_sol) {
952: if (mat) {
953: VecDestroy(&vecr);
954: } else if (ksp->dm) {
955: DMRestoreGlobalVector(ksp->dm,&vecr);
956: }
957: }
958: }
959: if (leftn) {
960: if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
961: if (ksp->vec_rhs) vecl = ksp->vec_rhs;
962: else {
963: if (ksp->pc) {
964: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
965: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
966: if (matset) {
967: PCGetOperators(ksp->pc,&mat,NULL);
968: MatCreateVecs(mat,NULL,&vecl);
969: } else if (pmatset) {
970: PCGetOperators(ksp->pc,NULL,&mat);
971: MatCreateVecs(mat,NULL,&vecl);
972: }
973: }
974: if (!vecl) {
975: if (ksp->dm) {
976: DMGetGlobalVector(ksp->dm,&vecl);
977: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
978: }
979: }
980: VecDuplicateVecs(vecl,leftn,left);
981: if (!ksp->vec_rhs) {
982: if (mat) {
983: VecDestroy(&vecl);
984: } else if (ksp->dm) {
985: DMRestoreGlobalVector(ksp->dm,&vecl);
986: }
987: }
988: }
989: return(0);
990: }
992: /*@C
993: KSPSetWorkVecs - Sets a number of work vectors into a KSP object
995: Collective on KSP
997: Input Parameters:
998: + ksp - iterative context
999: - nw - number of work vectors to allocate
1001: Level: developer
1003: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1005: @*/
1006: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1007: {
1011: VecDestroyVecs(ksp->nwork,&ksp->work);
1012: ksp->nwork = nw;
1013: KSPCreateVecs(ksp,nw,&ksp->work,0,NULL);
1014: PetscLogObjectParents(ksp,nw,ksp->work);
1015: return(0);
1016: }
1018: /*
1019: KSPDestroyDefault - Destroys a iterative context variable for methods with
1020: no separate context. Preferred calling sequence KSPDestroy().
1022: Input Parameter:
1023: . ksp - the iterative context
1025: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1027: */
1028: PetscErrorCode KSPDestroyDefault(KSP ksp)
1029: {
1034: PetscFree(ksp->data);
1035: return(0);
1036: }
1038: /*@
1039: KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.
1041: Not Collective
1043: Input Parameter:
1044: . ksp - the KSP context
1046: Output Parameter:
1047: . reason - negative value indicates diverged, positive value converged, see KSPConvergedReason
1049: Possible values for reason: See also manual page for each reason
1050: $ KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1051: $ KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1052: $ KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence test routine is set.
1053: $ KSP_CONVERGED_CG_NEG_CURVE (see note below)
1054: $ KSP_CONVERGED_CG_CONSTRAINED (see note below)
1055: $ KSP_CONVERGED_STEP_LENGTH (see note below)
1056: $ KSP_CONVERGED_ITERATING (returned if the solver is not yet finished)
1057: $ KSP_DIVERGED_ITS (required more than its to reach convergence)
1058: $ KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1059: $ KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1060: $ KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1061: $ KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial residual. Try a different preconditioner, or a different initial Level.)
1063: Options Database:
1064: . -ksp_converged_reason - prints the reason to standard out
1066: Notes:
1067: If this routine is called before or doing the KSPSolve() the value of KSP_CONVERGED_ITERATING is returned
1069: The values KSP_CONVERGED_CG_NEG_CURVE, KSP_CONVERGED_CG_CONSTRAINED, and KSP_CONVERGED_STEP_LENGTH are returned only by the special KSPCGNASH, KSPCGSTCG, and KSPCGGLTR
1070: solvers which are used by the SNESNEWTONTR (trust region) solver.
1072: Level: intermediate
1074: .keywords: KSP, nonlinear, set, convergence, test
1076: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason
1077: @*/
1078: PetscErrorCode KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1079: {
1083: *reason = ksp->reason;
1084: return(0);
1085: }
1087: #include <petsc/private/dmimpl.h>
1088: /*@
1089: KSPSetDM - Sets the DM that may be used by some preconditioners
1091: Logically Collective on KSP
1093: Input Parameters:
1094: + ksp - the preconditioner context
1095: - dm - the dm, cannot be NULL
1097: Notes:
1098: If this is used then the KSP will attempt to use the DM to create the matrix and use the routine set with
1099: DMKSPSetComputeOperators(). Use KSPSetDMActive(ksp,PETSC_FALSE) to instead use the matrix you've provided with
1100: KSPSetOperators().
1102: A DM can only be used for solving one problem at a time because information about the problem is stored on the DM,
1103: even when not using interfaces like DMKSPSetComputeOperators(). Use DMClone() to get a distinct DM when solving
1104: different problems using the same function space.
1106: Level: intermediate
1108: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1109: @*/
1110: PetscErrorCode KSPSetDM(KSP ksp,DM dm)
1111: {
1113: PC pc;
1118: PetscObjectReference((PetscObject)dm);
1119: if (ksp->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
1120: if (ksp->dm->dmksp && !dm->dmksp) {
1121: DMKSP kdm;
1122: DMCopyDMKSP(ksp->dm,dm);
1123: DMGetDMKSP(ksp->dm,&kdm);
1124: if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1125: }
1126: DMDestroy(&ksp->dm);
1127: }
1128: ksp->dm = dm;
1129: ksp->dmAuto = PETSC_FALSE;
1130: KSPGetPC(ksp,&pc);
1131: PCSetDM(pc,dm);
1132: ksp->dmActive = PETSC_TRUE;
1133: return(0);
1134: }
1136: /*@
1137: KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side
1139: Logically Collective on KSP
1141: Input Parameters:
1142: + ksp - the preconditioner context
1143: - flg - use the DM
1145: Notes:
1146: By default KSPSetDM() sets the DM as active, call KSPSetDMActive(ksp,PETSC_FALSE); after KSPSetDM(ksp,dm) to not have the KSP object use the DM to generate the matrices.
1148: Level: intermediate
1150: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1151: @*/
1152: PetscErrorCode KSPSetDMActive(KSP ksp,PetscBool flg)
1153: {
1157: ksp->dmActive = flg;
1158: return(0);
1159: }
1161: /*@
1162: KSPGetDM - Gets the DM that may be used by some preconditioners
1164: Not Collective
1166: Input Parameter:
1167: . ksp - the preconditioner context
1169: Output Parameter:
1170: . dm - the dm
1172: Level: intermediate
1175: .seealso: KSPSetDM(), KSPSetDMActive()
1176: @*/
1177: PetscErrorCode KSPGetDM(KSP ksp,DM *dm)
1178: {
1183: if (!ksp->dm) {
1184: DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1185: ksp->dmAuto = PETSC_TRUE;
1186: }
1187: *dm = ksp->dm;
1188: return(0);
1189: }
1191: /*@
1192: KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.
1194: Logically Collective on KSP
1196: Input Parameters:
1197: + ksp - the KSP context
1198: - usrP - optional user context
1200: Fortran Notes:
1201: To use this from Fortran you must write a Fortran interface definition for this
1202: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1204: Level: intermediate
1206: .keywords: KSP, set, application, context
1208: .seealso: KSPGetApplicationContext()
1209: @*/
1210: PetscErrorCode KSPSetApplicationContext(KSP ksp,void *usrP)
1211: {
1213: PC pc;
1217: ksp->user = usrP;
1218: KSPGetPC(ksp,&pc);
1219: PCSetApplicationContext(pc,usrP);
1220: return(0);
1221: }
1223: /*@
1224: KSPGetApplicationContext - Gets the user-defined context for the linear solver.
1226: Not Collective
1228: Input Parameter:
1229: . ksp - KSP context
1231: Output Parameter:
1232: . usrP - user context
1234: Fortran Notes:
1235: To use this from Fortran you must write a Fortran interface definition for this
1236: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1238: Level: intermediate
1240: .keywords: KSP, get, application, context
1242: .seealso: KSPSetApplicationContext()
1243: @*/
1244: PetscErrorCode KSPGetApplicationContext(KSP ksp,void *usrP)
1245: {
1248: *(void**)usrP = ksp->user;
1249: return(0);
1250: }