Actual source code: iterativ.c
petsc-3.9.4 2018-09-11
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: Use KSPGetIterationNumber() to get the count for the most recent solve only
82: If this is called within a linear solve (such as in a KSPMonitor routine) then it does not include iterations within that current solve
84: .keywords: KSP, get, residual norm
86: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetIterationNumber()
87: @*/
88: PetscErrorCode KSPGetTotalIterations(KSP ksp,PetscInt *its)
89: {
93: *its = ksp->totalits;
94: return(0);
95: }
97: /*@C
98: KSPMonitorSingularValue - Prints the two norm of the true residual and
99: estimation of the extreme singular values of the preconditioned problem
100: at each iteration.
102: Logically Collective on KSP
104: Input Parameters:
105: + ksp - the iterative context
106: . n - the iteration
107: - rnorm - the two norm of the residual
109: Options Database Key:
110: . -ksp_monitor_singular_value - Activates KSPMonitorSingularValue()
112: Notes:
113: The CG solver uses the Lanczos technique for eigenvalue computation,
114: while GMRES uses the Arnoldi technique; other iterative methods do
115: not currently compute singular values.
117: Level: intermediate
119: .keywords: KSP, CG, default, monitor, extreme, singular values, Lanczos, Arnoldi
121: .seealso: KSPComputeExtremeSingularValues()
122: @*/
123: PetscErrorCode KSPMonitorSingularValue(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
124: {
125: PetscReal emin,emax,c;
127: PetscViewer viewer = dummy->viewer;
132: PetscViewerPushFormat(viewer,dummy->format);
133: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
134: if (!ksp->calc_sings) {
135: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
136: } else {
137: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
138: c = emax/emin;
139: 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);
140: }
141: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
142: PetscViewerPopFormat(viewer);
143: return(0);
144: }
146: /*@C
147: KSPMonitorSolution - Monitors progress of the KSP solvers by calling
148: VecView() for the approximate solution at each iteration.
150: Collective on KSP
152: Input Parameters:
153: + ksp - the KSP context
154: . its - iteration number
155: . fgnorm - 2-norm of residual (or gradient)
156: - dummy - a viewer
158: Level: intermediate
160: Notes:
161: For some Krylov methods such as GMRES constructing the solution at
162: each iteration is expensive, hence using this will slow the code.
164: .keywords: KSP, nonlinear, vector, monitor, view
166: .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView()
167: @*/
168: PetscErrorCode KSPMonitorSolution(KSP ksp,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *dummy)
169: {
171: Vec x;
172: PetscViewer viewer = dummy->viewer;
176: KSPBuildSolution(ksp,NULL,&x);
177: PetscViewerPushFormat(viewer,dummy->format);
178: VecView(x,viewer);
179: PetscViewerPopFormat(viewer);
180: return(0);
181: }
183: /*@C
184: KSPMonitorDefault - Print the residual norm at each iteration of an
185: iterative solver.
187: Collective on KSP
189: Input Parameters:
190: + ksp - iterative context
191: . n - iteration number
192: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
193: - dummy - an ASCII PetscViewer
195: Level: intermediate
197: .keywords: KSP, default, monitor, residual
199: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
200: @*/
201: PetscErrorCode KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
202: {
204: PetscViewer viewer = dummy->viewer;
208: PetscViewerPushFormat(viewer,dummy->format);
209: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
210: if (n == 0 && ((PetscObject)ksp)->prefix) {
211: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
212: }
213: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
214: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
215: PetscViewerPopFormat(viewer);
216: return(0);
217: }
219: /*@C
220: KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
221: residual norm at each iteration of an iterative solver.
223: Collective on KSP
225: Input Parameters:
226: + ksp - iterative context
227: . n - iteration number
228: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
229: - dummy - an ASCII PetscViewer
231: Options Database Key:
232: . -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()
234: Notes:
235: When using right preconditioning, these values are equivalent.
237: Level: intermediate
239: .keywords: KSP, default, monitor, residual
241: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
242: @*/
243: PetscErrorCode KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
244: {
246: Vec resid;
247: PetscReal truenorm,bnorm;
248: PetscViewer viewer = dummy->viewer;
249: char normtype[256];
253: PetscViewerPushFormat(viewer,dummy->format);
254: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
255: if (n == 0 && ((PetscObject)ksp)->prefix) {
256: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
257: }
258: KSPBuildResidual(ksp,NULL,NULL,&resid);
259: VecNorm(resid,NORM_2,&truenorm);
260: VecDestroy(&resid);
261: VecNorm(ksp->vec_rhs,NORM_2,&bnorm);
262: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
263: PetscStrtolower(normtype);
264: 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));
265: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
266: PetscViewerPopFormat(viewer);
267: return(0);
268: }
270: /*@C
271: KSPMonitorTrueResidualMaxNorm - Prints the true residual max norm each iteration of an iterative solver.
273: Collective on KSP
275: Input Parameters:
276: + ksp - iterative context
277: . n - iteration number
278: . rnorm - norm (preconditioned) residual value (may be estimated).
279: - dummy - an ASCII viewer
281: Options Database Key:
282: . -ksp_monitor_max - Activates KSPMonitorTrueResidualMaxNorm()
284: Notes:
285: This could be implemented (better) with a flag in ksp.
287: Level: intermediate
289: .keywords: KSP, default, monitor, residual
291: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualNorm()
292: @*/
293: PetscErrorCode KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
294: {
296: Vec resid;
297: PetscReal truenorm,bnorm;
298: PetscViewer viewer = dummy->viewer;
299: char normtype[256];
303: PetscViewerPushFormat(viewer,dummy->format);
304: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
305: if (n == 0 && ((PetscObject)ksp)->prefix) {
306: PetscViewerASCIIPrintf(viewer," Residual norms (max) for %s solve.\n",((PetscObject)ksp)->prefix);
307: }
308: KSPBuildResidual(ksp,NULL,NULL,&resid);
309: VecNorm(resid,NORM_INFINITY,&truenorm);
310: VecDestroy(&resid);
311: VecNorm(ksp->vec_rhs,NORM_INFINITY,&bnorm);
312: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
313: PetscStrtolower(normtype);
314: PetscViewerASCIIPrintf(viewer,"%3D KSP true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n",n,(double)truenorm,(double)(truenorm/bnorm));
315: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
316: PetscViewerPopFormat(viewer);
317: return(0);
318: }
320: PetscErrorCode KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
321: {
323: Vec resid;
324: PetscReal rmax,pwork;
325: PetscInt i,n,N;
326: const PetscScalar *r;
329: KSPBuildResidual(ksp,NULL,NULL,&resid);
330: VecNorm(resid,NORM_INFINITY,&rmax);
331: VecGetLocalSize(resid,&n);
332: VecGetSize(resid,&N);
333: VecGetArrayRead(resid,&r);
334: pwork = 0.0;
335: for (i=0; i<n; i++) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
336: VecRestoreArrayRead(resid,&r);
337: VecDestroy(&resid);
338: MPIU_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
339: *per = *per/N;
340: return(0);
341: }
343: /*@C
344: KSPMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.
346: Collective on KSP
348: Input Parameters:
349: + ksp - iterative context
350: . it - iteration number
351: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
352: - dummy - an ASCII viewer
354: Options Database Key:
355: . -ksp_monitor_range - Activates KSPMonitorRange()
357: Level: intermediate
359: .keywords: KSP, default, monitor, residual
361: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
362: @*/
363: PetscErrorCode KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,PetscViewerAndFormat *dummy)
364: {
365: PetscErrorCode ierr;
366: PetscReal perc,rel;
367: PetscViewer viewer = dummy->viewer;
368: /* should be in a MonitorRangeContext */
369: static PetscReal prev;
373: PetscViewerPushFormat(viewer,dummy->format);
374: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
375: if (!it) prev = rnorm;
376: if (it == 0 && ((PetscObject)ksp)->prefix) {
377: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
378: }
379: KSPMonitorRange_Private(ksp,it,&perc);
381: rel = (prev - rnorm)/prev;
382: prev = rnorm;
383: 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));
384: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
385: PetscViewerPopFormat(viewer);
386: return(0);
387: }
389: /*@C
390: KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
391: outer iteration in an adaptive way.
393: Collective on KSP
395: Input Parameters:
396: + ksp - iterative context
397: . n - iteration number (not used)
398: . fnorm - the current residual norm
399: . dummy - some context as a C struct. fields:
400: coef: a scaling coefficient. default 1.0. can be passed through
401: -sub_ksp_dynamic_tolerance_param
402: bnrm: norm of the right-hand side. store it to avoid repeated calculation
404: Notes:
405: This may be useful for a flexibly preconditioner Krylov method to
406: control the accuracy of the inner solves needed to gaurantee the
407: convergence of the outer iterations.
409: Level: advanced
411: .keywords: KSP, inner tolerance
413: .seealso: KSPMonitorDynamicToleranceDestroy()
414: @*/
415: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
416: {
418: PC pc;
419: PetscReal outer_rtol, outer_abstol, outer_dtol, inner_rtol;
420: PetscInt outer_maxits,nksp,first,i;
421: KSPDynTolCtx *scale = (KSPDynTolCtx*)dummy;
422: KSP kspinner = NULL, *subksp = NULL;
425: KSPGetPC(ksp, &pc);
427: /* compute inner_rtol */
428: if (scale->bnrm < 0.0) {
429: Vec b;
430: KSPGetRhs(ksp, &b);
431: VecNorm(b, NORM_2, &(scale->bnrm));
432: }
433: KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
434: inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
435: /*PetscPrintf(PETSC_COMM_WORLD, " Inner rtol = %g\n", (double)inner_rtol);*/
437: /* if pc is ksp */
438: PCKSPGetKSP(pc, &kspinner);
439: if (kspinner) {
440: KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
441: return(0);
442: }
444: /* if pc is bjacobi */
445: PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
446: if (subksp) {
447: for (i=0; i<nksp; i++) {
448: KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
449: }
450: return(0);
451: }
453: /* todo: dynamic tolerance may apply to other types of pc too */
454: return(0);
455: }
457: /*
458: Destroy the dummy context used in KSPMonitorDynamicTolerance()
459: */
460: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
461: {
465: PetscFree(*dummy);
466: return(0);
467: }
469: /*
470: Default (short) KSP Monitor, same as KSPMonitorDefault() except
471: it prints fewer digits of the residual as the residual gets smaller.
472: This is because the later digits are meaningless and are often
473: different on different machines; by using this routine different
474: machines will usually generate the same output.
475: */
476: PetscErrorCode KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *dummy)
477: {
479: PetscViewer viewer = dummy->viewer;
483: PetscViewerPushFormat(viewer,dummy->format);
484: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
485: if (its == 0 && ((PetscObject)ksp)->prefix) {
486: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
487: }
489: if (fnorm > 1.e-9) {
490: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %g \n",its,(double)fnorm);
491: } else if (fnorm > 1.e-11) {
492: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,(double)fnorm);
493: } else {
494: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);
495: }
496: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
497: PetscViewerPopFormat(viewer);
498: return(0);
499: }
501: /*@C
502: KSPConvergedSkip - Convergence test that do not return as converged
503: until the maximum number of iterations is reached.
505: Collective on KSP
507: Input Parameters:
508: + ksp - iterative context
509: . n - iteration number
510: . rnorm - 2-norm residual value (may be estimated)
511: - dummy - unused convergence context
513: Returns:
514: . reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS
516: Notes:
517: This should be used as the convergence test with the option
518: KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
519: not computed. Convergence is then declared after the maximum number
520: of iterations have been reached. Useful when one is using CG or
521: BiCGStab as a smoother.
523: Level: advanced
525: .keywords: KSP, default, convergence, residual
527: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
528: @*/
529: PetscErrorCode KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
530: {
534: *reason = KSP_CONVERGED_ITERATING;
535: if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
536: return(0);
537: }
540: /*@C
541: KSPConvergedDefaultCreate - Creates and initializes the space used by the KSPConvergedDefault() function context
543: Collective on KSP
545: Output Parameter:
546: . ctx - convergence context
548: Level: intermediate
550: .keywords: KSP, default, convergence, residual
552: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
553: KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
554: @*/
555: PetscErrorCode KSPConvergedDefaultCreate(void **ctx)
556: {
557: PetscErrorCode ierr;
558: KSPConvergedDefaultCtx *cctx;
561: PetscNew(&cctx);
562: *ctx = cctx;
563: return(0);
564: }
566: /*@
567: KSPConvergedDefaultSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
568: instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
569: is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.
571: Collective on KSP
573: Input Parameters:
574: . ksp - iterative context
576: Options Database:
577: . -ksp_converged_use_initial_residual_norm
579: Notes:
580: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
582: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
583: are defined in petscksp.h.
585: If the convergence test is not KSPConvergedDefault() then this is ignored.
587: If right preconditioning is being used then B does not appear in the above formula.
590: Level: intermediate
592: .keywords: KSP, default, convergence, residual
594: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm()
595: @*/
596: PetscErrorCode KSPConvergedDefaultSetUIRNorm(KSP ksp)
597: {
598: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
602: if (ksp->converged != KSPConvergedDefault) return(0);
603: if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
604: ctx->initialrtol = PETSC_TRUE;
605: return(0);
606: }
608: /*@
609: KSPConvergedDefaultSetUMIRNorm - makes the default convergence test use min(|| B*(b - A*(initial guess))||,|| B*b ||)
610: In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
611: is used there is no B in the above formula. UMIRNorm is short for Use Minimum Initial Residual Norm.
613: Collective on KSP
615: Input Parameters:
616: . ksp - iterative context
618: Options Database:
619: . -ksp_converged_use_min_initial_residual_norm
621: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
623: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
624: are defined in petscksp.h.
626: Level: intermediate
628: .keywords: KSP, default, convergence, residual
630: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm()
631: @*/
632: PetscErrorCode KSPConvergedDefaultSetUMIRNorm(KSP ksp)
633: {
634: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
638: if (ksp->converged != KSPConvergedDefault) return(0);
639: if (ctx->initialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
640: ctx->mininitialrtol = PETSC_TRUE;
641: return(0);
642: }
644: /*@C
645: KSPConvergedDefault - Determines convergence of the linear iterative solvers by default
647: Collective on KSP
649: Input Parameters:
650: + ksp - iterative context
651: . n - iteration number
652: . rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
653: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
655: Output Parameter:
656: + positive - if the iteration has converged;
657: . negative - if residual norm exceeds divergence threshold;
658: - 0 - otherwise.
660: Notes:
661: KSPConvergedDefault() reaches convergence when rnorm < MAX (rtol * rnorm_0, abstol);
662: Divergence is detected if rnorm > dtol * rnorm_0,
664: where:
665: + rtol = relative tolerance,
666: . abstol = absolute tolerance.
667: . dtol = divergence tolerance,
668: - rnorm_0 is the two norm of the right hand side (or the preconditioned norm, depending on what was set with
669: KSPSetNormType(). When initial guess is non-zero you
670: can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
671: as the starting point for relative norm convergence testing, that is as rnorm_0
673: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
675: Use KSPSetNormType() (or -ksp_norm_type <none,preconditioned,unpreconditioned,natural>) to change the norm used for computing rnorm
677: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
679: This routine is used by KSP by default so the user generally never needs call it directly.
681: Use KSPSetConvergenceTest() to provide your own test instead of using this one.
683: Level: intermediate
685: .keywords: KSP, default, convergence, residual
687: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
688: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
689: @*/
690: PetscErrorCode KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
691: {
692: PetscErrorCode ierr;
693: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
694: KSPNormType normtype;
699: *reason = KSP_CONVERGED_ITERATING;
701: KSPGetNormType(ksp,&normtype);
702: if (normtype == KSP_NORM_NONE) return(0);
704: if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
705: if (!n) {
706: /* if user gives initial guess need to compute norm of b */
707: if (!ksp->guess_zero && !cctx->initialrtol) {
708: PetscReal snorm = 0.0;
709: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
710: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
711: VecNorm(ksp->vec_rhs,NORM_2,&snorm); /* <- b'*b */
712: } else {
713: Vec z;
714: /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
715: VecDuplicate(ksp->vec_rhs,&z);
716: KSP_PCApply(ksp,ksp->vec_rhs,z);
717: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
718: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
719: VecNorm(z,NORM_2,&snorm); /* dp <- b'*B'*B*b */
720: } else if (ksp->normtype == KSP_NORM_NATURAL) {
721: PetscScalar norm;
722: PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
723: VecDot(ksp->vec_rhs,z,&norm);
724: snorm = PetscSqrtReal(PetscAbsScalar(norm)); /* dp <- b'*B*b */
725: }
726: VecDestroy(&z);
727: }
728: /* handle special case of zero RHS and nonzero guess */
729: if (!snorm) {
730: PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
731: snorm = rnorm;
732: }
733: if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
734: else ksp->rnorm0 = snorm;
735: } else {
736: ksp->rnorm0 = rnorm;
737: }
738: ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
739: }
741: if (n <= ksp->chknorm) return(0);
743: if (PetscIsInfOrNanReal(rnorm)) {
744: PCFailedReason pcreason;
745: PetscInt sendbuf,pcreason_max;
746: PCGetSetUpFailedReason(ksp->pc,&pcreason);
747: sendbuf = (PetscInt)pcreason;
748: MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp));
749: if (pcreason_max) {
750: *reason = KSP_DIVERGED_PCSETUP_FAILED;
751: VecSetInf(ksp->vec_sol);
752: PetscInfo(ksp,"Linear solver pcsetup fails, declaring divergence \n");
753: } else {
754: *reason = KSP_DIVERGED_NANORINF;
755: PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
756: }
757: } else if (rnorm <= ksp->ttol) {
758: if (rnorm < ksp->abstol) {
759: 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);
760: *reason = KSP_CONVERGED_ATOL;
761: } else {
762: if (cctx->initialrtol) {
763: 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);
764: } else {
765: 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);
766: }
767: *reason = KSP_CONVERGED_RTOL;
768: }
769: } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
770: 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);
771: *reason = KSP_DIVERGED_DTOL;
772: }
773: return(0);
774: }
776: /*@C
777: KSPConvergedDefaultDestroy - Frees the space used by the KSPConvergedDefault() function context
779: Collective on KSP
781: Input Parameters:
782: . ctx - convergence context
784: Level: intermediate
786: .keywords: KSP, default, convergence, residual
788: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
789: KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
790: @*/
791: PetscErrorCode KSPConvergedDefaultDestroy(void *ctx)
792: {
793: PetscErrorCode ierr;
794: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
797: VecDestroy(&cctx->work);
798: PetscFree(ctx);
799: return(0);
800: }
802: /*
803: KSPBuildSolutionDefault - Default code to create/move the solution.
805: Input Parameters:
806: + ksp - iterative context
807: - v - pointer to the user's vector
809: Output Parameter:
810: . V - pointer to a vector containing the solution
812: Level: advanced
814: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
816: .keywords: KSP, build, solution, default
818: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
819: */
820: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
821: {
825: if (ksp->pc_side == PC_RIGHT) {
826: if (ksp->pc) {
827: if (v) {
828: KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
829: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
830: } else {
831: if (v) {
832: VecCopy(ksp->vec_sol,v); *V = v;
833: } else *V = ksp->vec_sol;
834: }
835: } else if (ksp->pc_side == PC_SYMMETRIC) {
836: if (ksp->pc) {
837: if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
838: if (v) {
839: PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
840: *V = v;
841: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
842: } else {
843: if (v) {
844: VecCopy(ksp->vec_sol,v); *V = v;
845: } else *V = ksp->vec_sol;
846: }
847: } else {
848: if (v) {
849: VecCopy(ksp->vec_sol,v); *V = v;
850: } else *V = ksp->vec_sol;
851: }
852: return(0);
853: }
855: /*
856: KSPBuildResidualDefault - Default code to compute the residual.
858: Input Parameters:
859: . ksp - iterative context
860: . t - pointer to temporary vector
861: . v - pointer to user vector
863: Output Parameter:
864: . V - pointer to a vector containing the residual
866: Level: advanced
868: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
870: .keywords: KSP, build, residual, default
872: .seealso: KSPBuildSolutionDefault()
873: */
874: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
875: {
877: Mat Amat,Pmat;
880: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
881: PCGetOperators(ksp->pc,&Amat,&Pmat);
882: KSPBuildSolution(ksp,t,NULL);
883: KSP_MatMult(ksp,Amat,t,v);
884: VecAYPX(v,-1.0,ksp->vec_rhs);
885: *V = v;
886: return(0);
887: }
889: /*@C
890: KSPCreateVecs - Gets a number of work vectors.
892: Input Parameters:
893: + ksp - iterative context
894: . rightn - number of right work vectors
895: - leftn - number of left work vectors to allocate
897: Output Parameter:
898: + right - the array of vectors created
899: - left - the array of left vectors
901: Note: The right vector has as many elements as the matrix has columns. The left
902: vector has as many elements as the matrix has rows.
904: The vectors are new vectors that are not owned by the KSP, they should be destroyed with calls to VecDestroyVecs() when no longer needed.
906: 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
907: that does not exist tries to get them from the DM (if it is provided).
909: Level: advanced
911: .seealso: MatCreateVecs(), VecDestroyVecs()
913: @*/
914: PetscErrorCode KSPCreateVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
915: {
917: Vec vecr = NULL,vecl = NULL;
918: PetscBool matset,pmatset;
919: Mat mat = NULL;
922: if (rightn) {
923: if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
924: if (ksp->vec_sol) vecr = ksp->vec_sol;
925: else {
926: if (ksp->pc) {
927: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
928: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
929: if (matset) {
930: PCGetOperators(ksp->pc,&mat,NULL);
931: MatCreateVecs(mat,&vecr,NULL);
932: } else if (pmatset) {
933: PCGetOperators(ksp->pc,NULL,&mat);
934: MatCreateVecs(mat,&vecr,NULL);
935: }
936: }
937: if (!vecr) {
938: if (ksp->dm) {
939: DMGetGlobalVector(ksp->dm,&vecr);
940: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
941: }
942: }
943: VecDuplicateVecs(vecr,rightn,right);
944: if (!ksp->vec_sol) {
945: if (mat) {
946: VecDestroy(&vecr);
947: } else if (ksp->dm) {
948: DMRestoreGlobalVector(ksp->dm,&vecr);
949: }
950: }
951: }
952: if (leftn) {
953: if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
954: if (ksp->vec_rhs) vecl = ksp->vec_rhs;
955: else {
956: if (ksp->pc) {
957: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
958: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
959: if (matset) {
960: PCGetOperators(ksp->pc,&mat,NULL);
961: MatCreateVecs(mat,NULL,&vecl);
962: } else if (pmatset) {
963: PCGetOperators(ksp->pc,NULL,&mat);
964: MatCreateVecs(mat,NULL,&vecl);
965: }
966: }
967: if (!vecl) {
968: if (ksp->dm) {
969: DMGetGlobalVector(ksp->dm,&vecl);
970: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
971: }
972: }
973: VecDuplicateVecs(vecl,leftn,left);
974: if (!ksp->vec_rhs) {
975: if (mat) {
976: VecDestroy(&vecl);
977: } else if (ksp->dm) {
978: DMRestoreGlobalVector(ksp->dm,&vecl);
979: }
980: }
981: }
982: return(0);
983: }
985: /*@C
986: KSPSetWorkVecs - Sets a number of work vectors into a KSP object
988: Collective on KSP
990: Input Parameters:
991: + ksp - iterative context
992: - nw - number of work vectors to allocate
994: Level: developer
996: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
998: @*/
999: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1000: {
1004: VecDestroyVecs(ksp->nwork,&ksp->work);
1005: ksp->nwork = nw;
1006: KSPCreateVecs(ksp,nw,&ksp->work,0,NULL);
1007: PetscLogObjectParents(ksp,nw,ksp->work);
1008: return(0);
1009: }
1011: /*
1012: KSPDestroyDefault - Destroys a iterative context variable for methods with
1013: no separate context. Preferred calling sequence KSPDestroy().
1015: Input Parameter:
1016: . ksp - the iterative context
1018: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1020: */
1021: PetscErrorCode KSPDestroyDefault(KSP ksp)
1022: {
1027: PetscFree(ksp->data);
1028: return(0);
1029: }
1031: /*@
1032: KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.
1034: Not Collective
1036: Input Parameter:
1037: . ksp - the KSP context
1039: Output Parameter:
1040: . reason - negative value indicates diverged, positive value converged, see KSPConvergedReason
1042: Possible values for reason: See also manual page for each reason
1043: $ KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1044: $ KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1045: $ KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence test routine is set.
1046: $ KSP_CONVERGED_CG_NEG_CURVE (see note below)
1047: $ KSP_CONVERGED_CG_CONSTRAINED (see note below)
1048: $ KSP_CONVERGED_STEP_LENGTH (see note below)
1049: $ KSP_CONVERGED_ITERATING (returned if the solver is not yet finished)
1050: $ KSP_DIVERGED_ITS (required more than its to reach convergence)
1051: $ KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1052: $ KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1053: $ KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1054: $ KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial residual. Try a different preconditioner, or a different initial Level.)
1056: Options Database:
1057: . -ksp_converged_reason - prints the reason to standard out
1059: Notes: If this routine is called before or doing the KSPSolve() the value of KSP_CONVERGED_ITERATING is returned
1061: 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
1062: solvers which are used by the SNESNEWTONTR (trust region) solver.
1064: Level: intermediate
1066: .keywords: KSP, nonlinear, set, convergence, test
1068: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason
1069: @*/
1070: PetscErrorCode KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1071: {
1075: *reason = ksp->reason;
1076: return(0);
1077: }
1079: #include <petsc/private/dmimpl.h>
1080: /*@
1081: KSPSetDM - Sets the DM that may be used by some preconditioners
1083: Logically Collective on KSP
1085: Input Parameters:
1086: + ksp - the preconditioner context
1087: - dm - the dm, cannot be NULL
1089: Notes: If this is used then the KSP will attempt to use the DM to create the matrix and use the routine
1090: set with DMKSPSetComputeOperators(). Use KSPSetDMActive(ksp,PETSC_FALSE) to instead use the matrix
1091: you've provided with KSPSetOperators().
1093: Level: intermediate
1095: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1096: @*/
1097: PetscErrorCode KSPSetDM(KSP ksp,DM dm)
1098: {
1100: PC pc;
1105: PetscObjectReference((PetscObject)dm);
1106: if (ksp->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
1107: if (ksp->dm->dmksp && !dm->dmksp) {
1108: DMKSP kdm;
1109: DMCopyDMKSP(ksp->dm,dm);
1110: DMGetDMKSP(ksp->dm,&kdm);
1111: if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1112: }
1113: DMDestroy(&ksp->dm);
1114: }
1115: ksp->dm = dm;
1116: ksp->dmAuto = PETSC_FALSE;
1117: KSPGetPC(ksp,&pc);
1118: PCSetDM(pc,dm);
1119: ksp->dmActive = PETSC_TRUE;
1120: return(0);
1121: }
1123: /*@
1124: KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side
1126: Logically Collective on KSP
1128: Input Parameters:
1129: + ksp - the preconditioner context
1130: - flg - use the DM
1132: Notes:
1133: 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.
1135: Level: intermediate
1137: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1138: @*/
1139: PetscErrorCode KSPSetDMActive(KSP ksp,PetscBool flg)
1140: {
1144: ksp->dmActive = flg;
1145: return(0);
1146: }
1148: /*@
1149: KSPGetDM - Gets the DM that may be used by some preconditioners
1151: Not Collective
1153: Input Parameter:
1154: . ksp - the preconditioner context
1156: Output Parameter:
1157: . dm - the dm
1159: Level: intermediate
1162: .seealso: KSPSetDM(), KSPSetDMActive()
1163: @*/
1164: PetscErrorCode KSPGetDM(KSP ksp,DM *dm)
1165: {
1170: if (!ksp->dm) {
1171: DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1172: ksp->dmAuto = PETSC_TRUE;
1173: }
1174: *dm = ksp->dm;
1175: return(0);
1176: }
1178: /*@
1179: KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.
1181: Logically Collective on KSP
1183: Input Parameters:
1184: + ksp - the KSP context
1185: - usrP - optional user context
1187: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
1188: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1190: Level: intermediate
1192: .keywords: KSP, set, application, context
1194: .seealso: KSPGetApplicationContext()
1195: @*/
1196: PetscErrorCode KSPSetApplicationContext(KSP ksp,void *usrP)
1197: {
1199: PC pc;
1203: ksp->user = usrP;
1204: KSPGetPC(ksp,&pc);
1205: PCSetApplicationContext(pc,usrP);
1206: return(0);
1207: }
1209: /*@
1210: KSPGetApplicationContext - Gets the user-defined context for the linear solver.
1212: Not Collective
1214: Input Parameter:
1215: . ksp - KSP context
1217: Output Parameter:
1218: . usrP - user context
1220: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
1221: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1223: Level: intermediate
1225: .keywords: KSP, get, application, context
1227: .seealso: KSPSetApplicationContext()
1228: @*/
1229: PetscErrorCode KSPGetApplicationContext(KSP ksp,void *usrP)
1230: {
1233: *(void**)usrP = ksp->user;
1234: return(0);
1235: }