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