Actual source code: iterativ.c
petsc-3.3-p7 2013-05-11
1: /*
2: This file contains some simple default routines.
3: These routines should be SHORT, since they will be included in every
4: executable image that uses the iterative routines (note that, through
5: the registry system, we provide a way to load only the truely necessary
6: files)
7: */
8: #include <petsc-private/kspimpl.h> /*I "petscksp.h" I*/
9: #include <petscdmshell.h>
13: /*@
14: KSPGetResidualNorm - Gets the last (approximate preconditioned)
15: residual norm that has been computed.
16:
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.
46:
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.
78:
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;
103: PetscErrorCode ierr;
104: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm);
108: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
109: if (!ksp->calc_sings) {
110: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
111: } else {
112: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
113: c = emax/emin;
114: 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);
115: }
116: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
117: return(0);
118: }
122: /*@C
123: KSPMonitorSolution - Monitors progress of the KSP solvers by calling
124: VecView() for the approximate solution at each iteration.
126: Collective on KSP
128: Input Parameters:
129: + ksp - the KSP context
130: . its - iteration number
131: . fgnorm - 2-norm of residual (or gradient)
132: - dummy - either a viewer or PETSC_NULL
134: Level: intermediate
136: Notes:
137: For some Krylov methods such as GMRES constructing the solution at
138: each iteration is expensive, hence using this will slow the code.
140: .keywords: KSP, nonlinear, vector, monitor, view
142: .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView()
143: @*/
144: PetscErrorCode KSPMonitorSolution(KSP ksp,PetscInt its,PetscReal fgnorm,void *dummy)
145: {
147: Vec x;
148: PetscViewer viewer = (PetscViewer) dummy;
151: KSPBuildSolution(ksp,PETSC_NULL,&x);
152: if (!viewer) {
153: MPI_Comm comm;
154: PetscObjectGetComm((PetscObject)ksp,&comm);
155: viewer = PETSC_VIEWER_DRAW_(comm);
156: }
157: VecView(x,viewer);
159: return(0);
160: }
164: /*@C
165: KSPMonitorDefault - Print the residual norm at each iteration of an
166: iterative solver.
168: Collective on KSP
170: Input Parameters:
171: + ksp - iterative context
172: . n - iteration number
173: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
174: - dummy - unused monitor context
176: Level: intermediate
178: .keywords: KSP, default, monitor, residual
180: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGCreate()
181: @*/
182: PetscErrorCode KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
183: {
185: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm);
188: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
189: if (n == 0 && ((PetscObject)ksp)->prefix) {
190: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
191: }
192: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
193: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
194: return(0);
195: }
199: /*@C
200: KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
201: residual norm at each iteration of an iterative solver.
203: Collective on KSP
205: Input Parameters:
206: + ksp - iterative context
207: . n - iteration number
208: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
209: - dummy - unused monitor context
211: Options Database Key:
212: . -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()
214: Notes:
215: When using right preconditioning, these values are equivalent.
217: Level: intermediate
219: .keywords: KSP, default, monitor, residual
221: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGCreate()
222: @*/
223: PetscErrorCode KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
224: {
225: PetscErrorCode ierr;
226: Vec resid,work;
227: PetscReal scnorm,bnorm;
228: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm);
229: char normtype[256];
232: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
233: if (n == 0 && ((PetscObject)ksp)->prefix) {
234: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
235: }
236: VecDuplicate(ksp->vec_rhs,&work);
237: KSPBuildResidual(ksp,0,work,&resid);
239: /*
240: Unscale the residual but only if both matrices are the same matrix, since only then would
241: they be scaled.
242: */
243: VecCopy(resid,work);
244: VecNorm(work,NORM_2,&scnorm);
245: VecDestroy(&work);
246: VecNorm(ksp->vec_rhs,NORM_2,&bnorm);
247: PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof normtype);
248: PetscStrtolower(normtype);
249: 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)scnorm,(double)(scnorm/bnorm));
250: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
251: return(0);
252: }
256: PetscErrorCode KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
257: {
258: PetscErrorCode ierr;
259: Vec resid,work;
260: PetscReal rmax,pwork;
261: PetscInt i,n,N;
262: PetscScalar *r;
265: VecDuplicate(ksp->vec_rhs,&work);
266: KSPBuildResidual(ksp,0,work,&resid);
268: /*
269: Unscale the residual if the matrix is, but only if both matrices are the same matrix, since only then would
270: they be scaled.
271: */
272: VecCopy(resid,work);
273: VecNorm(work,NORM_INFINITY,&rmax);
274: VecGetLocalSize(work,&n);
275: VecGetSize(work,&N);
276: VecGetArray(work,&r);
277: pwork = 0.0;
278: for (i=0; i<n; i++) {
279: pwork += (PetscAbsScalar(r[i]) > .20*rmax);
280: }
281: MPI_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,((PetscObject)ksp)->comm);
282: VecRestoreArray(work,&r);
283: VecDestroy(&work);
284: *per = *per/N;
285: return(0);
286: }
290: /*@C
291: KSPMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.
293: Collective on KSP
295: Input Parameters:
296: + ksp - iterative context
297: . it - iteration number
298: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
299: - dummy - unused monitor context
301: Options Database Key:
302: . -ksp_monitor_range - Activates KSPMonitorRange()
304: Level: intermediate
306: .keywords: KSP, default, monitor, residual
308: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGCreate()
309: @*/
310: PetscErrorCode KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,void *dummy)
311: {
312: PetscErrorCode ierr;
313: PetscReal perc,rel;
314: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm);
315: /* should be in a MonitorRangeContext */
316: static PetscReal prev;
319: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
320: if (!it) prev = rnorm;
321: if (it == 0 && ((PetscObject)ksp)->prefix) {
322: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
323: }
324: KSPMonitorRange_Private(ksp,it,&perc);
326: rel = (prev - rnorm)/prev;
327: prev = rnorm;
328: 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);
329: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
330: return(0);
331: }
335: /*
336: Default (short) KSP Monitor, same as KSPMonitorDefault() except
337: it prints fewer digits of the residual as the residual gets smaller.
338: This is because the later digits are meaningless and are often
339: different on different machines; by using this routine different
340: machines will usually generate the same output.
341: */
342: PetscErrorCode KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
343: {
345: PetscViewer viewer = dummy ? (PetscViewer) dummy : PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm);
348: PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
349: if (its == 0 && ((PetscObject)ksp)->prefix) {
350: PetscViewerASCIIPrintf(viewer," Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
351: }
353: if (fnorm > 1.e-9) {
354: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %G \n",its,fnorm);
355: } else if (fnorm > 1.e-11){
356: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,(double)fnorm);
357: } else {
358: PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);
359: }
360: PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
361: return(0);
362: }
366: /*@C
367: KSPSkipConverged - Convergence test that do not return as converged
368: until the maximum number of iterations is reached.
370: Collective on KSP
372: Input Parameters:
373: + ksp - iterative context
374: . n - iteration number
375: . rnorm - 2-norm residual value (may be estimated)
376: - dummy - unused convergence context
378: Returns:
379: . reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS
381: Notes:
382: This should be used as the convergence test with the option
383: KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
384: not computed. Convergence is then declared after the maximum number
385: of iterations have been reached. Useful when one is using CG or
386: BiCGStab as a smoother.
387:
388: Level: advanced
390: .keywords: KSP, default, convergence, residual
392: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
393: @*/
394: PetscErrorCode KSPSkipConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
395: {
399: *reason = KSP_CONVERGED_ITERATING;
400: if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
401: return(0);
402: }
407: /*@C
408: KSPDefaultConvergedCreate - Creates and initializes the space used by the KSPDefaultConverged() function context
410: Collective on KSP
412: Output Parameter:
413: . ctx - convergence context
415: Level: intermediate
417: .keywords: KSP, default, convergence, residual
419: .seealso: KSPDefaultConverged(), KSPDefaultConvergedDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
420: KSPSkipConverged(), KSPConvergedReason, KSPGetConvergedReason(), KSPDefaultConvergedSetUIRNorm(), KSPDefaultConvergedSetUMIRNorm()
421: @*/
422: PetscErrorCode KSPDefaultConvergedCreate(void **ctx)
423: {
424: PetscErrorCode ierr;
425: KSPDefaultConvergedCtx *cctx;
428: PetscNew(KSPDefaultConvergedCtx,&cctx);
429: *ctx = cctx;
430: return(0);
431: }
435: /*@
436: KSPDefaultConvergedSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
437: instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
438: is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.
440: Collective on KSP
442: Input Parameters:
443: . ksp - iterative context
445: Options Database:
446: . -ksp_converged_use_initial_residual_norm
448: Notes:
449: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
451: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
452: are defined in petscksp.h.
454: If the convergence test is not KSPDefaultConverged() then this is ignored.
456: If right preconditioning is being used then B does not appear in the above formula.
458:
459: Level: intermediate
461: .keywords: KSP, default, convergence, residual
463: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSkipConverged(), KSPConvergedReason, KSPGetConvergedReason(), KSPDefaultConvergedSetUMIRNorm()
464: @*/
465: PetscErrorCode KSPDefaultConvergedSetUIRNorm(KSP ksp)
466: {
467: KSPDefaultConvergedCtx *ctx = (KSPDefaultConvergedCtx*) ksp->cnvP;
471: if (ksp->converged != KSPDefaultConverged) return(0);
472: if (ctx->mininitialrtol) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPDefaultConvergedSetUIRNorm() and KSPDefaultConvergedSetUMIRNorm() together");
473: ctx->initialrtol = PETSC_TRUE;
474: return(0);
475: }
479: /*@
480: KSPDefaultConvergedSetUMIRNorm - makes the default convergence test use min(|| B*(b - A*(initial guess))||,|| B*b ||)
481: In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
482: is used there is no B in the above formula. UMIRNorm is short for Use Minimum Initial Residual Norm.
484: Collective on KSP
486: Input Parameters:
487: . ksp - iterative context
489: Options Database:
490: . -ksp_converged_use_min_initial_residual_norm
492: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
494: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
495: are defined in petscksp.h.
497: Level: intermediate
499: .keywords: KSP, default, convergence, residual
501: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSkipConverged(), KSPConvergedReason, KSPGetConvergedReason(), KSPDefaultConvergedSetUIRNorm()
502: @*/
503: PetscErrorCode KSPDefaultConvergedSetUMIRNorm(KSP ksp)
504: {
505: KSPDefaultConvergedCtx *ctx = (KSPDefaultConvergedCtx*) ksp->cnvP;
509: if (ksp->converged != KSPDefaultConverged) return(0);
510: if (ctx->initialrtol) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPDefaultConvergedSetUIRNorm() and KSPDefaultConvergedSetUMIRNorm() together");
511: ctx->mininitialrtol = PETSC_TRUE;
512: return(0);
513: }
517: /*@C
518: KSPDefaultConverged - Determines convergence of
519: the iterative solvers (default code).
521: Collective on KSP
523: Input Parameters:
524: + ksp - iterative context
525: . n - iteration number
526: . rnorm - 2-norm residual value (may be estimated)
527: - ctx - convergence context which must be created by KSPDefaultConvergedCreate()
529: reason is set to:
530: + positive - if the iteration has converged;
531: . negative - if residual norm exceeds divergence threshold;
532: - 0 - otherwise.
534: Notes:
535: KSPDefaultConverged() reaches convergence when
536: $ rnorm < MAX (rtol * rnorm_0, abstol);
537: Divergence is detected if
538: $ rnorm > dtol * rnorm_0,
540: where
541: + rtol = relative tolerance,
542: . abstol = absolute tolerance.
543: . dtol = divergence tolerance,
544: - rnorm_0 is the two norm of the right hand side. When initial guess is non-zero you
545: can call KSPDefaultConvergedSetUIRNorm() to use the norm of (b - A*(initial guess))
546: as the starting point for relative norm convergence testing.
548: Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.
550: The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
551: are defined in petscksp.h.
553: Level: intermediate
555: .keywords: KSP, default, convergence, residual
557: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSkipConverged(), KSPConvergedReason, KSPGetConvergedReason(),
558: KSPDefaultConvergedSetUIRNorm(), KSPDefaultConvergedSetUMIRNorm(), KSPDefaultConvergedCreate(), KSPDefaultConvergedDestroy()
559: @*/
560: PetscErrorCode KSPDefaultConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
561: {
562: PetscErrorCode ierr;
563: KSPDefaultConvergedCtx *cctx = (KSPDefaultConvergedCtx*) ctx;
564: KSPNormType normtype;
569: *reason = KSP_CONVERGED_ITERATING;
570:
571: KSPGetNormType(ksp,&normtype);
572: if (normtype == KSP_NORM_NONE) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_WRONGSTATE,"Use KSPSkipConverged() with KSPNormType of KSP_NORM_NONE");
574: if (!cctx) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPDefaultConvergedCreate()");
575: if (!n) {
576: /* if user gives initial guess need to compute norm of b */
577: if (!ksp->guess_zero && !cctx->initialrtol) {
578: PetscReal snorm;
579: if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
580: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
581: VecNorm(ksp->vec_rhs,NORM_2,&snorm); /* <- b'*b */
582: } else {
583: Vec z;
584: /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
585: VecDuplicate(ksp->vec_rhs,&z);
586: KSP_PCApply(ksp,ksp->vec_rhs,z);
587: if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
588: PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
589: VecNorm(z,NORM_2,&snorm); /* dp <- b'*B'*B*b */
590: } else if (ksp->normtype == KSP_NORM_NATURAL) {
591: PetscScalar norm;
592: PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
593: VecDot(ksp->vec_rhs,z,&norm);
594: snorm = PetscSqrtReal(PetscAbsScalar(norm)); /* dp <- b'*B*b */
595: }
596: VecDestroy(&z);
597: }
598: /* handle special case of zero RHS and nonzero guess */
599: if (!snorm) {
600: PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
601: snorm = rnorm;
602: }
603: if (cctx->mininitialrtol) {
604: ksp->rnorm0 = PetscMin(snorm,rnorm);
605: } else {
606: ksp->rnorm0 = snorm;
607: }
608: } else {
609: ksp->rnorm0 = rnorm;
610: }
611: ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
612: }
614: if (n <= ksp->chknorm) return(0);
616: if (PetscIsInfOrNanScalar(rnorm)) {
617: PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
618: *reason = KSP_DIVERGED_NAN;
619: } else if (rnorm <= ksp->ttol) {
620: if (rnorm < ksp->abstol) {
621: 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);
622: *reason = KSP_CONVERGED_ATOL;
623: } else {
624: if (cctx->initialrtol) {
625: 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);
626: } else {
627: 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);
628: }
629: *reason = KSP_CONVERGED_RTOL;
630: }
631: } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
632: 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);
633: *reason = KSP_DIVERGED_DTOL;
634: }
635: return(0);
636: }
640: /*@C
641: KSPDefaultConvergedDestroy - Frees the space used by the KSPDefaultConverged() function context
643: Collective on KSP
645: Input Parameters:
646: . ctx - convergence context
648: Level: intermediate
650: .keywords: KSP, default, convergence, residual
652: .seealso: KSPDefaultConverged(), KSPDefaultConvergedCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPSkipConverged(),
653: KSPConvergedReason, KSPGetConvergedReason(), KSPDefaultConvergedSetUIRNorm(), KSPDefaultConvergedSetUMIRNorm()
654: @*/
655: PetscErrorCode KSPDefaultConvergedDestroy(void *ctx)
656: {
657: PetscErrorCode ierr;
658: KSPDefaultConvergedCtx *cctx = (KSPDefaultConvergedCtx*) ctx;
661: VecDestroy(&cctx->work);
662: PetscFree(ctx);
663: return(0);
664: }
668: /*
669: KSPDefaultBuildSolution - Default code to create/move the solution.
671: Input Parameters:
672: + ksp - iterative context
673: - v - pointer to the user's vector
675: Output Parameter:
676: . V - pointer to a vector containing the solution
678: Level: advanced
680: .keywords: KSP, build, solution, default
682: .seealso: KSPGetSolution(), KSPDefaultBuildResidual()
683: */
684: PetscErrorCode KSPDefaultBuildSolution(KSP ksp,Vec v,Vec *V)
685: {
688: if (ksp->pc_side == PC_RIGHT) {
689: if (ksp->pc) {
690: if (v) {KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;}
691: else SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Not working with right preconditioner");
692: } else {
693: if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
694: else { *V = ksp->vec_sol;}
695: }
696: } else if (ksp->pc_side == PC_SYMMETRIC) {
697: if (ksp->pc) {
698: if (ksp->transpose_solve) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
699: if (v) {PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v); *V = v;}
700: else SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"Not working with symmetric preconditioner");
701: } else {
702: if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
703: else { *V = ksp->vec_sol;}
704: }
705: } else {
706: if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
707: else { *V = ksp->vec_sol; }
708: }
709: return(0);
710: }
714: /*
715: KSPDefaultBuildResidual - Default code to compute the residual.
717: Input Parameters:
718: . ksp - iterative context
719: . t - pointer to temporary vector
720: . v - pointer to user vector
722: Output Parameter:
723: . V - pointer to a vector containing the residual
725: Level: advanced
727: .keywords: KSP, build, residual, default
729: .seealso: KSPDefaultBuildSolution()
730: */
731: PetscErrorCode KSPDefaultBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
732: {
734: MatStructure pflag;
735: Mat Amat,Pmat;
738: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
739: PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
740: KSPBuildSolution(ksp,t,PETSC_NULL);
741: KSP_MatMult(ksp,Amat,t,v);
742: VecAYPX(v,-1.0,ksp->vec_rhs);
743: *V = v;
744: return(0);
745: }
749: /*@C
750: KSPGetVecs - Gets a number of work vectors.
752: Input Parameters:
753: + ksp - iterative context
754: . rightn - number of right work vectors
755: - leftn - number of left work vectors to allocate
757: Output Parameter:
758: + right - the array of vectors created
759: - left - the array of left vectors
761: Note: The right vector has as many elements as the matrix has columns. The left
762: vector has as many elements as the matrix has rows.
764: Level: advanced
766: .seealso: MatGetVecs()
768: @*/
769: PetscErrorCode KSPGetVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
770: {
772: Vec vecr,vecl;
775: if (rightn) {
776: if (!right) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
777: if (ksp->vec_sol) vecr = ksp->vec_sol;
778: else {
779: if (ksp->dm) {
780: DMGetGlobalVector(ksp->dm,&vecr);
781: } else {
782: Mat mat;
783: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
784: PCGetOperators(ksp->pc,&mat,PETSC_NULL,PETSC_NULL);
785: MatGetVecs(mat,&vecr,PETSC_NULL);
786: }
787: }
788: VecDuplicateVecs(vecr,rightn,right);
789: if (!ksp->vec_sol) {
790: if (ksp->dm) {
791: DMRestoreGlobalVector(ksp->dm,&vecr);
792: } else {
793: VecDestroy(&vecr);
794: }
795: }
796: }
797: if (leftn) {
798: if (!left) SETERRQ(((PetscObject)ksp)->comm,PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
799: if (ksp->vec_rhs) vecl = ksp->vec_rhs;
800: else {
801: if (ksp->dm) {
802: DMGetGlobalVector(ksp->dm,&vecl);
803: } else {
804: Mat mat;
805: if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
806: PCGetOperators(ksp->pc,&mat,PETSC_NULL,PETSC_NULL);
807: MatGetVecs(mat,PETSC_NULL,&vecl);
808: }
809: }
810: VecDuplicateVecs(vecl,leftn,left);
811: if (!ksp->vec_rhs) {
812: if (ksp->dm) {
813: DMRestoreGlobalVector(ksp->dm,&vecl);
814: } else {
815: VecDestroy(&vecl);
816: }
817: }
818: }
819: return(0);
820: }
824: /*
825: KSPDefaultGetWork - Gets a number of work vectors.
827: Input Parameters:
828: . ksp - iterative context
829: . nw - number of work vectors to allocate
831: Notes:
832: Call this only if no work vectors have been allocated
833: */
834: PetscErrorCode KSPDefaultGetWork(KSP ksp,PetscInt nw)
835: {
839: VecDestroyVecs(ksp->nwork,&ksp->work);
840: ksp->nwork = nw;
841: KSPGetVecs(ksp,nw,&ksp->work,0,PETSC_NULL);
842: PetscLogObjectParents(ksp,nw,ksp->work);
843: return(0);
844: }
848: /*
849: KSPDefaultDestroy - Destroys a iterative context variable for methods with
850: no separate context. Preferred calling sequence KSPDestroy().
852: Input Parameter:
853: . ksp - the iterative context
854: */
855: PetscErrorCode KSPDefaultDestroy(KSP ksp)
856: {
861: PetscFree(ksp->data);
862: return(0);
863: }
867: /*@
868: KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.
870: Not Collective
872: Input Parameter:
873: . ksp - the KSP context
875: Output Parameter:
876: . reason - negative value indicates diverged, positive value converged, see KSPConvergedReason
878: Possible values for reason:
879: + KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
880: . KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
881: . KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPSkipConverged() convergence
882: test routine is set.
883: . KSP_CONVERGED_CG_NEG_CURVE
884: . KSP_CONVERGED_CG_CONSTRAINED
885: . KSP_CONVERGED_STEP_LENGTH
886: . KSP_DIVERGED_ITS (required more than its to reach convergence)
887: . KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
888: . KSP_DIVERGED_NAN (residual norm became Not-a-number likely due to 0/0)
889: . KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
890: - KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial
891: residual. Try a different preconditioner, or a different initial Level.)
892:
893: See also manual page for each reason.
895: guess: beginner
897: Notes: Can only be called after the call the KSPSolve() is complete.
899: Level: intermediate
900:
901: .keywords: KSP, nonlinear, set, convergence, test
903: .seealso: KSPSetConvergenceTest(), KSPDefaultConverged(), KSPSetTolerances(), KSPConvergedReason
904: @*/
905: PetscErrorCode KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
906: {
910: *reason = ksp->reason;
911: return(0);
912: }
916: /*@
917: KSPSetDM - Sets the DM that may be used by some preconditioners
919: Logically Collective on KSP
921: Input Parameters:
922: + ksp - the preconditioner context
923: - dm - the dm
925: Level: intermediate
928: .seealso: KSPGetDM(), KSPSetDM(), KSPGetDM()
929: @*/
930: PetscErrorCode KSPSetDM(KSP ksp,DM dm)
931: {
933: PC pc;
937: if (dm) {PetscObjectReference((PetscObject)dm);}
938: if (ksp->dm) { /* Move the SNESDM context over to the new DM unless the new DM already has one */
939: PetscContainer oldcontainer,container;
940: KSPDM kdm;
941: PetscObjectQuery((PetscObject)ksp->dm,"KSPDM",(PetscObject*)&oldcontainer);
942: PetscObjectQuery((PetscObject)dm,"KSPDM",(PetscObject*)&container);
943: if (oldcontainer && !container) {
944: DMKSPCopyContext(ksp->dm,dm);
945: DMKSPGetContext(ksp->dm,&kdm);
946: if (kdm->originaldm == ksp->dm) { /* Grant write privileges to the replacement DM */
947: kdm->originaldm = dm;
948: }
949: }
950: DMDestroy(&ksp->dm);
951: }
952: ksp->dm = dm;
953: KSPGetPC(ksp,&pc);
954: PCSetDM(pc,dm);
955: ksp->dmActive = PETSC_TRUE;
956: return(0);
957: }
961: /*@
962: KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side
964: Logically Collective on KSP
966: Input Parameters:
967: + ksp - the preconditioner context
968: - flg - use the DM
970: Level: intermediate
972: Notes:
973: By default KSPSetDM() sets the DM as active, call KSPSetDMActive(dm,PETSC_FALSE); after KSPSetDM(dm) to not have the KSP object use the DM to generate the matrices
975: .seealso: KSPGetDM(), KSPSetDM(), KSPGetDM()
976: @*/
977: PetscErrorCode KSPSetDMActive(KSP ksp,PetscBool flg)
978: {
982: ksp->dmActive = flg;
983: return(0);
984: }
988: /*@
989: KSPGetDM - Gets the DM that may be used by some preconditioners
991: Not Collective
993: Input Parameter:
994: . ksp - the preconditioner context
996: Output Parameter:
997: . dm - the dm
999: Level: intermediate
1002: .seealso: KSPSetDM(), KSPSetDM(), KSPGetDM()
1003: @*/
1004: PetscErrorCode KSPGetDM(KSP ksp,DM *dm)
1005: {
1010: if (!ksp->dm) {
1011: DMShellCreate(((PetscObject)ksp)->comm,&ksp->dm);
1012: }
1013: *dm = ksp->dm;
1014: return(0);
1015: }
1019: /*@
1020: KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.
1022: Logically Collective on KSP
1024: Input Parameters:
1025: + ksp - the KSP context
1026: - usrP - optional user context
1028: Level: intermediate
1030: .keywords: KSP, set, application, context
1032: .seealso: KSPGetApplicationContext()
1033: @*/
1034: PetscErrorCode KSPSetApplicationContext(KSP ksp,void *usrP)
1035: {
1037: PC pc;
1041: ksp->user = usrP;
1042: KSPGetPC(ksp,&pc);
1043: PCSetApplicationContext(pc,usrP);
1044: return(0);
1045: }
1049: /*@
1050: KSPGetApplicationContext - Gets the user-defined context for the linear solver.
1052: Not Collective
1054: Input Parameter:
1055: . ksp - KSP context
1057: Output Parameter:
1058: . usrP - user context
1060: Level: intermediate
1062: .keywords: KSP, get, application, context
1064: .seealso: KSPSetApplicationContext()
1065: @*/
1066: PetscErrorCode KSPGetApplicationContext(KSP ksp,void *usrP)
1067: {
1070: *(void**)usrP = ksp->user;
1071: return(0);
1072: }