Actual source code: iterativ.c
petsc-3.6.4 2016-04-12
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,void *dummy)
132: {
133: PetscReal emin,emax,c;
135: PetscViewer viewer = (PetscViewer) dummy;
139: if (!viewer) {
140: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
141: }
142: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
143: if (!ksp->calc_sings) {
144: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
145: } else {
146: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
147: c = emax/emin;
148: 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);
149: }
150: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
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 - either a viewer or NULL
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,void *dummy)
179: {
181: Vec x;
182: PetscViewer viewer = (PetscViewer) dummy;
185: KSPBuildSolution(ksp,NULL,&x);
186: if (!viewer) {
187: MPI_Comm comm;
188: PetscObjectGetComm((PetscObject)ksp,&comm);
189: viewer = PETSC_VIEWER_DRAW_(comm);
190: }
191: VecView(x,viewer);
192: return(0);
193: }
195: #if defined(PETSC_HAVE_THREADSAFETY)
196: #define KSPMonitorDefault KSPMonitorDefaultUnsafe
197: #endif
201: /*@C
202: KSPMonitorDefault - Print the residual norm at each iteration of an
203: iterative solver.
205: Collective on KSP
207: Input Parameters:
208: + ksp - iterative context
209: . n - iteration number
210: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
211: - dummy - unused monitor context
213: Level: intermediate
215: .keywords: KSP, default, monitor, residual
217: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
218: @*/
219: PetscErrorCode KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
220: {
222: PetscViewer viewer = (PetscViewer) dummy;
225: if (!viewer) {
226: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
227: }
228: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
229: if (n == 0 && ((PetscObject)ksp)->prefix) {
230: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
231: }
232: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
233: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
234: return(0);
235: }
237: #if defined(PETSC_HAVE_THREADSAFETY)
238: #undef KSPMonitorDefault
239: PetscErrorCode KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
240: {
242: #pragma omp critical
243: KSPMonitorDefaultUnsafe(ksp,n,rnorm,dummy);
244: return ierr;
245: }
246: #endif
250: /*@C
251: KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
252: residual norm at each iteration of an iterative solver.
254: Collective on KSP
256: Input Parameters:
257: + ksp - iterative context
258: . n - iteration number
259: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
260: - dummy - unused monitor context
262: Options Database Key:
263: . -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()
265: Notes:
266: When using right preconditioning, these values are equivalent.
268: Level: intermediate
270: .keywords: KSP, default, monitor, residual
272: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
273: @*/
274: PetscErrorCode KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
275: {
277: Vec resid;
278: PetscReal truenorm,bnorm;
279: PetscViewer viewer = (PetscViewer)dummy;
280: char normtype[256];
283: if (!viewer) {
284: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
285: }
286: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
287: if (n == 0 && ((PetscObject)ksp)->prefix) {
288: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
289: }
290: KSPBuildResidual(ksp,NULL,NULL,&resid);
291: VecNorm(resid,NORM_2,&truenorm);
292: VecDestroy(&resid);
293: VecNorm(ksp->vec_rhs,NORM_2,&bnorm);
294: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
295: PetscStrtolower(normtype);
296: 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));
297: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
298: return(0);
299: }
303: /*@C
304: KSPMonitorTrueResidualMaxNorm - Prints the true residual max norm as well as the preconditioned
305: residual norm at each iteration of an iterative solver.
307: Collective on KSP
309: Input Parameters:
310: + ksp - iterative context
311: . n - iteration number
312: . rnorm - norm (preconditioned) residual value (may be estimated).
313: - dummy - unused monitor context
315: Options Database Key:
316: . -ksp_monitor_max - Activates KSPMonitorTrueResidualMaxNorm()
318: Notes:
319: This could be implemented (better) with a flag in ksp.
321: Level: intermediate
323: .keywords: KSP, default, monitor, residual
325: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualNorm()
326: @*/
327: PetscErrorCode KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
328: {
330: Vec resid;
331: PetscReal truenorm,bnorm;
332: PetscViewer viewer = (PetscViewer)dummy;
333: char normtype[256];
336: if (!viewer) {
337: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
338: }
339: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
340: if (n == 0 && ((PetscObject)ksp)->prefix) {
341: PetscViewerASCIIPrintf(viewer," Residual norms (max) for %s solve.\n",((PetscObject)ksp)->prefix);
342: }
343: KSPBuildResidual(ksp,NULL,NULL,&resid);
344: VecNorm(resid,NORM_INFINITY,&truenorm);
345: VecDestroy(&resid);
346: VecNorm(ksp->vec_rhs,NORM_INFINITY,&bnorm);
347: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
348: PetscStrtolower(normtype);
349: /* 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)); */
350: PetscViewerASCIIPrintf(viewer,"%3D KSP true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n",n,(double)truenorm,(double)(truenorm/bnorm));
351: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
352: return(0);
353: }
357: PetscErrorCode KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
358: {
360: Vec resid;
361: PetscReal rmax,pwork;
362: PetscInt i,n,N;
363: const PetscScalar *r;
366: KSPBuildResidual(ksp,NULL,NULL,&resid);
367: VecNorm(resid,NORM_INFINITY,&rmax);
368: VecGetLocalSize(resid,&n);
369: VecGetSize(resid,&N);
370: VecGetArrayRead(resid,&r);
371: pwork = 0.0;
372: for (i=0; i<n; i++) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
373: VecRestoreArrayRead(resid,&r);
374: VecDestroy(&resid);
375: MPI_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
376: *per = *per/N;
377: return(0);
378: }
382: /*@C
383: KSPMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.
385: Collective on KSP
387: Input Parameters:
388: + ksp - iterative context
389: . it - iteration number
390: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
391: - dummy - unused monitor context
393: Options Database Key:
394: . -ksp_monitor_range - Activates KSPMonitorRange()
396: Level: intermediate
398: .keywords: KSP, default, monitor, residual
400: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
401: @*/
402: PetscErrorCode KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,void *dummy)
403: {
404: PetscErrorCode ierr;
405: PetscReal perc,rel;
406: PetscViewer viewer = (PetscViewer)dummy;
407: /* should be in a MonitorRangeContext */
408: static PetscReal prev;
411: if (!viewer) {
412: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
413: }
414: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
415: if (!it) prev = rnorm;
416: if (it == 0 && ((PetscObject)ksp)->prefix) {
417: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
418: }
419: KSPMonitorRange_Private(ksp,it,&perc);
421: rel = (prev - rnorm)/prev;
422: prev = rnorm;
423: 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));
424: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
425: return(0);
426: }
430: /*@C
431: KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
432: outer iteration in an adaptive way.
434: Collective on KSP
436: Input Parameters:
437: + ksp - iterative context
438: . n - iteration number (not used)
439: . fnorm - the current residual norm
440: . dummy - some context as a C struct. fields:
441: coef: a scaling coefficient. default 1.0. can be passed through
442: -sub_ksp_dynamic_tolerance_param
443: bnrm: norm of the right-hand side. store it to avoid repeated calculation
445: Notes:
446: This may be useful for a flexibly preconditioner Krylov method to
447: control the accuracy of the inner solves needed to gaurantee the
448: convergence of the outer iterations.
450: Level: advanced
452: .keywords: KSP, inner tolerance
454: .seealso: KSPMonitorDynamicToleranceDestroy()
455: @*/
456: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
457: {
459: PC pc;
460: PetscReal outer_rtol, outer_abstol, outer_dtol, inner_rtol;
461: PetscInt outer_maxits,nksp,first,i;
462: KSPDynTolCtx *scale = (KSPDynTolCtx*)dummy;
463: KSP kspinner = NULL, *subksp = NULL;
466: KSPGetPC(ksp, &pc);
468: /* compute inner_rtol */
469: if (scale->bnrm < 0.0) {
470: Vec b;
471: KSPGetRhs(ksp, &b);
472: VecNorm(b, NORM_2, &(scale->bnrm));
473: }
474: KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
475: inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
476: /*PetscPrintf(PETSC_COMM_WORLD, " Inner rtol = %g\n", (double)inner_rtol);*/
478: /* if pc is ksp */
479: PCKSPGetKSP(pc, &kspinner);
480: if (kspinner) {
481: KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
482: return(0);
483: }
485: /* if pc is bjacobi */
486: PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
487: if (subksp) {
488: for (i=0; i<nksp; i++) {
489: KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
490: }
491: return(0);
492: }
494: /* todo: dynamic tolerance may apply to other types of pc too */
495: return(0);
496: }
500: /*
501: Destroy the dummy context used in KSPMonitorDynamicTolerance()
502: */
503: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
504: {
508: PetscFree(*dummy);
509: return(0);
510: }
514: /*
515: Default (short) KSP Monitor, same as KSPMonitorDefault() except
516: it prints fewer digits of the residual as the residual gets smaller.
517: This is because the later digits are meaningless and are often
518: different on different machines; by using this routine different
519: machines will usually generate the same output.
520: */
521: PetscErrorCode KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
522: {
524: PetscViewer viewer = (PetscViewer)dummy;
527: if (!viewer) {
528: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
529: }
530: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
531: if (its == 0 && ((PetscObject)ksp)->prefix) {
532: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
533: }
535: if (fnorm > 1.e-9) {
536: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %g \n",its,(double)fnorm);
537: } else if (fnorm > 1.e-11) {
538: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,(double)fnorm);
539: } else {
540: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);
541: }
542: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
543: return(0);
544: }
548: /*@C
549: KSPConvergedSkip - Convergence test that do not return as converged
550: until the maximum number of iterations is reached.
552: Collective on KSP
554: Input Parameters:
555: + ksp - iterative context
556: . n - iteration number
557: . rnorm - 2-norm residual value (may be estimated)
558: - dummy - unused convergence context
560: Returns:
561: . reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS
563: Notes:
564: This should be used as the convergence test with the option
565: KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
566: not computed. Convergence is then declared after the maximum number
567: of iterations have been reached. Useful when one is using CG or
568: BiCGStab as a smoother.
570: Level: advanced
572: .keywords: KSP, default, convergence, residual
574: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
575: @*/
576: PetscErrorCode KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
577: {
581: *reason = KSP_CONVERGED_ITERATING;
582: if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
583: return(0);
584: }
589: /*@C
590: KSPConvergedDefaultCreate - Creates and initializes the space used by the KSPConvergedDefault() function context
592: Collective on KSP
594: Output Parameter:
595: . ctx - convergence context
597: Level: intermediate
599: .keywords: KSP, default, convergence, residual
601: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
602: KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
603: @*/
604: PetscErrorCode KSPConvergedDefaultCreate(void **ctx)
605: {
606: PetscErrorCode ierr;
607: KSPConvergedDefaultCtx *cctx;
610: PetscNew(&cctx);
611: *ctx = cctx;
612: return(0);
613: }
617: /*@
618: KSPConvergedDefaultSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
619: instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
620: is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.
622: Collective on KSP
624: Input Parameters:
625: . ksp - iterative context
627: Options Database:
628: . -ksp_converged_use_initial_residual_norm
630: Notes:
631: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
633: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
634: are defined in petscksp.h.
636: If the convergence test is not KSPConvergedDefault() then this is ignored.
638: If right preconditioning is being used then B does not appear in the above formula.
641: Level: intermediate
643: .keywords: KSP, default, convergence, residual
645: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm()
646: @*/
647: PetscErrorCode KSPConvergedDefaultSetUIRNorm(KSP ksp)
648: {
649: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
653: if (ksp->converged != KSPConvergedDefault) return(0);
654: if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
655: ctx->initialrtol = PETSC_TRUE;
656: return(0);
657: }
661: /*@
662: KSPConvergedDefaultSetUMIRNorm - makes the default convergence test use min(|| B*(b - A*(initial guess))||,|| B*b ||)
663: In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
664: is used there is no B in the above formula. UMIRNorm is short for Use Minimum Initial Residual Norm.
666: Collective on KSP
668: Input Parameters:
669: . ksp - iterative context
671: Options Database:
672: . -ksp_converged_use_min_initial_residual_norm
674: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
676: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
677: are defined in petscksp.h.
679: Level: intermediate
681: .keywords: KSP, default, convergence, residual
683: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm()
684: @*/
685: PetscErrorCode KSPConvergedDefaultSetUMIRNorm(KSP ksp)
686: {
687: KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;
691: if (ksp->converged != KSPConvergedDefault) return(0);
692: if (ctx->initialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
693: ctx->mininitialrtol = PETSC_TRUE;
694: return(0);
695: }
699: /*@C
700: KSPConvergedDefault - Determines convergence of the linear iterative solvers by default
702: Collective on KSP
704: Input Parameters:
705: + ksp - iterative context
706: . n - iteration number
707: . rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
708: - ctx - convergence context which must be created by KSPConvergedDefaultCreate()
710: Output Parameter:
711: + positive - if the iteration has converged;
712: . negative - if residual norm exceeds divergence threshold;
713: - 0 - otherwise.
715: Notes:
716: KSPConvergedDefault() reaches convergence when rnorm < MAX (rtol * rnorm_0, abstol);
717: Divergence is detected if rnorm > dtol * rnorm_0,
719: where:
720: + rtol = relative tolerance,
721: . abstol = absolute tolerance.
722: . dtol = divergence tolerance,
723: - rnorm_0 is the two norm of the right hand side. When initial guess is non-zero you
724: can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
725: as the starting point for relative norm convergence testing, that is as rnorm_0
727: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
729: Use KSPSetNormType() (or -ksp_norm_type <none,preconditioned,unpreconditioned,natural>) to change the norm used for computing rnorm
731: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.
733: This routine is used by KSP by default so the user generally never needs call it directly.
735: Use KSPSetConvergenceTest() to provide your own test instead of using this one.
737: Level: intermediate
739: .keywords: KSP, default, convergence, residual
741: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
742: KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
743: @*/
744: PetscErrorCode KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
745: {
746: PetscErrorCode ierr;
747: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
748: KSPNormType normtype;
753: *reason = KSP_CONVERGED_ITERATING;
755: KSPGetNormType(ksp,&normtype);
756: if (normtype == KSP_NORM_NONE) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Use KSPConvergedSkip() with KSPNormType of KSP_NORM_NONE");
758: if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
759: if (!n) {
760: /* if user gives initial guess need to compute norm of b */
761: if (!ksp->guess_zero && !cctx->initialrtol) {
762: PetscReal snorm;
763: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
764: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
765: VecNorm(ksp->vec_rhs,NORM_2,&snorm); /* <- b'*b */
766: } else {
767: Vec z;
768: /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
769: VecDuplicate(ksp->vec_rhs,&z);
770: KSP_PCApply(ksp,ksp->vec_rhs,z);
771: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
772: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
773: VecNorm(z,NORM_2,&snorm); /* dp <- b'*B'*B*b */
774: } else if (ksp->normtype == KSP_NORM_NATURAL) {
775: PetscScalar norm;
776: PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
777: VecDot(ksp->vec_rhs,z,&norm);
778: snorm = PetscSqrtReal(PetscAbsScalar(norm)); /* dp <- b'*B*b */
779: }
780: VecDestroy(&z);
781: }
782: /* handle special case of zero RHS and nonzero guess */
783: if (!snorm) {
784: PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
785: snorm = rnorm;
786: }
787: if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
788: else ksp->rnorm0 = snorm;
789: } else {
790: ksp->rnorm0 = rnorm;
791: }
792: ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
793: }
795: if (n <= ksp->chknorm) return(0);
797: if (PetscIsInfOrNanReal(rnorm)) {
798: PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
799: *reason = KSP_DIVERGED_NANORINF;
800: } else if (rnorm <= ksp->ttol) {
801: if (rnorm < ksp->abstol) {
802: 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);
803: *reason = KSP_CONVERGED_ATOL;
804: } else {
805: if (cctx->initialrtol) {
806: 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);
807: } else {
808: 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);
809: }
810: *reason = KSP_CONVERGED_RTOL;
811: }
812: } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
813: 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);
814: *reason = KSP_DIVERGED_DTOL;
815: }
816: return(0);
817: }
821: /*@C
822: KSPConvergedDefaultDestroy - Frees the space used by the KSPConvergedDefault() function context
824: Collective on KSP
826: Input Parameters:
827: . ctx - convergence context
829: Level: intermediate
831: .keywords: KSP, default, convergence, residual
833: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
834: KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
835: @*/
836: PetscErrorCode KSPConvergedDefaultDestroy(void *ctx)
837: {
838: PetscErrorCode ierr;
839: KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
842: VecDestroy(&cctx->work);
843: PetscFree(ctx);
844: return(0);
845: }
849: /*
850: KSPBuildSolutionDefault - Default code to create/move the solution.
852: Input Parameters:
853: + ksp - iterative context
854: - v - pointer to the user's vector
856: Output Parameter:
857: . V - pointer to a vector containing the solution
859: Level: advanced
861: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
863: .keywords: KSP, build, solution, default
865: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
866: */
867: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
868: {
872: if (ksp->pc_side == PC_RIGHT) {
873: if (ksp->pc) {
874: if (v) {
875: KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
876: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
877: } else {
878: if (v) {
879: VecCopy(ksp->vec_sol,v); *V = v;
880: } else *V = ksp->vec_sol;
881: }
882: } else if (ksp->pc_side == PC_SYMMETRIC) {
883: if (ksp->pc) {
884: if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
885: if (v) {
886: PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
887: *V = v;
888: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
889: } else {
890: if (v) {
891: VecCopy(ksp->vec_sol,v); *V = v;
892: } else *V = ksp->vec_sol;
893: }
894: } else {
895: if (v) {
896: VecCopy(ksp->vec_sol,v); *V = v;
897: } else *V = ksp->vec_sol;
898: }
899: return(0);
900: }
904: /*
905: KSPBuildResidualDefault - Default code to compute the residual.
907: Input Parameters:
908: . ksp - iterative context
909: . t - pointer to temporary vector
910: . v - pointer to user vector
912: Output Parameter:
913: . V - pointer to a vector containing the residual
915: Level: advanced
917: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
919: .keywords: KSP, build, residual, default
921: .seealso: KSPBuildSolutionDefault()
922: */
923: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
924: {
926: Mat Amat,Pmat;
929: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
930: PCGetOperators(ksp->pc,&Amat,&Pmat);
931: KSPBuildSolution(ksp,t,NULL);
932: KSP_MatMult(ksp,Amat,t,v);
933: VecAYPX(v,-1.0,ksp->vec_rhs);
934: *V = v;
935: return(0);
936: }
940: /*@C
941: KSPCreateVecs - Gets a number of work vectors.
943: Input Parameters:
944: + ksp - iterative context
945: . rightn - number of right work vectors
946: - leftn - number of left work vectors to allocate
948: Output Parameter:
949: + right - the array of vectors created
950: - left - the array of left vectors
952: Note: The right vector has as many elements as the matrix has columns. The left
953: vector has as many elements as the matrix has rows.
955: The vectors are new vectors that are not owned by the KSP, they should be destroyed with calls to VecDestroyVecs() when no longer needed.
957: 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
958: that does not exist tries to get them from the DM (if it is provided).
960: Level: advanced
962: .seealso: MatCreateVecs(), VecDestroyVecs()
964: @*/
965: PetscErrorCode KSPCreateVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
966: {
968: Vec vecr = NULL,vecl = NULL;
969: PetscBool matset,pmatset;
970: Mat mat = NULL;
973: if (rightn) {
974: if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
975: if (ksp->vec_sol) vecr = ksp->vec_sol;
976: else {
977: if (ksp->pc) {
978: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
979: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
980: if (matset) {
981: PCGetOperators(ksp->pc,&mat,NULL);
982: MatCreateVecs(mat,&vecr,NULL);
983: } else if (pmatset) {
984: PCGetOperators(ksp->pc,NULL,&mat);
985: MatCreateVecs(mat,&vecr,NULL);
986: }
987: }
988: if (!vecr) {
989: if (ksp->dm) {
990: DMGetGlobalVector(ksp->dm,&vecr);
991: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
992: }
993: }
994: VecDuplicateVecs(vecr,rightn,right);
995: if (!ksp->vec_sol) {
996: if (mat) {
997: VecDestroy(&vecr);
998: } else if (ksp->dm) {
999: DMRestoreGlobalVector(ksp->dm,&vecr);
1000: }
1001: }
1002: }
1003: if (leftn) {
1004: if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
1005: if (ksp->vec_rhs) vecl = ksp->vec_rhs;
1006: else {
1007: if (ksp->pc) {
1008: PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
1009: /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
1010: if (matset) {
1011: PCGetOperators(ksp->pc,&mat,NULL);
1012: MatCreateVecs(mat,NULL,&vecl);
1013: } else if (pmatset) {
1014: PCGetOperators(ksp->pc,NULL,&mat);
1015: MatCreateVecs(mat,NULL,&vecl);
1016: }
1017: }
1018: if (!vecl) {
1019: if (ksp->dm) {
1020: DMGetGlobalVector(ksp->dm,&vecl);
1021: } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
1022: }
1023: }
1024: VecDuplicateVecs(vecl,leftn,left);
1025: if (!ksp->vec_rhs) {
1026: if (mat) {
1027: VecDestroy(&vecl);
1028: } else if (ksp->dm) {
1029: DMRestoreGlobalVector(ksp->dm,&vecl);
1030: }
1031: }
1032: }
1033: return(0);
1034: }
1038: /*
1039: KSPSetWorkVecs - Sets a number of work vectors into a KSP object
1041: Input Parameters:
1042: . ksp - iterative context
1043: . nw - number of work vectors to allocate
1045: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1047: */
1048: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1049: {
1053: VecDestroyVecs(ksp->nwork,&ksp->work);
1054: ksp->nwork = nw;
1055: KSPCreateVecs(ksp,nw,&ksp->work,0,NULL);
1056: PetscLogObjectParents(ksp,nw,ksp->work);
1057: return(0);
1058: }
1062: /*
1063: KSPDestroyDefault - Destroys a iterative context variable for methods with
1064: no separate context. Preferred calling sequence KSPDestroy().
1066: Input Parameter:
1067: . ksp - the iterative context
1069: Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations
1071: */
1072: PetscErrorCode KSPDestroyDefault(KSP ksp)
1073: {
1078: PetscFree(ksp->data);
1079: return(0);
1080: }
1084: /*@
1085: KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.
1087: Not Collective
1089: Input Parameter:
1090: . ksp - the KSP context
1092: Output Parameter:
1093: . reason - negative value indicates diverged, positive value converged, see KSPConvergedReason
1095: Possible values for reason:
1096: + KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1097: . KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1098: . KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence
1099: test routine is set.
1100: . KSP_CONVERGED_CG_NEG_CURVE
1101: . KSP_CONVERGED_CG_CONSTRAINED
1102: . KSP_CONVERGED_STEP_LENGTH
1103: . KSP_DIVERGED_ITS (required more than its to reach convergence)
1104: . KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1105: . KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1106: . KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1107: - KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial
1108: residual. Try a different preconditioner, or a different initial Level.)
1110: See also manual page for each reason.
1112: guess: beginner
1114: Notes: Can only be called after the call the KSPSolve() is complete.
1116: Level: intermediate
1118: .keywords: KSP, nonlinear, set, convergence, test
1120: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason
1121: @*/
1122: PetscErrorCode KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1123: {
1127: *reason = ksp->reason;
1128: return(0);
1129: }
1131: #include <petsc/private/dmimpl.h>
1134: /*@
1135: KSPSetDM - Sets the DM that may be used by some preconditioners
1137: Logically Collective on KSP
1139: Input Parameters:
1140: + ksp - the preconditioner context
1141: - dm - the dm
1143: Notes: If this is used then the KSP will attempt to use the DM to create the matrix and use the routine
1144: set with DMKSPSetComputeOperators(). Use KSPSetDMActive(ksp,PETSC_FALSE) to instead use the matrix
1145: you've provided with KSPSetOperators().
1147: Level: intermediate
1149: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1150: @*/
1151: PetscErrorCode KSPSetDM(KSP ksp,DM dm)
1152: {
1154: PC pc;
1158: if (dm) {PetscObjectReference((PetscObject)dm);}
1159: if (ksp->dm) { /* Move the DMSNES context over to the new DM unless the new DM already has one */
1160: if (ksp->dm->dmksp && ksp->dmAuto && !dm->dmksp) {
1161: DMKSP kdm;
1162: DMCopyDMKSP(ksp->dm,dm);
1163: DMGetDMKSP(ksp->dm,&kdm);
1164: if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1165: }
1166: DMDestroy(&ksp->dm);
1167: }
1168: ksp->dm = dm;
1169: ksp->dmAuto = PETSC_FALSE;
1170: KSPGetPC(ksp,&pc);
1171: PCSetDM(pc,dm);
1172: ksp->dmActive = PETSC_TRUE;
1173: return(0);
1174: }
1178: /*@
1179: KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side
1181: Logically Collective on KSP
1183: Input Parameters:
1184: + ksp - the preconditioner context
1185: - flg - use the DM
1187: Notes:
1188: 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.
1190: Level: intermediate
1192: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1193: @*/
1194: PetscErrorCode KSPSetDMActive(KSP ksp,PetscBool flg)
1195: {
1199: ksp->dmActive = flg;
1200: return(0);
1201: }
1205: /*@
1206: KSPGetDM - Gets the DM that may be used by some preconditioners
1208: Not Collective
1210: Input Parameter:
1211: . ksp - the preconditioner context
1213: Output Parameter:
1214: . dm - the dm
1216: Level: intermediate
1219: .seealso: KSPSetDM(), KSPSetDMActive()
1220: @*/
1221: PetscErrorCode KSPGetDM(KSP ksp,DM *dm)
1222: {
1227: if (!ksp->dm) {
1228: DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1229: ksp->dmAuto = PETSC_TRUE;
1230: }
1231: *dm = ksp->dm;
1232: return(0);
1233: }
1237: /*@
1238: KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.
1240: Logically Collective on KSP
1242: Input Parameters:
1243: + ksp - the KSP context
1244: - usrP - optional user context
1246: Level: intermediate
1248: .keywords: KSP, set, application, context
1250: .seealso: KSPGetApplicationContext()
1251: @*/
1252: PetscErrorCode KSPSetApplicationContext(KSP ksp,void *usrP)
1253: {
1255: PC pc;
1259: ksp->user = usrP;
1260: KSPGetPC(ksp,&pc);
1261: PCSetApplicationContext(pc,usrP);
1262: return(0);
1263: }
1267: /*@
1268: KSPGetApplicationContext - Gets the user-defined context for the linear solver.
1270: Not Collective
1272: Input Parameter:
1273: . ksp - KSP context
1275: Output Parameter:
1276: . usrP - user context
1278: Level: intermediate
1280: .keywords: KSP, get, application, context
1282: .seealso: KSPSetApplicationContext()
1283: @*/
1284: PetscErrorCode KSPGetApplicationContext(KSP ksp,void *usrP)
1285: {
1288: *(void**)usrP = ksp->user;
1289: return(0);
1290: }