Actual source code: iterativ.c
1: /*$Id: iterativ.c,v 1.108 2001/08/07 03:03:45 balay Exp $*/
3: /*
4: This file contains some simple default routines.
5: These routines should be SHORT, since they will be included in every
6: executable image that uses the iterative routines (note that, through
7: the registry system, we provide a way to load only the truely necessary
8: files)
9: */
10: #include src/sles/ksp/kspimpl.h
12: #undef __FUNCT__
14: /*
15: KSPDefaultFreeWork - Free work vectors
17: Input Parameters:
18: . ksp - iterative context
19: */
20: int KSPDefaultFreeWork(KSP ksp)
21: {
25: if (ksp->work) {
26: VecDestroyVecs(ksp->work,ksp->nwork);
27: }
28: return(0);
29: }
31: #undef __FUNCT__
33: /*@C
34: KSPGetResidualNorm - Gets the last (approximate preconditioned)
35: residual norm that has been computed.
36:
37: Not Collective
39: Input Parameters:
40: . ksp - the iterative context
42: Output Parameters:
43: . rnorm - residual norm
45: Level: intermediate
47: .keywords: KSP, get, residual norm
49: .seealso: KSPComputeResidual()
50: @*/
51: int KSPGetResidualNorm(KSP ksp,PetscReal *rnorm)
52: {
55: *rnorm = ksp->rnorm;
56: return(0);
57: }
59: #undef __FUNCT__
61: /*@
62: KSPGetIterationNumber - Gets the current iteration number (if the
63: KSPSolve() (SLESSolve()) is complete, returns the number of iterations
64: used.
65:
66: Not Collective
68: Input Parameters:
69: . ksp - the iterative context
71: Output Parameters:
72: . its - number of iterations
74: Level: intermediate
76: Notes:
77: During the ith iteration this returns i-1
78: .keywords: KSP, get, residual norm
80: .seealso: KSPComputeResidual(), KSPGetResidualNorm()
81: @*/
82: int KSPGetIterationNumber(KSP ksp,int *its)
83: {
86: *its = ksp->its;
87: return(0);
88: }
90: #undef __FUNCT__
92: /*@C
93: KSPSingularValueMonitor - Prints the two norm of the true residual and
94: estimation of the extreme singular values of the preconditioned problem
95: at each iteration.
96:
97: Collective on KSP
99: Input Parameters:
100: + ksp - the iterative context
101: . n - the iteration
102: - rnorm - the two norm of the residual
104: Options Database Key:
105: . -ksp_singmonitor - Activates KSPSingularValueMonitor()
107: Notes:
108: The CG solver uses the Lanczos technique for eigenvalue computation,
109: while GMRES uses the Arnoldi technique; other iterative methods do
110: not currently compute singular values.
112: Level: intermediate
114: .keywords: KSP, CG, default, monitor, extreme, singular values, Lanczos, Arnoldi
116: .seealso: KSPComputeExtremeSingularValues()
117: @*/
118: int KSPSingularValueMonitor(KSP ksp,int n,PetscReal rnorm,void *dummy)
119: {
120: PetscReal emin,emax,c;
121: int ierr;
125: if (!ksp->calc_sings) {
126: PetscPrintf(ksp->comm,"%3d KSP Residual norm %14.12e n",n,rnorm);
127: } else {
128: KSPComputeExtremeSingularValues(ksp,&emax,&emin);
129: c = emax/emin;
130: PetscPrintf(ksp->comm,"%3d KSP Residual norm %14.12e %% max %g min %g max/min %gn",n,rnorm,emax,emin,c);
131: }
132: return(0);
133: }
135: #undef __FUNCT__
137: /*@C
138: KSPVecViewMonitor - Monitors progress of the KSP solvers by calling
139: VecView() for the approximate solution at each iteration.
141: Collective on KSP
143: Input Parameters:
144: + ksp - the KSP context
145: . its - iteration number
146: . fgnorm - 2-norm of residual (or gradient)
147: - dummy - either a viewer or PETSC_NULL
149: Level: intermediate
151: Notes:
152: For some Krylov methods such as GMRES constructing the solution at
153: each iteration is expensive, hence using this will slow the code.
155: .keywords: KSP, nonlinear, vector, monitor, view
157: .seealso: KSPSetMonitor(), KSPDefaultMonitor(), VecView()
158: @*/
159: int KSPVecViewMonitor(KSP ksp,int its,PetscReal fgnorm,void *dummy)
160: {
161: int ierr;
162: Vec x;
163: PetscViewer viewer = (PetscViewer) dummy;
166: KSPBuildSolution(ksp,PETSC_NULL,&x);
167: if (!viewer) {
168: MPI_Comm comm;
169: ierr = PetscObjectGetComm((PetscObject)ksp,&comm);
170: viewer = PETSC_VIEWER_DRAW_(comm);
171: }
172: VecView(x,viewer);
174: return(0);
175: }
177: #undef __FUNCT__
179: /*@C
180: KSPDefaultMonitor - Print the residual norm at each iteration of an
181: iterative solver.
183: Collective on KSP
185: Input Parameters:
186: + ksp - iterative context
187: . n - iteration number
188: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
189: - dummy - unused monitor context
191: Level: intermediate
193: .keywords: KSP, default, monitor, residual
195: .seealso: KSPSetMonitor(), KSPTrueMonitor(), KSPLGMonitorCreate()
196: @*/
197: int KSPDefaultMonitor(KSP ksp,int n,PetscReal rnorm,void *dummy)
198: {
199: int ierr;
200: PetscViewer viewer = (PetscViewer) dummy;
203: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(ksp->comm);
204: PetscViewerASCIIPrintf(viewer,"%3d KSP Residual norm %14.12e n",n,rnorm);
205: return(0);
206: }
208: #undef __FUNCT__
210: /*@C
211: KSPTrueMonitor - Prints the true residual norm as well as the preconditioned
212: residual norm at each iteration of an iterative solver.
214: Collective on KSP
216: Input Parameters:
217: + ksp - iterative context
218: . n - iteration number
219: . rnorm - 2-norm (preconditioned) residual value (may be estimated).
220: - dummy - unused monitor context
222: Options Database Key:
223: . -ksp_truemonitor - Activates KSPTrueMonitor()
225: Notes:
226: When using right preconditioning, these values are equivalent.
228: When using either ICC or ILU preconditioners in BlockSolve95
229: (via MATMPIROWBS matrix format), then use this monitor will
230: print both the residual norm associated with the original
231: (unscaled) matrix.
233: Level: intermediate
235: .keywords: KSP, default, monitor, residual
237: .seealso: KSPSetMonitor(), KSPDefaultMonitor(), KSPLGMonitorCreate()
238: @*/
239: int KSPTrueMonitor(KSP ksp,int n,PetscReal rnorm,void *dummy)
240: {
241: int ierr;
242: Vec resid,work;
243: PetscReal scnorm;
244: PC pc;
245: Mat A,B;
246: PetscViewer viewer = (PetscViewer) dummy;
247:
249: VecDuplicate(ksp->vec_rhs,&work);
250: KSPBuildResidual(ksp,0,work,&resid);
252: /*
253: Unscale the residual if the matrix is, for example, a BlockSolve matrix
254: but only if both matrices are the same matrix, since only then would
255: they be scaled.
256: */
257: VecCopy(resid,work);
258: KSPGetPC(ksp,&pc);
259: PCGetOperators(pc,&A,&B,PETSC_NULL);
260: if (A == B) {
261: MatUnScaleSystem(A,PETSC_NULL,work);
262: }
263: VecNorm(work,NORM_2,&scnorm);
264: VecDestroy(work);
265: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(ksp->comm);
266: PetscViewerASCIIPrintf(viewer,"%3d KSP preconditioned resid norm %14.12e true resid norm %14.12en",n,rnorm,scnorm);
267: return(0);
268: }
270: #undef __FUNCT__
272: /*
273: Default (short) KSP Monitor, same as KSPDefaultMonitor() except
274: it prints fewer digits of the residual as the residual gets smaller.
275: This is because the later digits are meaningless and are often
276: different on different machines; by using this routine different
277: machines will usually generate the same output.
278: */
279: int KSPDefaultSMonitor(KSP ksp,int its,PetscReal fnorm,void *dummy)
280: {
281: int ierr;
282: PetscViewer viewer = (PetscViewer) dummy;
285: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(ksp->comm);
287: if (fnorm > 1.e-9) {
288: PetscViewerASCIIPrintf(viewer,"%3d KSP Residual norm %g n",its,fnorm);
289: } else if (fnorm > 1.e-11){
290: PetscViewerASCIIPrintf(viewer,"%3d KSP Residual norm %5.3e n",its,fnorm);
291: } else {
292: PetscViewerASCIIPrintf(viewer,"%3d KSP Residual norm < 1.e-11n",its);
293: }
294: return(0);
295: }
297: #undef __FUNCT__
299: /*@C
300: KSPSkipConverged - Convergence test that NEVER returns as converged.
302: Collective on KSP
304: Input Parameters:
305: + ksp - iterative context
306: . n - iteration number
307: . rnorm - 2-norm residual value (may be estimated)
308: - dummy - unused convergence context
310: Returns:
311: . 0 - always
313: Notes:
314: This is used as the convergence test with the option KSPSetNormType(ksp,KSP_NO_NORM),
315: since norms of the residual are not computed. Convergence is then declared
316: after a fixed number of iterations have been used. Useful when one is
317: using CG or Bi-CG-stab as a smoother.
318:
319: Level: advanced
321: .keywords: KSP, default, convergence, residual
323: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
324: @*/
325: int KSPSkipConverged(KSP ksp,int n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
326: {
329: return(0);
330: }
332: #undef __FUNCT__
334: /*@C
335: KSPDefaultConverged - Determines convergence of
336: the iterative solvers (default code).
338: Collective on KSP
340: Input Parameters:
341: + ksp - iterative context
342: . n - iteration number
343: . rnorm - 2-norm residual value (may be estimated)
344: - dummy - unused convergence context
346: Returns:
347: + 1 - if the iteration has converged;
348: . -1 - if residual norm exceeds divergence threshold;
349: - 0 - otherwise.
351: Notes:
352: KSPDefaultConverged() reaches convergence when
353: $ rnorm < MAX (rtol * rnorm_0, atol);
354: Divergence is detected if
355: $ rnorm > dtol * rnorm_0,
357: where
358: + rtol = relative tolerance,
359: . atol = absolute tolerance.
360: . dtol = divergence tolerance,
361: - rnorm_0 = initial residual norm
363: Use KSPSetTolerances() to alter the defaults for rtol, atol, dtol.
365: Level: intermediate
367: .keywords: KSP, default, convergence, residual
369: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSkipConverged()
370: @*/
371: int KSPDefaultConverged(KSP ksp,int n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
372: {
375: *reason = KSP_CONVERGED_ITERATING;
377: if (!n) {
378: ksp->ttol = PetscMax(ksp->rtol*rnorm,ksp->atol);
379: ksp->rnorm0 = rnorm;
380: }
381: if (rnorm <= ksp->ttol) {
382: if (rnorm < ksp->atol) {
383: PetscLogInfo(ksp,"Linear solver has converged. Residual norm %g is less than absolute tolerance %g at iteration %dn",rnorm,ksp->atol,n);
384: *reason = KSP_CONVERGED_ATOL;
385: } else {
386: PetscLogInfo(ksp,"Linear solver has converged. Residual norm %g is less than relative tolerance %g times initial residual norm %g at iteration %dn",rnorm,ksp->rtol,ksp->rnorm0,n);
387: *reason = KSP_CONVERGED_RTOL;
388: }
389: } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
390: PetscLogInfo(ksp,"Linear solver is diverging. Initial residual norm %g, current residual norm %g at iteration %dn",ksp->rnorm0,rnorm,n);
391: *reason = KSP_DIVERGED_DTOL;
392: } else if (rnorm != rnorm) {
393: PetscLogInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergencen");
394: *reason = KSP_DIVERGED_DTOL;
395: }
396: return(0);
397: }
399: #undef __FUNCT__
401: /*
402: KSPDefaultBuildSolution - Default code to create/move the solution.
404: Input Parameters:
405: + ksp - iterative context
406: - v - pointer to the user's vector
408: Output Parameter:
409: . V - pointer to a vector containing the solution
411: Level: advanced
413: .keywords: KSP, build, solution, default
415: .seealso: KSPGetSolution(), KSPDefaultBuildResidual()
416: */
417: int KSPDefaultBuildSolution(KSP ksp,Vec v,Vec *V)
418: {
421: if (ksp->pc_side == PC_RIGHT) {
422: if (ksp->B) {
423: if (v) {KSP_PCApply(ksp,ksp->B,ksp->vec_sol,v); *V = v;}
424: else {SETERRQ(PETSC_ERR_SUP,"Not working with right preconditioner");}
425: } else {
426: if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
427: else { *V = ksp->vec_sol;}
428: }
429: } else if (ksp->pc_side == PC_SYMMETRIC) {
430: if (ksp->B) {
431: if (ksp->transpose_solve) SETERRQ(PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
432: if (v) {PCApplySymmetricRight(ksp->B,ksp->vec_sol,v); *V = v;}
433: else {SETERRQ(PETSC_ERR_SUP,"Not working with symmetric preconditioner");}
434: } else {
435: if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
436: else { *V = ksp->vec_sol;}
437: }
438: } else {
439: if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
440: else { *V = ksp->vec_sol; }
441: }
442: return(0);
443: }
445: #undef __FUNCT__
447: /*
448: KSPDefaultBuildResidual - Default code to compute the residual.
450: Input Parameters:
451: . ksp - iterative context
452: . t - pointer to temporary vector
453: . v - pointer to user vector
455: Output Parameter:
456: . V - pointer to a vector containing the residual
458: Level: advanced
460: .keywords: KSP, build, residual, default
462: .seealso: KSPDefaultBuildSolution()
463: */
464: int KSPDefaultBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
465: {
466: int ierr;
467: MatStructure pflag;
468: Vec T;
469: PetscScalar mone = -1.0;
470: Mat Amat,Pmat;
473: PCGetOperators(ksp->B,&Amat,&Pmat,&pflag);
474: KSPBuildSolution(ksp,t,&T);
475: KSP_MatMult(ksp,Amat,t,v);
476: VecAYPX(&mone,ksp->vec_rhs,v);
477: *V = v;
478: return(0);
479: }
481: #undef __FUNCT__
483: /*
484: KSPDefaultGetWork - Gets a number of work vectors.
486: Input Parameters:
487: . ksp - iterative context
488: . nw - number of work vectors to allocate
490: Notes:
491: Call this only if no work vectors have been allocated
492: */
493: int KSPDefaultGetWork(KSP ksp,int nw)
494: {
498: if (ksp->work) {KSPDefaultFreeWork(ksp);}
499: ksp->nwork = nw;
500: VecDuplicateVecs(ksp->vec_rhs,nw,&ksp->work);
501: PetscLogObjectParents(ksp,nw,ksp->work);
502: return(0);
503: }
505: #undef __FUNCT__
507: /*
508: KSPDefaultDestroy - Destroys a iterative context variable for methods with
509: no separate context. Preferred calling sequence KSPDestroy().
511: Input Parameter:
512: . ksp - the iterative context
513: */
514: int KSPDefaultDestroy(KSP ksp)
515: {
520: if (ksp->data) {PetscFree(ksp->data);}
522: /* free work vectors */
523: KSPDefaultFreeWork(ksp);
524: return(0);
525: }
527: #undef __FUNCT__
529: /*@C
530: KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.
532: Not Collective
534: Input Parameter:
535: . ksp - the KSP context
537: Output Parameter:
538: . reason - negative value indicates diverged, positive value converged, see KSPConvergedReason
540: Possible values for reason:
541: + KSP_CONVERGED_RTOL (residual norm decreased by a factor of rtol)
542: . KSP_CONVERGED_ATOL (residual norm less than atol)
543: . KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration)
544: . KSP_CONVERGED_QCG_NEG_CURVE
545: . KSP_CONVERGED_QCG_CONSTRAINED
546: . KSP_CONVERGED_STEP_LENGTH
547: . KSP_DIVERGED_ITS (required more than its to reach convergence)
548: . KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
549: . KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
550: - KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial
551: residual. Try a different preconditioner, or a different initial guess.
552:
554: Level: intermediate
556: Notes: Can only be called after the call the KSPSolve() is complete.
558: .keywords: KSP, nonlinear, set, convergence, test
560: .seealso: KSPSetConvergenceTest(), KSPDefaultConverged(), KSPSetTolerances(), KSPConvergedReason
561: @*/
562: int KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
563: {
566: *reason = ksp->reason;
567: return(0);
568: }