Actual source code: iterativ.c
petsc-3.5.4 2015-05-23
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()
62: @*/
63: PetscErrorCode KSPGetIterationNumber(KSP ksp,PetscInt *its)
64: {
68: *its = ksp->its;
69: return(0);
70: }
74: /*@C
75: KSPMonitorSingularValue - Prints the two norm of the true residual and
76: estimation of the extreme singular values of the preconditioned problem
77: at each iteration.
79: Logically Collective on KSP
81: Input Parameters:
82: + ksp - the iterative context
83: . n - the iteration
84: - rnorm - the two norm of the residual
86: Options Database Key:
87: . -ksp_monitor_singular_value - Activates KSPMonitorSingularValue()
89: Notes:
90: The CG solver uses the Lanczos technique for eigenvalue computation,
91: while GMRES uses the Arnoldi technique; other iterative methods do
92: not currently compute singular values.
94: Level: intermediate
96: .keywords: KSP, CG, default, monitor, extreme, singular values, Lanczos, Arnoldi
98: .seealso: KSPComputeExtremeSingularValues()
99: @*/
100: PetscErrorCode KSPMonitorSingularValue(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
101: {
102: PetscReal emin,emax,c;
104: PetscViewer viewer = (PetscViewer) dummy;
108: if (!viewer) {
109: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
110: }
111: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
112: if (!ksp->calc_sings) {
113: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
114: } else {
115: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
116: c = emax/emin;
117: 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);
118: }
119: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
120: return(0);
121: }
125: /*@C
126: KSPMonitorSolution - Monitors progress of the KSP solvers by calling
127: VecView() for the approximate solution at each iteration.
129: Collective on KSP
131: Input Parameters:
132: + ksp - the KSP context
133: . its - iteration number
134: . fgnorm - 2-norm of residual (or gradient)
135: - dummy - either a viewer or NULL
137: Level: intermediate
139: Notes:
140: For some Krylov methods such as GMRES constructing the solution at
141: each iteration is expensive, hence using this will slow the code.
143: .keywords: KSP, nonlinear, vector, monitor, view
145: .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView()
146: @*/
147: PetscErrorCode KSPMonitorSolution(KSP ksp,PetscInt its,PetscReal fgnorm,void *dummy)
148: {
150: Vec x;
151: PetscViewer viewer = (PetscViewer) dummy;
154: KSPBuildSolution(ksp,NULL,&x);
155: if (!viewer) {
156: MPI_Comm comm;
157: PetscObjectGetComm((PetscObject)ksp,&comm);
158: viewer = PETSC_VIEWER_DRAW_(comm);
159: }
160: VecView(x,viewer);
161: return(0);
162: }
166: /*@C
167: KSPMonitorDefault - Print the residual norm at each iteration of an
168: iterative solver.
170: Collective on KSP
172: Input Parameters:
173: + ksp - iterative context
174: . n - iteration number
175: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
176: - dummy - unused monitor context
178: Level: intermediate
180: .keywords: KSP, default, monitor, residual
182: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
183: @*/
184: PetscErrorCode KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
185: {
187: PetscViewer viewer = (PetscViewer) dummy;
190: if (!viewer) {
191: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
192: }
193: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
194: if (n == 0 && ((PetscObject)ksp)->prefix) {
195: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
196: }
197: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
198: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
199: return(0);
200: }
204: /*@C
205: KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
206: residual norm at each iteration of an iterative solver.
208: Collective on KSP
210: Input Parameters:
211: + ksp - iterative context
212: . n - iteration number
213: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
214: - dummy - unused monitor context
216: Options Database Key:
217: . -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()
219: Notes:
220: When using right preconditioning, these values are equivalent.
222: Level: intermediate
224: .keywords: KSP, default, monitor, residual
226: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
227: @*/
228: PetscErrorCode KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
229: {
231: Vec resid;
232: PetscReal truenorm,bnorm;
233: PetscViewer viewer = (PetscViewer)dummy;
234: char normtype[256];
237: if (!viewer) {
238: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
239: }
240: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
241: if (n == 0 && ((PetscObject)ksp)->prefix) {
242: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
243: }
244: KSPBuildResidual(ksp,NULL,NULL,&resid);
245: VecNorm(resid,NORM_2,&truenorm);
246: VecDestroy(&resid);
247: VecNorm(ksp->vec_rhs,NORM_2,&bnorm);
248: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
249: PetscStrtolower(normtype);
250: 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));
251: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
252: return(0);
253: }
257: /*@C
258: KSPMonitorTrueResidualMaxNorm - Prints the true residual max norm as well as the preconditioned
259: residual norm at each iteration of an iterative solver.
261: Collective on KSP
263: Input Parameters:
264: + ksp - iterative context
265: . n - iteration number
266: . rnorm - norm (preconditioned) residual value (may be estimated).
267: - dummy - unused monitor context
269: Options Database Key:
270: . -ksp_monitor_max - Activates KSPMonitorTrueResidualMaxNorm()
272: Notes:
273: This could be implemented (better) with a flag in ksp.
275: Level: intermediate
277: .keywords: KSP, default, monitor, residual
279: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualNorm()
280: @*/
281: PetscErrorCode KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
282: {
284: Vec resid;
285: PetscReal truenorm,bnorm;
286: PetscViewer viewer = (PetscViewer)dummy;
287: char normtype[256];
290: if (!viewer) {
291: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
292: }
293: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
294: if (n == 0 && ((PetscObject)ksp)->prefix) {
295: PetscViewerASCIIPrintf(viewer," Residual norms (max) for %s solve.\n",((PetscObject)ksp)->prefix);
296: }
297: KSPBuildResidual(ksp,NULL,NULL,&resid);
298: VecNorm(resid,NORM_INFINITY,&truenorm);
299: VecDestroy(&resid);
300: VecNorm(ksp->vec_rhs,NORM_INFINITY,&bnorm);
301: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
302: PetscStrtolower(normtype);
303: /* PetscViewerASCIIPrintf(viewer,"%3D KSP %s resid norm %14.12e true resid norm %14.12e ||r(i)||_inf/||b||_inf %14.12e\n",n,normtype,(double)rnorm,(double)truenorm,(double)(truenorm/bnorm)); */
304: PetscViewerASCIIPrintf(viewer,"%3D KSP true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n",n,(double)truenorm,(double)(truenorm/bnorm));
305: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
306: return(0);
307: }
311: PetscErrorCode KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
312: {
314: Vec resid;
315: PetscReal rmax,pwork;
316: PetscInt i,n,N;
317: const PetscScalar *r;
320: KSPBuildResidual(ksp,NULL,NULL,&resid);
321: VecNorm(resid,NORM_INFINITY,&rmax);
322: VecGetLocalSize(resid,&n);
323: VecGetSize(resid,&N);
324: VecGetArrayRead(resid,&r);
325: pwork = 0.0;
326: for (i=0; i<n; i++) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
327: VecRestoreArrayRead(resid,&r);
328: VecDestroy(&resid);
329: MPI_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
330: *per = *per/N;
331: return(0);
332: }
336: /*@C
337: KSPMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.
339: Collective on KSP
341: Input Parameters:
342: + ksp - iterative context
343: . it - iteration number
344: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
345: - dummy - unused monitor context
347: Options Database Key:
348: . -ksp_monitor_range - Activates KSPMonitorRange()
350: Level: intermediate
352: .keywords: KSP, default, monitor, residual
354: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
355: @*/
356: PetscErrorCode KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,void *dummy)
357: {
358: PetscErrorCode ierr;
359: PetscReal perc,rel;
360: PetscViewer viewer = (PetscViewer)dummy;
361: /* should be in a MonitorRangeContext */
362: static PetscReal prev;
365: if (!viewer) {
366: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
367: }
368: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
369: if (!it) prev = rnorm;
370: if (it == 0 && ((PetscObject)ksp)->prefix) {
371: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
372: }
373: KSPMonitorRange_Private(ksp,it,&perc);
375: rel = (prev - rnorm)/prev;
376: prev = rnorm;
377: 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));
378: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
379: return(0);
380: }
384: /*@C
385: KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
386: outer iteration in an adaptive way.
388: Collective on KSP
390: Input Parameters:
391: + ksp - iterative context
392: . n - iteration number (not used)
393: . fnorm - the current residual norm
394: . dummy - some context as a C struct. fields:
395: coef: a scaling coefficient. default 1.0. can be passed through
396: -sub_ksp_dynamic_tolerance_param
397: bnrm: norm of the right-hand side. store it to avoid repeated calculation
399: Notes:
400: This may be useful for a flexibly preconditioner Krylov method to
401: control the accuracy of the inner solves needed to gaurantee the
402: convergence of the outer iterations.
404: Level: advanced
406: .keywords: KSP, inner tolerance
408: .seealso: KSPMonitorDynamicToleranceDestroy()
409: @*/
410: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
411: {
413: PC pc;
414: PetscReal outer_rtol, outer_abstol, outer_dtol, inner_rtol;
415: PetscInt outer_maxits,nksp,first,i;
416: KSPDynTolCtx *scale = (KSPDynTolCtx*)dummy;
417: KSP kspinner = NULL, *subksp = NULL;
420: KSPGetPC(ksp, &pc);
422: /* compute inner_rtol */
423: if (scale->bnrm < 0.0) {
424: Vec b;
425: KSPGetRhs(ksp, &b);
426: VecNorm(b, NORM_2, &(scale->bnrm));
427: }
428: KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
429: inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
430: /*PetscPrintf(PETSC_COMM_WORLD, " Inner rtol = %g\n", (double)inner_rtol);*/
432: /* if pc is ksp */
433: PCKSPGetKSP(pc, &kspinner);
434: if (kspinner) {
435: KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
436: return(0);
437: }
439: /* if pc is bjacobi */
440: PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
441: if (subksp) {
442: for (i=0; i<nksp; i++) {
443: KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
444: }
445: return(0);
446: }
448: /* todo: dynamic tolerance may apply to other types of pc too */
449: return(0);
450: }
454: /*
455: Destroy the dummy context used in KSPMonitorDynamicTolerance()
456: */
457: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
458: {
462: PetscFree(*dummy);
463: return(0);
464: }
468: /*
469: Default (short) KSP Monitor, same as KSPMonitorDefault() except
470: it prints fewer digits of the residual as the residual gets smaller.
471: This is because the later digits are meaningless and are often
472: different on different machines; by using this routine different
473: machines will usually generate the same output.
474: */
475: PetscErrorCode KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
476: {
478: PetscViewer viewer = (PetscViewer)dummy;
481: if (!viewer) {
482: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
483: }
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: return(0);
498: }
502: /*@C
503: KSPConvergedSkip - Convergence test that do not return as converged
504: until the maximum number of iterations is reached.
506: Collective on KSP
508: Input Parameters:
509: + ksp - iterative context
510: . n - iteration number
511: . rnorm - 2-norm residual value (may be estimated)
512: - dummy - unused convergence context
514: Returns:
515: . reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS
517: Notes:
518: This should be used as the convergence test with the option
519: KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
520: not computed. Convergence is then declared after the maximum number
521: of iterations have been reached. Useful when one is using CG or
522: BiCGStab as a smoother.
524: Level: advanced
526: .keywords: KSP, default, convergence, residual
528: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
529: @*/
530: PetscErrorCode KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
531: {
535: *reason = KSP_CONVERGED_ITERATING;
536: if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
537: return(0);
538: }
543: /*@C
544: KSPConvergedDefaultCreate - Creates and initializes the space used by the KSPConvergedDefault() function context
546: Collective on KSP
548: Output Parameter:
549: . ctx - convergence context
551: Level: intermediate
553: .keywords: KSP, default, convergence, residual
555: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
556: KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
557: @*/
558: PetscErrorCode KSPConvergedDefaultCreate(void **ctx)
559: {
560: PetscErrorCode ierr;
561: KSPConvergedDefaultCtx *cctx;
564: PetscNew(&cctx);
565: *ctx = cctx;
566: return(0);
567: }
571: /*@
572: KSPConvergedDefaultSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
573: instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
574: is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.
576: Collective on KSP
578: Input Parameters:
579: . ksp - iterative context
581: Options Database:
582: . -ksp_converged_use_initial_residual_norm
584: Notes:
585: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
587: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
588: are defined in petscksp.h.
590: If the convergence test is not KSPConvergedDefault() then this is ignored.
592: If right preconditioning is being used then B does not appear in the above formula.
595: Level: intermediate
597: .keywords: KSP, default, convergence, residual
599: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm()
600: @*/
601: PetscErrorCode KSPConvergedDefaultSetUIRNorm(KSP ksp)
602: {
603: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
607: if (ksp->converged != KSPConvergedDefault) return(0);
608: if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
609: ctx->initialrtol = PETSC_TRUE;
610: return(0);
611: }
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: }
653: /*@C
654: KSPConvergedDefault - Determines convergence of the linear iterative solvers by default
656: Collective on KSP
658: Input Parameters:
659: + ksp - iterative context
660: . n - iteration number
661: . rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
662: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
664: Output Parameter:
665: + positive - if the iteration has converged;
666: . negative - if residual norm exceeds divergence threshold;
667: - 0 - otherwise.
669: Notes:
670: KSPConvergedDefault() reaches convergence when rnorm < MAX (rtol * rnorm_0, abstol);
671: Divergence is detected if rnorm > dtol * rnorm_0,
673: where:
674: + rtol = relative tolerance,
675: . abstol = absolute tolerance.
676: . dtol = divergence tolerance,
677: - rnorm_0 is the two norm of the right hand side. When initial guess is non-zero you
678: can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
679: as the starting point for relative norm convergence testing, that is as rnorm_0
681: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
683: Use KSPSetNormType() (or -ksp_norm_type <none,preconditioned,unpreconditioned,natural>) to change the norm used for computing rnorm
685: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
687: This routine is used by KSP by default so the user generally never needs call it directly.
689: Use KSPSetConvergenceTest() to provide your own test instead of using this one.
691: Level: intermediate
693: .keywords: KSP, default, convergence, residual
695: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
696: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
697: @*/
698: PetscErrorCode KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
699: {
700: PetscErrorCode ierr;
701: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
702: KSPNormType normtype;
707: *reason = KSP_CONVERGED_ITERATING;
709: KSPGetNormType(ksp,&normtype);
710: if (normtype == KSP_NORM_NONE) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Use KSPConvergedSkip() with KSPNormType of KSP_NORM_NONE");
712: if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
713: if (!n) {
714: /* if user gives initial guess need to compute norm of b */
715: if (!ksp->guess_zero && !cctx->initialrtol) {
716: PetscReal snorm;
717: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
718: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
719: VecNorm(ksp->vec_rhs,NORM_2,&snorm); /* <- b'*b */
720: } else {
721: Vec z;
722: /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
723: VecDuplicate(ksp->vec_rhs,&z);
724: KSP_PCApply(ksp,ksp->vec_rhs,z);
725: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
726: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
727: VecNorm(z,NORM_2,&snorm); /* dp <- b'*B'*B*b */
728: } else if (ksp->normtype == KSP_NORM_NATURAL) {
729: PetscScalar norm;
730: PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
731: VecDot(ksp->vec_rhs,z,&norm);
732: snorm = PetscSqrtReal(PetscAbsScalar(norm)); /* dp <- b'*B*b */
733: }
734: VecDestroy(&z);
735: }
736: /* handle special case of zero RHS and nonzero guess */
737: if (!snorm) {
738: PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
739: snorm = rnorm;
740: }
741: if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
742: else ksp->rnorm0 = snorm;
743: } else {
744: ksp->rnorm0 = rnorm;
745: }
746: ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
747: }
749: if (n <= ksp->chknorm) return(0);
751: if (PetscIsInfOrNanScalar(rnorm)) {
752: PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
753: *reason = KSP_DIVERGED_NANORINF;
754: } else if (rnorm <= ksp->ttol) {
755: if (rnorm < ksp->abstol) {
756: 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);
757: *reason = KSP_CONVERGED_ATOL;
758: } else {
759: if (cctx->initialrtol) {
760: 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);
761: } else {
762: 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);
763: }
764: *reason = KSP_CONVERGED_RTOL;
765: }
766: } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
767: 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);
768: *reason = KSP_DIVERGED_DTOL;
769: }
770: return(0);
771: }
775: /*@C
776: KSPConvergedDefaultDestroy - Frees the space used by the KSPConvergedDefault() function context
778: Collective on KSP
780: Input Parameters:
781: . ctx - convergence context
783: Level: intermediate
785: .keywords: KSP, default, convergence, residual
787: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
788: KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
789: @*/
790: PetscErrorCode KSPConvergedDefaultDestroy(void *ctx)
791: {
792: PetscErrorCode ierr;
793: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
796: VecDestroy(&cctx->work);
797: PetscFree(ctx);
798: return(0);
799: }
803: /*
804: KSPBuildSolutionDefault - Default code to create/move the solution.
806: Input Parameters:
807: + ksp - iterative context
808: - v - pointer to the user's vector
810: Output Parameter:
811: . V - pointer to a vector containing the solution
813: Level: advanced
815: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
817: .keywords: KSP, build, solution, default
819: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
820: */
821: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
822: {
826: if (ksp->pc_side == PC_RIGHT) {
827: if (ksp->pc) {
828: if (v) {
829: KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
830: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
831: } else {
832: if (v) {
833: VecCopy(ksp->vec_sol,v); *V = v;
834: } else *V = ksp->vec_sol;
835: }
836: } else if (ksp->pc_side == PC_SYMMETRIC) {
837: if (ksp->pc) {
838: if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
839: if (v) {
840: PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
841: *V = v;
842: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
843: } else {
844: if (v) {
845: VecCopy(ksp->vec_sol,v); *V = v;
846: } else *V = ksp->vec_sol;
847: }
848: } else {
849: if (v) {
850: VecCopy(ksp->vec_sol,v); *V = v;
851: } else *V = ksp->vec_sol;
852: }
853: return(0);
854: }
858: /*
859: KSPBuildResidualDefault - Default code to compute the residual.
861: Input Parameters:
862: . ksp - iterative context
863: . t - pointer to temporary vector
864: . v - pointer to user vector
866: Output Parameter:
867: . V - pointer to a vector containing the residual
869: Level: advanced
871: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
873: .keywords: KSP, build, residual, default
875: .seealso: KSPBuildSolutionDefault()
876: */
877: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
878: {
880: Mat Amat,Pmat;
883: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
884: PCGetOperators(ksp->pc,&Amat,&Pmat);
885: KSPBuildSolution(ksp,t,NULL);
886: KSP_MatMult(ksp,Amat,t,v);
887: VecAYPX(v,-1.0,ksp->vec_rhs);
888: *V = v;
889: return(0);
890: }
894: /*@C
895: KSPGetVecs - Gets a number of work vectors.
897: Input Parameters:
898: + ksp - iterative context
899: . rightn - number of right work vectors
900: - leftn - number of left work vectors to allocate
902: Output Parameter:
903: + right - the array of vectors created
904: - left - the array of left vectors
906: Note: The right vector has as many elements as the matrix has columns. The left
907: vector has as many elements as the matrix has rows.
909: Level: advanced
911: .seealso: MatGetVecs()
913: @*/
914: PetscErrorCode KSPGetVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
915: {
917: Vec vecr,vecl;
920: if (rightn) {
921: if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
922: if (ksp->vec_sol) vecr = ksp->vec_sol;
923: else {
924: 0;
925: if (ksp->dm) {
926: DMGetGlobalVector(ksp->dm,&vecr); /* don't check for errors -- if any errors, pass down to next block */
927: }
928: if (!ksp->dm || ierr) {
929: Mat mat;
930: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
931: PCGetOperators(ksp->pc,&mat,NULL);
932: MatGetVecs(mat,&vecr,NULL);
933: }
934: }
935: VecDuplicateVecs(vecr,rightn,right);
936: if (!ksp->vec_sol) {
937: if (ksp->dm) {
938: DMRestoreGlobalVector(ksp->dm,&vecr);
939: } else {
940: VecDestroy(&vecr);
941: }
942: }
943: }
944: if (leftn) {
945: if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
946: if (ksp->vec_rhs) vecl = ksp->vec_rhs;
947: else {
948: 0;
949: if (ksp->dm) {
950: DMGetGlobalVector(ksp->dm,&vecl); /* don't check for errors -- if any errors, pass down to next block */
951: }
952: if (!ksp->dm || ierr) {
953: Mat mat;
954: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
955: PCGetOperators(ksp->pc,&mat,NULL);
956: MatGetVecs(mat,NULL,&vecl);
957: }
958: }
959: VecDuplicateVecs(vecl,leftn,left);
960: if (!ksp->vec_rhs) {
961: if (ksp->dm) {
962: DMRestoreGlobalVector(ksp->dm,&vecl);
963: } else {
964: VecDestroy(&vecl);
965: }
966: }
967: }
968: return(0);
969: }
973: /*
974: KSPSetWorkVecs - Sets a number of work vectors into a KSP object
976: Input Parameters:
977: . ksp - iterative context
978: . nw - number of work vectors to allocate
980: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
982: */
983: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
984: {
988: VecDestroyVecs(ksp->nwork,&ksp->work);
989: ksp->nwork = nw;
990: KSPGetVecs(ksp,nw,&ksp->work,0,NULL);
991: PetscLogObjectParents(ksp,nw,ksp->work);
992: return(0);
993: }
997: /*
998: KSPDestroyDefault - Destroys a iterative context variable for methods with
999: no separate context. Preferred calling sequence KSPDestroy().
1001: Input Parameter:
1002: . ksp - the iterative context
1004: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1006: */
1007: PetscErrorCode KSPDestroyDefault(KSP ksp)
1008: {
1013: PetscFree(ksp->data);
1014: return(0);
1015: }
1019: /*@
1020: KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.
1022: Not Collective
1024: Input Parameter:
1025: . ksp - the KSP context
1027: Output Parameter:
1028: . reason - negative value indicates diverged, positive value converged, see KSPConvergedReason
1030: Possible values for reason:
1031: + KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1032: . KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1033: . KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence
1034: test routine is set.
1035: . KSP_CONVERGED_CG_NEG_CURVE
1036: . KSP_CONVERGED_CG_CONSTRAINED
1037: . KSP_CONVERGED_STEP_LENGTH
1038: . KSP_DIVERGED_ITS (required more than its to reach convergence)
1039: . KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1040: . KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1041: . KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1042: - KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial
1043: residual. Try a different preconditioner, or a different initial Level.)
1045: See also manual page for each reason.
1047: guess: beginner
1049: Notes: Can only be called after the call the KSPSolve() is complete.
1051: Level: intermediate
1053: .keywords: KSP, nonlinear, set, convergence, test
1055: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason
1056: @*/
1057: PetscErrorCode KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1058: {
1062: *reason = ksp->reason;
1063: return(0);
1064: }
1066: #include <petsc-private/dmimpl.h>
1069: /*@
1070: KSPSetDM - Sets the DM that may be used by some preconditioners
1072: Logically Collective on KSP
1074: Input Parameters:
1075: + ksp - the preconditioner context
1076: - dm - the dm
1078: Notes: If this is used then the KSP will attempt to use the DM to create the matrix and use the routine
1079: set with DMKSPSetComputeOperators(). Use KSPSetDMActive(ksp,PETSC_FALSE) to instead use the matrix
1080: you've provided with KSPSetOperators().
1082: Level: intermediate
1084: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1085: @*/
1086: PetscErrorCode KSPSetDM(KSP ksp,DM dm)
1087: {
1089: PC pc;
1093: if (dm) {PetscObjectReference((PetscObject)dm);}
1094: if (ksp->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
1095: if (ksp->dm->dmksp && ksp->dmAuto && !dm->dmksp) {
1096: DMKSP kdm;
1097: DMCopyDMKSP(ksp->dm,dm);
1098: DMGetDMKSP(ksp->dm,&kdm);
1099: if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1100: }
1101: DMDestroy(&ksp->dm);
1102: }
1103: ksp->dm = dm;
1104: ksp->dmAuto = PETSC_FALSE;
1105: KSPGetPC(ksp,&pc);
1106: PCSetDM(pc,dm);
1107: ksp->dmActive = PETSC_TRUE;
1108: return(0);
1109: }
1113: /*@
1114: KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side
1116: Logically Collective on KSP
1118: Input Parameters:
1119: + ksp - the preconditioner context
1120: - flg - use the DM
1122: Notes:
1123: 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.
1125: Level: intermediate
1127: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1128: @*/
1129: PetscErrorCode KSPSetDMActive(KSP ksp,PetscBool flg)
1130: {
1134: ksp->dmActive = flg;
1135: return(0);
1136: }
1140: /*@
1141: KSPGetDM - Gets the DM that may be used by some preconditioners
1143: Not Collective
1145: Input Parameter:
1146: . ksp - the preconditioner context
1148: Output Parameter:
1149: . dm - the dm
1151: Level: intermediate
1154: .seealso: KSPSetDM(), KSPSetDMActive()
1155: @*/
1156: PetscErrorCode KSPGetDM(KSP ksp,DM *dm)
1157: {
1162: if (!ksp->dm) {
1163: DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1164: ksp->dmAuto = PETSC_TRUE;
1165: }
1166: *dm = ksp->dm;
1167: return(0);
1168: }
1172: /*@
1173: KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.
1175: Logically Collective on KSP
1177: Input Parameters:
1178: + ksp - the KSP context
1179: - usrP - optional user context
1181: Level: intermediate
1183: .keywords: KSP, set, application, context
1185: .seealso: KSPGetApplicationContext()
1186: @*/
1187: PetscErrorCode KSPSetApplicationContext(KSP ksp,void *usrP)
1188: {
1190: PC pc;
1194: ksp->user = usrP;
1195: KSPGetPC(ksp,&pc);
1196: PCSetApplicationContext(pc,usrP);
1197: return(0);
1198: }
1202: /*@
1203: KSPGetApplicationContext - Gets the user-defined context for the linear solver.
1205: Not Collective
1207: Input Parameter:
1208: . ksp - KSP context
1210: Output Parameter:
1211: . usrP - user context
1213: Level: intermediate
1215: .keywords: KSP, get, application, context
1217: .seealso: KSPSetApplicationContext()
1218: @*/
1219: PetscErrorCode KSPGetApplicationContext(KSP ksp,void *usrP)
1220: {
1223: *(void**)usrP = ksp->user;
1224: return(0);
1225: }