Actual source code: iterativ.c
petsc-3.7.7 2017-09-25
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> /*I "petscksp.h" I*/
9: #include <petscdmshell.h>
13: /*@
14: KSPGetResidualNorm - Gets the last (approximate preconditioned)
15: residual norm that has been computed.
17: Not Collective
19: Input Parameters:
20: . ksp - the iterative context
22: Output Parameters:
23: . rnorm - residual norm
25: Level: intermediate
27: .keywords: KSP, get, residual norm
29: .seealso: KSPBuildResidual()
30: @*/
31: PetscErrorCode KSPGetResidualNorm(KSP ksp,PetscReal *rnorm)
32: {
36: *rnorm = ksp->rnorm;
37: return(0);
38: }
42: /*@
43: KSPGetIterationNumber - Gets the current iteration number; if the
44: KSPSolve() is complete, returns the number of iterations
45: used.
47: Not Collective
49: Input Parameters:
50: . ksp - the iterative context
52: Output Parameters:
53: . its - number of iterations
55: Level: intermediate
57: Notes:
58: During the ith iteration this returns i-1
59: .keywords: KSP, get, residual norm
61: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetTotalIterations()
62: @*/
63: PetscErrorCode KSPGetIterationNumber(KSP ksp,PetscInt *its)
64: {
68: *its = ksp->its;
69: return(0);
70: }
74: /*@
75: KSPGetTotalIterations - Gets the total number of iterations this KSP object has performed since was created, counted over all linear solves
77: Not Collective
79: Input Parameters:
80: . ksp - the iterative context
82: Output Parameters:
83: . its - total number of iterations
85: Level: intermediate
87: Notes: Use KSPGetIterationNumber() to get the count for the most recent solve only
88: If this is called within a linear solve (such as in a KSPMonitor routine) then it does not include iterations within that current solve
90: .keywords: KSP, get, residual norm
92: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetIterationNumber()
93: @*/
94: PetscErrorCode KSPGetTotalIterations(KSP ksp,PetscInt *its)
95: {
99: *its = ksp->totalits;
100: return(0);
101: }
105: /*@C
106: KSPMonitorSingularValue - Prints the two norm of the true residual and
107: estimation of the extreme singular values of the preconditioned problem
108: at each iteration.
110: Logically Collective on KSP
112: Input Parameters:
113: + ksp - the iterative context
114: . n - the iteration
115: - rnorm - the two norm of the residual
117: Options Database Key:
118: . -ksp_monitor_singular_value - Activates KSPMonitorSingularValue()
120: Notes:
121: The CG solver uses the Lanczos technique for eigenvalue computation,
122: while GMRES uses the Arnoldi technique; other iterative methods do
123: not currently compute singular values.
125: Level: intermediate
127: .keywords: KSP, CG, default, monitor, extreme, singular values, Lanczos, Arnoldi
129: .seealso: KSPComputeExtremeSingularValues()
130: @*/
131: PetscErrorCode KSPMonitorSingularValue(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
132: {
133: PetscReal emin,emax,c;
135: PetscViewer viewer = dummy->viewer;
140: PetscViewerPushFormat(viewer,dummy->format);
141: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
142: if (!ksp->calc_sings) {
143: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
144: } else {
145: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
146: c = emax/emin;
147: 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);
148: }
149: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
150: PetscViewerPopFormat(viewer);
151: return(0);
152: }
156: /*@C
157: KSPMonitorSolution - Monitors progress of the KSP solvers by calling
158: VecView() for the approximate solution at each iteration.
160: Collective on KSP
162: Input Parameters:
163: + ksp - the KSP context
164: . its - iteration number
165: . fgnorm - 2-norm of residual (or gradient)
166: - dummy - a viewer
168: Level: intermediate
170: Notes:
171: For some Krylov methods such as GMRES constructing the solution at
172: each iteration is expensive, hence using this will slow the code.
174: .keywords: KSP, nonlinear, vector, monitor, view
176: .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView()
177: @*/
178: PetscErrorCode KSPMonitorSolution(KSP ksp,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *dummy)
179: {
181: Vec x;
182: PetscViewer viewer = dummy->viewer;
186: KSPBuildSolution(ksp,NULL,&x);
187: PetscViewerPushFormat(viewer,dummy->format);
188: VecView(x,viewer);
189: PetscViewerPopFormat(viewer);
190: return(0);
191: }
195: /*@C
196: KSPMonitorDefault - Print the residual norm at each iteration of an
197: iterative solver.
199: Collective on KSP
201: Input Parameters:
202: + ksp - iterative context
203: . n - iteration number
204: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
205: - dummy - an ASCII PetscViewer
207: Level: intermediate
209: .keywords: KSP, default, monitor, residual
211: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
212: @*/
213: PetscErrorCode KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
214: {
216: PetscViewer viewer = dummy->viewer;
220: PetscViewerPushFormat(viewer,dummy->format);
221: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
222: if (n == 0 && ((PetscObject)ksp)->prefix) {
223: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
224: }
225: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
226: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
227: PetscViewerPopFormat(viewer);
228: return(0);
229: }
233: /*@C
234: KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
235: residual norm at each iteration of an iterative solver.
237: Collective on KSP
239: Input Parameters:
240: + ksp - iterative context
241: . n - iteration number
242: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
243: - dummy - an ASCII PetscViewer
245: Options Database Key:
246: . -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()
248: Notes:
249: When using right preconditioning, these values are equivalent.
251: Level: intermediate
253: .keywords: KSP, default, monitor, residual
255: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
256: @*/
257: PetscErrorCode KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
258: {
260: Vec resid;
261: PetscReal truenorm,bnorm;
262: PetscViewer viewer = dummy->viewer;
263: char normtype[256];
267: PetscViewerPushFormat(viewer,dummy->format);
268: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
269: if (n == 0 && ((PetscObject)ksp)->prefix) {
270: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
271: }
272: KSPBuildResidual(ksp,NULL,NULL,&resid);
273: VecNorm(resid,NORM_2,&truenorm);
274: VecDestroy(&resid);
275: VecNorm(ksp->vec_rhs,NORM_2,&bnorm);
276: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
277: PetscStrtolower(normtype);
278: 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));
279: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
280: PetscViewerPopFormat(viewer);
281: return(0);
282: }
286: /*@C
287: KSPMonitorTrueResidualMaxNorm - Prints the true residual max norm each iteration of an iterative solver.
289: Collective on KSP
291: Input Parameters:
292: + ksp - iterative context
293: . n - iteration number
294: . rnorm - norm (preconditioned) residual value (may be estimated).
295: - dummy - an ASCII viewer
297: Options Database Key:
298: . -ksp_monitor_max - Activates KSPMonitorTrueResidualMaxNorm()
300: Notes:
301: This could be implemented (better) with a flag in ksp.
303: Level: intermediate
305: .keywords: KSP, default, monitor, residual
307: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualNorm()
308: @*/
309: PetscErrorCode KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
310: {
312: Vec resid;
313: PetscReal truenorm,bnorm;
314: PetscViewer viewer = dummy->viewer;
315: char normtype[256];
319: PetscViewerPushFormat(viewer,dummy->format);
320: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
321: if (n == 0 && ((PetscObject)ksp)->prefix) {
322: PetscViewerASCIIPrintf(viewer," Residual norms (max) for %s solve.\n",((PetscObject)ksp)->prefix);
323: }
324: KSPBuildResidual(ksp,NULL,NULL,&resid);
325: VecNorm(resid,NORM_INFINITY,&truenorm);
326: VecDestroy(&resid);
327: VecNorm(ksp->vec_rhs,NORM_INFINITY,&bnorm);
328: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
329: PetscStrtolower(normtype);
330: PetscViewerASCIIPrintf(viewer,"%3D KSP true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n",n,(double)truenorm,(double)(truenorm/bnorm));
331: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
332: PetscViewerPopFormat(viewer);
333: return(0);
334: }
338: PetscErrorCode KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
339: {
341: Vec resid;
342: PetscReal rmax,pwork;
343: PetscInt i,n,N;
344: const PetscScalar *r;
347: KSPBuildResidual(ksp,NULL,NULL,&resid);
348: VecNorm(resid,NORM_INFINITY,&rmax);
349: VecGetLocalSize(resid,&n);
350: VecGetSize(resid,&N);
351: VecGetArrayRead(resid,&r);
352: pwork = 0.0;
353: for (i=0; i<n; i++) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
354: VecRestoreArrayRead(resid,&r);
355: VecDestroy(&resid);
356: MPIU_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
357: *per = *per/N;
358: return(0);
359: }
363: /*@C
364: KSPMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.
366: Collective on KSP
368: Input Parameters:
369: + ksp - iterative context
370: . it - iteration number
371: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
372: - dummy - an ASCII viewer
374: Options Database Key:
375: . -ksp_monitor_range - Activates KSPMonitorRange()
377: Level: intermediate
379: .keywords: KSP, default, monitor, residual
381: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
382: @*/
383: PetscErrorCode KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,PetscViewerAndFormat *dummy)
384: {
385: PetscErrorCode ierr;
386: PetscReal perc,rel;
387: PetscViewer viewer = dummy->viewer;
388: /* should be in a MonitorRangeContext */
389: static PetscReal prev;
393: PetscViewerPushFormat(viewer,dummy->format);
394: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
395: if (!it) prev = rnorm;
396: if (it == 0 && ((PetscObject)ksp)->prefix) {
397: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
398: }
399: KSPMonitorRange_Private(ksp,it,&perc);
401: rel = (prev - rnorm)/prev;
402: prev = rnorm;
403: 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));
404: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
405: PetscViewerPopFormat(viewer);
406: return(0);
407: }
411: /*@C
412: KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
413: outer iteration in an adaptive way.
415: Collective on KSP
417: Input Parameters:
418: + ksp - iterative context
419: . n - iteration number (not used)
420: . fnorm - the current residual norm
421: . dummy - some context as a C struct. fields:
422: coef: a scaling coefficient. default 1.0. can be passed through
423: -sub_ksp_dynamic_tolerance_param
424: bnrm: norm of the right-hand side. store it to avoid repeated calculation
426: Notes:
427: This may be useful for a flexibly preconditioner Krylov method to
428: control the accuracy of the inner solves needed to gaurantee the
429: convergence of the outer iterations.
431: Level: advanced
433: .keywords: KSP, inner tolerance
435: .seealso: KSPMonitorDynamicToleranceDestroy()
436: @*/
437: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
438: {
440: PC pc;
441: PetscReal outer_rtol, outer_abstol, outer_dtol, inner_rtol;
442: PetscInt outer_maxits,nksp,first,i;
443: KSPDynTolCtx *scale = (KSPDynTolCtx*)dummy;
444: KSP kspinner = NULL, *subksp = NULL;
447: KSPGetPC(ksp, &pc);
449: /* compute inner_rtol */
450: if (scale->bnrm < 0.0) {
451: Vec b;
452: KSPGetRhs(ksp, &b);
453: VecNorm(b, NORM_2, &(scale->bnrm));
454: }
455: KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
456: inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
457: /*PetscPrintf(PETSC_COMM_WORLD, " Inner rtol = %g\n", (double)inner_rtol);*/
459: /* if pc is ksp */
460: PCKSPGetKSP(pc, &kspinner);
461: if (kspinner) {
462: KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
463: return(0);
464: }
466: /* if pc is bjacobi */
467: PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
468: if (subksp) {
469: for (i=0; i<nksp; i++) {
470: KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
471: }
472: return(0);
473: }
475: /* todo: dynamic tolerance may apply to other types of pc too */
476: return(0);
477: }
481: /*
482: Destroy the dummy context used in KSPMonitorDynamicTolerance()
483: */
484: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
485: {
489: PetscFree(*dummy);
490: return(0);
491: }
495: /*
496: Default (short) KSP Monitor, same as KSPMonitorDefault() except
497: it prints fewer digits of the residual as the residual gets smaller.
498: This is because the later digits are meaningless and are often
499: different on different machines; by using this routine different
500: machines will usually generate the same output.
501: */
502: PetscErrorCode KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *dummy)
503: {
505: PetscViewer viewer = dummy->viewer;
509: PetscViewerPushFormat(viewer,dummy->format);
510: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
511: if (its == 0 && ((PetscObject)ksp)->prefix) {
512: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
513: }
515: if (fnorm > 1.e-9) {
516: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %g \n",its,(double)fnorm);
517: } else if (fnorm > 1.e-11) {
518: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,(double)fnorm);
519: } else {
520: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);
521: }
522: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
523: PetscViewerPopFormat(viewer);
524: return(0);
525: }
529: /*@C
530: KSPConvergedSkip - Convergence test that do not return as converged
531: until the maximum number of iterations is reached.
533: Collective on KSP
535: Input Parameters:
536: + ksp - iterative context
537: . n - iteration number
538: . rnorm - 2-norm residual value (may be estimated)
539: - dummy - unused convergence context
541: Returns:
542: . reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS
544: Notes:
545: This should be used as the convergence test with the option
546: KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
547: not computed. Convergence is then declared after the maximum number
548: of iterations have been reached. Useful when one is using CG or
549: BiCGStab as a smoother.
551: Level: advanced
553: .keywords: KSP, default, convergence, residual
555: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
556: @*/
557: PetscErrorCode KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
558: {
562: *reason = KSP_CONVERGED_ITERATING;
563: if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
564: return(0);
565: }
570: /*@C
571: KSPConvergedDefaultCreate - Creates and initializes the space used by the KSPConvergedDefault() function context
573: Collective on KSP
575: Output Parameter:
576: . ctx - convergence context
578: Level: intermediate
580: .keywords: KSP, default, convergence, residual
582: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
583: KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
584: @*/
585: PetscErrorCode KSPConvergedDefaultCreate(void **ctx)
586: {
587: PetscErrorCode ierr;
588: KSPConvergedDefaultCtx *cctx;
591: PetscNew(&cctx);
592: *ctx = cctx;
593: return(0);
594: }
598: /*@
599: KSPConvergedDefaultSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
600: instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
601: is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.
603: Collective on KSP
605: Input Parameters:
606: . ksp - iterative context
608: Options Database:
609: . -ksp_converged_use_initial_residual_norm
611: Notes:
612: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
614: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
615: are defined in petscksp.h.
617: If the convergence test is not KSPConvergedDefault() then this is ignored.
619: If right preconditioning is being used then B does not appear in the above formula.
622: Level: intermediate
624: .keywords: KSP, default, convergence, residual
626: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm()
627: @*/
628: PetscErrorCode KSPConvergedDefaultSetUIRNorm(KSP ksp)
629: {
630: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
634: if (ksp->converged != KSPConvergedDefault) return(0);
635: if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
636: ctx->initialrtol = PETSC_TRUE;
637: return(0);
638: }
642: /*@
643: KSPConvergedDefaultSetUMIRNorm - makes the default convergence test use min(|| B*(b - A*(initial guess))||,|| B*b ||)
644: In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
645: is used there is no B in the above formula. UMIRNorm is short for Use Minimum Initial Residual Norm.
647: Collective on KSP
649: Input Parameters:
650: . ksp - iterative context
652: Options Database:
653: . -ksp_converged_use_min_initial_residual_norm
655: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
657: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
658: are defined in petscksp.h.
660: Level: intermediate
662: .keywords: KSP, default, convergence, residual
664: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm()
665: @*/
666: PetscErrorCode KSPConvergedDefaultSetUMIRNorm(KSP ksp)
667: {
668: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
672: if (ksp->converged != KSPConvergedDefault) return(0);
673: if (ctx->initialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
674: ctx->mininitialrtol = PETSC_TRUE;
675: return(0);
676: }
680: /*@C
681: KSPConvergedDefault - Determines convergence of the linear iterative solvers by default
683: Collective on KSP
685: Input Parameters:
686: + ksp - iterative context
687: . n - iteration number
688: . rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
689: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
691: Output Parameter:
692: + positive - if the iteration has converged;
693: . negative - if residual norm exceeds divergence threshold;
694: - 0 - otherwise.
696: Notes:
697: KSPConvergedDefault() reaches convergence when rnorm < MAX (rtol * rnorm_0, abstol);
698: Divergence is detected if rnorm > dtol * rnorm_0,
700: where:
701: + rtol = relative tolerance,
702: . abstol = absolute tolerance.
703: . dtol = divergence tolerance,
704: - rnorm_0 is the two norm of the right hand side. When initial guess is non-zero you
705: can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
706: as the starting point for relative norm convergence testing, that is as rnorm_0
708: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
710: Use KSPSetNormType() (or -ksp_norm_type <none,preconditioned,unpreconditioned,natural>) to change the norm used for computing rnorm
712: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
714: This routine is used by KSP by default so the user generally never needs call it directly.
716: Use KSPSetConvergenceTest() to provide your own test instead of using this one.
718: Level: intermediate
720: .keywords: KSP, default, convergence, residual
722: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
723: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
724: @*/
725: PetscErrorCode KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
726: {
727: PetscErrorCode ierr;
728: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
729: KSPNormType normtype;
734: *reason = KSP_CONVERGED_ITERATING;
736: KSPGetNormType(ksp,&normtype);
737: if (normtype == KSP_NORM_NONE) return(0);
739: if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
740: if (!n) {
741: /* if user gives initial guess need to compute norm of b */
742: if (!ksp->guess_zero && !cctx->initialrtol) {
743: PetscReal snorm = 0.0;
744: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
745: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
746: VecNorm(ksp->vec_rhs,NORM_2,&snorm); /* <- b'*b */
747: } else {
748: Vec z;
749: /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
750: VecDuplicate(ksp->vec_rhs,&z);
751: KSP_PCApply(ksp,ksp->vec_rhs,z);
752: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
753: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
754: VecNorm(z,NORM_2,&snorm); /* dp <- b'*B'*B*b */
755: } else if (ksp->normtype == KSP_NORM_NATURAL) {
756: PetscScalar norm;
757: PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
758: VecDot(ksp->vec_rhs,z,&norm);
759: snorm = PetscSqrtReal(PetscAbsScalar(norm)); /* dp <- b'*B*b */
760: }
761: VecDestroy(&z);
762: }
763: /* handle special case of zero RHS and nonzero guess */
764: if (!snorm) {
765: PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
766: snorm = rnorm;
767: }
768: if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
769: else ksp->rnorm0 = snorm;
770: } else {
771: ksp->rnorm0 = rnorm;
772: }
773: ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
774: }
776: if (n <= ksp->chknorm) return(0);
778: if (PetscIsInfOrNanReal(rnorm)) {
779: PCFailedReason pcreason;
780: PetscInt sendbuf,pcreason_max;
781: PCGetSetUpFailedReason(ksp->pc,&pcreason);
782: sendbuf = (PetscInt)pcreason;
783: MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp));
784: if (pcreason_max) {
785: *reason = KSP_DIVERGED_PCSETUP_FAILED;
786: VecSetInf(ksp->vec_sol);
787: PetscInfo(ksp,"Linear solver pcsetup fails, declaring divergence \n");
788: } else {
789: *reason = KSP_DIVERGED_NANORINF;
790: PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
791: }
792: } else if (rnorm <= ksp->ttol) {
793: if (rnorm < ksp->abstol) {
794: 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);
795: *reason = KSP_CONVERGED_ATOL;
796: } else {
797: if (cctx->initialrtol) {
798: 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);
799: } else {
800: 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);
801: }
802: *reason = KSP_CONVERGED_RTOL;
803: }
804: } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
805: 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);
806: *reason = KSP_DIVERGED_DTOL;
807: }
808: return(0);
809: }
813: /*@C
814: KSPConvergedDefaultDestroy - Frees the space used by the KSPConvergedDefault() function context
816: Collective on KSP
818: Input Parameters:
819: . ctx - convergence context
821: Level: intermediate
823: .keywords: KSP, default, convergence, residual
825: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
826: KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
827: @*/
828: PetscErrorCode KSPConvergedDefaultDestroy(void *ctx)
829: {
830: PetscErrorCode ierr;
831: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
834: VecDestroy(&cctx->work);
835: PetscFree(ctx);
836: return(0);
837: }
841: /*
842: KSPBuildSolutionDefault - Default code to create/move the solution.
844: Input Parameters:
845: + ksp - iterative context
846: - v - pointer to the user's vector
848: Output Parameter:
849: . V - pointer to a vector containing the solution
851: Level: advanced
853: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
855: .keywords: KSP, build, solution, default
857: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
858: */
859: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
860: {
864: if (ksp->pc_side == PC_RIGHT) {
865: if (ksp->pc) {
866: if (v) {
867: KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
868: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
869: } else {
870: if (v) {
871: VecCopy(ksp->vec_sol,v); *V = v;
872: } else *V = ksp->vec_sol;
873: }
874: } else if (ksp->pc_side == PC_SYMMETRIC) {
875: if (ksp->pc) {
876: if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
877: if (v) {
878: PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
879: *V = v;
880: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
881: } else {
882: if (v) {
883: VecCopy(ksp->vec_sol,v); *V = v;
884: } else *V = ksp->vec_sol;
885: }
886: } else {
887: if (v) {
888: VecCopy(ksp->vec_sol,v); *V = v;
889: } else *V = ksp->vec_sol;
890: }
891: return(0);
892: }
896: /*
897: KSPBuildResidualDefault - Default code to compute the residual.
899: Input Parameters:
900: . ksp - iterative context
901: . t - pointer to temporary vector
902: . v - pointer to user vector
904: Output Parameter:
905: . V - pointer to a vector containing the residual
907: Level: advanced
909: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
911: .keywords: KSP, build, residual, default
913: .seealso: KSPBuildSolutionDefault()
914: */
915: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
916: {
918: Mat Amat,Pmat;
921: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
922: PCGetOperators(ksp->pc,&Amat,&Pmat);
923: KSPBuildSolution(ksp,t,NULL);
924: KSP_MatMult(ksp,Amat,t,v);
925: VecAYPX(v,-1.0,ksp->vec_rhs);
926: *V = v;
927: return(0);
928: }
932: /*@C
933: KSPCreateVecs - Gets a number of work vectors.
935: Input Parameters:
936: + ksp - iterative context
937: . rightn - number of right work vectors
938: - leftn - number of left work vectors to allocate
940: Output Parameter:
941: + right - the array of vectors created
942: - left - the array of left vectors
944: Note: The right vector has as many elements as the matrix has columns. The left
945: vector has as many elements as the matrix has rows.
947: The vectors are new vectors that are not owned by the KSP, they should be destroyed with calls to VecDestroyVecs() when no longer needed.
949: 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
950: that does not exist tries to get them from the DM (if it is provided).
952: Level: advanced
954: .seealso: MatCreateVecs(), VecDestroyVecs()
956: @*/
957: PetscErrorCode KSPCreateVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
958: {
960: Vec vecr = NULL,vecl = NULL;
961: PetscBool matset,pmatset;
962: Mat mat = NULL;
965: if (rightn) {
966: if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
967: if (ksp->vec_sol) vecr = ksp->vec_sol;
968: else {
969: if (ksp->pc) {
970: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
971: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
972: if (matset) {
973: PCGetOperators(ksp->pc,&mat,NULL);
974: MatCreateVecs(mat,&vecr,NULL);
975: } else if (pmatset) {
976: PCGetOperators(ksp->pc,NULL,&mat);
977: MatCreateVecs(mat,&vecr,NULL);
978: }
979: }
980: if (!vecr) {
981: if (ksp->dm) {
982: DMGetGlobalVector(ksp->dm,&vecr);
983: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
984: }
985: }
986: VecDuplicateVecs(vecr,rightn,right);
987: if (!ksp->vec_sol) {
988: if (mat) {
989: VecDestroy(&vecr);
990: } else if (ksp->dm) {
991: DMRestoreGlobalVector(ksp->dm,&vecr);
992: }
993: }
994: }
995: if (leftn) {
996: if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
997: if (ksp->vec_rhs) vecl = ksp->vec_rhs;
998: else {
999: if (ksp->pc) {
1000: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
1001: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
1002: if (matset) {
1003: PCGetOperators(ksp->pc,&mat,NULL);
1004: MatCreateVecs(mat,NULL,&vecl);
1005: } else if (pmatset) {
1006: PCGetOperators(ksp->pc,NULL,&mat);
1007: MatCreateVecs(mat,NULL,&vecl);
1008: }
1009: }
1010: if (!vecl) {
1011: if (ksp->dm) {
1012: DMGetGlobalVector(ksp->dm,&vecl);
1013: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
1014: }
1015: }
1016: VecDuplicateVecs(vecl,leftn,left);
1017: if (!ksp->vec_rhs) {
1018: if (mat) {
1019: VecDestroy(&vecl);
1020: } else if (ksp->dm) {
1021: DMRestoreGlobalVector(ksp->dm,&vecl);
1022: }
1023: }
1024: }
1025: return(0);
1026: }
1030: /*
1031: KSPSetWorkVecs - Sets a number of work vectors into a KSP object
1033: Input Parameters:
1034: . ksp - iterative context
1035: . nw - number of work vectors to allocate
1037: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1039: */
1040: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1041: {
1045: VecDestroyVecs(ksp->nwork,&ksp->work);
1046: ksp->nwork = nw;
1047: KSPCreateVecs(ksp,nw,&ksp->work,0,NULL);
1048: PetscLogObjectParents(ksp,nw,ksp->work);
1049: return(0);
1050: }
1054: /*
1055: KSPDestroyDefault - Destroys a iterative context variable for methods with
1056: no separate context. Preferred calling sequence KSPDestroy().
1058: Input Parameter:
1059: . ksp - the iterative context
1061: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1063: */
1064: PetscErrorCode KSPDestroyDefault(KSP ksp)
1065: {
1070: PetscFree(ksp->data);
1071: return(0);
1072: }
1076: /*@
1077: KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.
1079: Not Collective
1081: Input Parameter:
1082: . ksp - the KSP context
1084: Output Parameter:
1085: . reason - negative value indicates diverged, positive value converged, see KSPConvergedReason
1087: Possible values for reason: See also manual page for each reason
1088: $ KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1089: $ KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1090: $ KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence test routine is set.
1091: $ KSP_CONVERGED_CG_NEG_CURVE
1092: $ KSP_CONVERGED_CG_CONSTRAINED
1093: $ KSP_CONVERGED_STEP_LENGTH
1094: $ KSP_DIVERGED_ITS (required more than its to reach convergence)
1095: $ KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1096: $ KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1097: $ KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1098: $ KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial residual. Try a different preconditioner, or a different initial Level.)
1100: Notes: Can only be called after the call the KSPSolve() is complete.
1102: Level: intermediate
1104: .keywords: KSP, nonlinear, set, convergence, test
1106: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason
1107: @*/
1108: PetscErrorCode KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1109: {
1113: *reason = ksp->reason;
1114: return(0);
1115: }
1117: #include <petsc/private/dmimpl.h>
1120: /*@
1121: KSPSetDM - Sets the DM that may be used by some preconditioners
1123: Logically Collective on KSP
1125: Input Parameters:
1126: + ksp - the preconditioner context
1127: - dm - the dm
1129: Notes: If this is used then the KSP will attempt to use the DM to create the matrix and use the routine
1130: set with DMKSPSetComputeOperators(). Use KSPSetDMActive(ksp,PETSC_FALSE) to instead use the matrix
1131: you've provided with KSPSetOperators().
1133: Level: intermediate
1135: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1136: @*/
1137: PetscErrorCode KSPSetDM(KSP ksp,DM dm)
1138: {
1140: PC pc;
1144: if (dm) {PetscObjectReference((PetscObject)dm);}
1145: if (ksp->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
1146: if (ksp->dm->dmksp && ksp->dmAuto && !dm->dmksp) {
1147: DMKSP kdm;
1148: DMCopyDMKSP(ksp->dm,dm);
1149: DMGetDMKSP(ksp->dm,&kdm);
1150: if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1151: }
1152: DMDestroy(&ksp->dm);
1153: }
1154: ksp->dm = dm;
1155: ksp->dmAuto = PETSC_FALSE;
1156: KSPGetPC(ksp,&pc);
1157: PCSetDM(pc,dm);
1158: ksp->dmActive = PETSC_TRUE;
1159: return(0);
1160: }
1164: /*@
1165: KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side
1167: Logically Collective on KSP
1169: Input Parameters:
1170: + ksp - the preconditioner context
1171: - flg - use the DM
1173: Notes:
1174: 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.
1176: Level: intermediate
1178: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1179: @*/
1180: PetscErrorCode KSPSetDMActive(KSP ksp,PetscBool flg)
1181: {
1185: ksp->dmActive = flg;
1186: return(0);
1187: }
1191: /*@
1192: KSPGetDM - Gets the DM that may be used by some preconditioners
1194: Not Collective
1196: Input Parameter:
1197: . ksp - the preconditioner context
1199: Output Parameter:
1200: . dm - the dm
1202: Level: intermediate
1205: .seealso: KSPSetDM(), KSPSetDMActive()
1206: @*/
1207: PetscErrorCode KSPGetDM(KSP ksp,DM *dm)
1208: {
1213: if (!ksp->dm) {
1214: DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1215: ksp->dmAuto = PETSC_TRUE;
1216: }
1217: *dm = ksp->dm;
1218: return(0);
1219: }
1223: /*@
1224: KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.
1226: Logically Collective on KSP
1228: Input Parameters:
1229: + ksp - the KSP context
1230: - usrP - optional user context
1232: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
1233: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1235: Level: intermediate
1237: .keywords: KSP, set, application, context
1239: .seealso: KSPGetApplicationContext()
1240: @*/
1241: PetscErrorCode KSPSetApplicationContext(KSP ksp,void *usrP)
1242: {
1244: PC pc;
1248: ksp->user = usrP;
1249: KSPGetPC(ksp,&pc);
1250: PCSetApplicationContext(pc,usrP);
1251: return(0);
1252: }
1256: /*@
1257: KSPGetApplicationContext - Gets the user-defined context for the linear solver.
1259: Not Collective
1261: Input Parameter:
1262: . ksp - KSP context
1264: Output Parameter:
1265: . usrP - user context
1267: Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
1268: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1270: Level: intermediate
1272: .keywords: KSP, get, application, context
1274: .seealso: KSPSetApplicationContext()
1275: @*/
1276: PetscErrorCode KSPGetApplicationContext(KSP ksp,void *usrP)
1277: {
1280: *(void**)usrP = ksp->user;
1281: return(0);
1282: }