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: }