Actual source code: iterativ.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  1: /*
  2:    This file contains some simple default routines.
  3:    These routines should be SHORT, since they will be included in every
  4:    executable image that uses the iterative routines (note that, through
  5:    the registry system, we provide a way to load only the truely necessary
  6:    files)
  7:  */
  8: #include <petsc-private/kspimpl.h>   /*I "petscksp.h" I*/
  9: #include <petscdmshell.h>

 13: /*@
 14:    KSPGetResidualNorm - Gets the last (approximate preconditioned)
 15:    residual norm that has been computed.

 17:    Not Collective

 19:    Input Parameters:
 20: .  ksp - the iterative context

 22:    Output Parameters:
 23: .  rnorm - residual norm

 25:    Level: intermediate

 27: .keywords: KSP, get, residual norm

 29: .seealso: KSPBuildResidual()
 30: @*/
 31: PetscErrorCode  KSPGetResidualNorm(KSP ksp,PetscReal *rnorm)
 32: {
 36:   *rnorm = ksp->rnorm;
 37:   return(0);
 38: }

 42: /*@
 43:    KSPGetIterationNumber - Gets the current iteration number; if the
 44:          KSPSolve() is complete, returns the number of iterations
 45:          used.

 47:    Not Collective

 49:    Input Parameters:
 50: .  ksp - the iterative context

 52:    Output Parameters:
 53: .  its - number of iterations

 55:    Level: intermediate

 57:    Notes:
 58:       During the ith iteration this returns i-1
 59: .keywords: KSP, get, residual norm

 61: .seealso: KSPBuildResidual(), KSPGetResidualNorm()
 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.

 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;
104:   PetscViewer    viewer = (PetscViewer) dummy;

108:   if (!viewer) {
109:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
110:   }
111:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
112:   if (!ksp->calc_sings) {
113:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
114:   } else {
115:     KSPComputeExtremeSingularValues(ksp,&emax,&emin);
116:     c    = emax/emin;
117:     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);
118:   }
119:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
120:   return(0);
121: }

125: /*@C
126:    KSPMonitorSolution - Monitors progress of the KSP solvers by calling
127:    VecView() for the approximate solution at each iteration.

129:    Collective on KSP

131:    Input Parameters:
132: +  ksp - the KSP context
133: .  its - iteration number
134: .  fgnorm - 2-norm of residual (or gradient)
135: -  dummy - either a viewer or NULL

137:    Level: intermediate

139:    Notes:
140:     For some Krylov methods such as GMRES constructing the solution at
141:   each iteration is expensive, hence using this will slow the code.

143: .keywords: KSP, nonlinear, vector, monitor, view

145: .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView()
146: @*/
147: PetscErrorCode  KSPMonitorSolution(KSP ksp,PetscInt its,PetscReal fgnorm,void *dummy)
148: {
150:   Vec            x;
151:   PetscViewer    viewer = (PetscViewer) dummy;

154:   KSPBuildSolution(ksp,NULL,&x);
155:   if (!viewer) {
156:     MPI_Comm comm;
157:     PetscObjectGetComm((PetscObject)ksp,&comm);
158:     viewer = PETSC_VIEWER_DRAW_(comm);
159:   }
160:   VecView(x,viewer);
161:   return(0);
162: }

166: /*@C
167:    KSPMonitorDefault - Print the residual norm at each iteration of an
168:    iterative solver.

170:    Collective on KSP

172:    Input Parameters:
173: +  ksp   - iterative context
174: .  n     - iteration number
175: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
176: -  dummy - unused monitor context

178:    Level: intermediate

180: .keywords: KSP, default, monitor, residual

182: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
183: @*/
184: PetscErrorCode  KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
185: {
187:   PetscViewer    viewer = (PetscViewer) dummy;

190:   if (!viewer) {
191:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
192:   }
193:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
194:   if (n == 0 && ((PetscObject)ksp)->prefix) {
195:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
196:   }
197:   PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
198:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
199:   return(0);
200: }

204: /*@C
205:    KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
206:    residual norm at each iteration of an iterative solver.

208:    Collective on KSP

210:    Input Parameters:
211: +  ksp   - iterative context
212: .  n     - iteration number
213: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
214: -  dummy - unused monitor context

216:    Options Database Key:
217: .  -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()

219:    Notes:
220:    When using right preconditioning, these values are equivalent.

222:    Level: intermediate

224: .keywords: KSP, default, monitor, residual

226: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
227: @*/
228: PetscErrorCode  KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
229: {
231:   Vec            resid;
232:   PetscReal      truenorm,bnorm;
233:   PetscViewer    viewer = (PetscViewer)dummy;
234:   char           normtype[256];

237:   if (!viewer) {
238:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
239:   }
240:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
241:   if (n == 0 && ((PetscObject)ksp)->prefix) {
242:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
243:   }
244:   KSPBuildResidual(ksp,NULL,NULL,&resid);
245:   VecNorm(resid,NORM_2,&truenorm);
246:   VecDestroy(&resid);
247:   VecNorm(ksp->vec_rhs,NORM_2,&bnorm);
248:   PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
249:   PetscStrtolower(normtype);
250:   PetscViewerASCIIPrintf(viewer,"%3D KSP %s resid norm %14.12e true resid norm %14.12e ||r(i)||/||b|| %14.12e\n",n,normtype,(double)rnorm,(double)truenorm,(double)(truenorm/bnorm));
251:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
252:   return(0);
253: }

257: /*@C
258:    KSPMonitorTrueResidualMaxNorm - Prints the true residual max norm as well as the preconditioned
259:    residual norm at each iteration of an iterative solver.

261:    Collective on KSP

263:    Input Parameters:
264: +  ksp   - iterative context
265: .  n     - iteration number
266: .  rnorm - norm (preconditioned) residual value (may be estimated).
267: -  dummy - unused monitor context

269:    Options Database Key:
270: .  -ksp_monitor_max - Activates KSPMonitorTrueResidualMaxNorm()

272:    Notes:
273:    This could be implemented (better) with a flag in ksp.

275:    Level: intermediate

277: .keywords: KSP, default, monitor, residual

279: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualNorm()
280: @*/
281: PetscErrorCode  KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
282: {
284:   Vec            resid;
285:   PetscReal      truenorm,bnorm;
286:   PetscViewer    viewer = (PetscViewer)dummy;
287:   char           normtype[256];

290:   if (!viewer) {
291:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
292:   }
293:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
294:   if (n == 0 && ((PetscObject)ksp)->prefix) {
295:     PetscViewerASCIIPrintf(viewer,"  Residual norms (max) for %s solve.\n",((PetscObject)ksp)->prefix);
296:   }
297:   KSPBuildResidual(ksp,NULL,NULL,&resid);
298:   VecNorm(resid,NORM_INFINITY,&truenorm);
299:   VecDestroy(&resid);
300:   VecNorm(ksp->vec_rhs,NORM_INFINITY,&bnorm);
301:   PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
302:   PetscStrtolower(normtype);
303:   /* PetscViewerASCIIPrintf(viewer,"%3D KSP %s resid norm %14.12e true resid norm %14.12e ||r(i)||_inf/||b||_inf %14.12e\n",n,normtype,(double)rnorm,(double)truenorm,(double)(truenorm/bnorm)); */
304:   PetscViewerASCIIPrintf(viewer,"%3D KSP true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n",n,(double)truenorm,(double)(truenorm/bnorm));
305:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
306:   return(0);
307: }

311: PetscErrorCode  KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
312: {
314:   Vec            resid;
315:   PetscReal      rmax,pwork;
316:   PetscInt       i,n,N;
317:   const PetscScalar *r;

320:   KSPBuildResidual(ksp,NULL,NULL,&resid);
321:   VecNorm(resid,NORM_INFINITY,&rmax);
322:   VecGetLocalSize(resid,&n);
323:   VecGetSize(resid,&N);
324:   VecGetArrayRead(resid,&r);
325:   pwork = 0.0;
326:   for (i=0; i<n; i++) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
327:   VecRestoreArrayRead(resid,&r);
328:   VecDestroy(&resid);
329:   MPI_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
330:   *per = *per/N;
331:   return(0);
332: }

336: /*@C
337:    KSPMonitorRange - Prints the percentage of residual elements that are more then 10 percent of the maximum value.

339:    Collective on KSP

341:    Input Parameters:
342: +  ksp   - iterative context
343: .  it    - iteration number
344: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
345: -  dummy - unused monitor context

347:    Options Database Key:
348: .  -ksp_monitor_range - Activates KSPMonitorRange()

350:    Level: intermediate

352: .keywords: KSP, default, monitor, residual

354: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
355: @*/
356: PetscErrorCode  KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,void *dummy)
357: {
358:   PetscErrorCode   ierr;
359:   PetscReal        perc,rel;
360:   PetscViewer      viewer = (PetscViewer)dummy;
361:   /* should be in a MonitorRangeContext */
362:   static PetscReal prev;

365:   if (!viewer) {
366:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
367:   }
368:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
369:   if (!it) prev = rnorm;
370:   if (it == 0 && ((PetscObject)ksp)->prefix) {
371:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
372:   }
373:   KSPMonitorRange_Private(ksp,it,&perc);

375:   rel  = (prev - rnorm)/prev;
376:   prev = rnorm;
377:   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));
378:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
379:   return(0);
380: }

384: /*@C
385:    KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
386:    outer iteration in an adaptive way.

388:    Collective on KSP

390:    Input Parameters:
391: +  ksp   - iterative context
392: .  n     - iteration number (not used)
393: .  fnorm - the current residual norm
394: .  dummy - some context as a C struct. fields:
395:              coef: a scaling coefficient. default 1.0. can be passed through
396:                    -sub_ksp_dynamic_tolerance_param
397:              bnrm: norm of the right-hand side. store it to avoid repeated calculation

399:    Notes:
400:    This may be useful for a flexibly preconditioner Krylov method to
401:    control the accuracy of the inner solves needed to gaurantee the
402:    convergence of the outer iterations.

404:    Level: advanced

406: .keywords: KSP, inner tolerance

408: .seealso: KSPMonitorDynamicToleranceDestroy()
409: @*/
410: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
411: {
413:   PC             pc;
414:   PetscReal      outer_rtol, outer_abstol, outer_dtol, inner_rtol;
415:   PetscInt       outer_maxits,nksp,first,i;
416:   KSPDynTolCtx   *scale   = (KSPDynTolCtx*)dummy;
417:   KSP            kspinner = NULL, *subksp = NULL;

420:   KSPGetPC(ksp, &pc);

422:   /* compute inner_rtol */
423:   if (scale->bnrm < 0.0) {
424:     Vec b;
425:     KSPGetRhs(ksp, &b);
426:     VecNorm(b, NORM_2, &(scale->bnrm));
427:   }
428:   KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
429:   inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
430:   /*PetscPrintf(PETSC_COMM_WORLD, "        Inner rtol = %g\n", (double)inner_rtol);*/

432:   /* if pc is ksp */
433:   PCKSPGetKSP(pc, &kspinner);
434:   if (kspinner) {
435:     KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
436:     return(0);
437:   }

439:   /* if pc is bjacobi */
440:   PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
441:   if (subksp) {
442:     for (i=0; i<nksp; i++) {
443:       KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
444:     }
445:     return(0);
446:   }

448:   /* todo: dynamic tolerance may apply to other types of pc too */
449:   return(0);
450: }

454: /*
455:   Destroy the dummy context used in KSPMonitorDynamicTolerance()
456: */
457: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
458: {

462:   PetscFree(*dummy);
463:   return(0);
464: }

468: /*
469:   Default (short) KSP Monitor, same as KSPMonitorDefault() except
470:   it prints fewer digits of the residual as the residual gets smaller.
471:   This is because the later digits are meaningless and are often
472:   different on different machines; by using this routine different
473:   machines will usually generate the same output.
474: */
475: PetscErrorCode  KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
476: {
478:   PetscViewer    viewer = (PetscViewer)dummy;

481:   if (!viewer) {
482:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
483:   }
484:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
485:   if (its == 0 && ((PetscObject)ksp)->prefix) {
486:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
487:   }

489:   if (fnorm > 1.e-9) {
490:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %g \n",its,(double)fnorm);
491:   } else if (fnorm > 1.e-11) {
492:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,(double)fnorm);
493:   } else {
494:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);
495:   }
496:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
497:   return(0);
498: }

502: /*@C
503:    KSPConvergedSkip - Convergence test that do not return as converged
504:    until the maximum number of iterations is reached.

506:    Collective on KSP

508:    Input Parameters:
509: +  ksp   - iterative context
510: .  n     - iteration number
511: .  rnorm - 2-norm residual value (may be estimated)
512: -  dummy - unused convergence context

514:    Returns:
515: .  reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS

517:    Notes:
518:    This should be used as the convergence test with the option
519:    KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
520:    not computed. Convergence is then declared after the maximum number
521:    of iterations have been reached. Useful when one is using CG or
522:    BiCGStab as a smoother.

524:    Level: advanced

526: .keywords: KSP, default, convergence, residual

528: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
529: @*/
530: PetscErrorCode  KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
531: {
535:   *reason = KSP_CONVERGED_ITERATING;
536:   if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
537:   return(0);
538: }


543: /*@C
544:    KSPConvergedDefaultCreate - Creates and initializes the space used by the KSPConvergedDefault() function context

546:    Collective on KSP

548:    Output Parameter:
549: .  ctx - convergence context

551:    Level: intermediate

553: .keywords: KSP, default, convergence, residual

555: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
556:           KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
557: @*/
558: PetscErrorCode  KSPConvergedDefaultCreate(void **ctx)
559: {
560:   PetscErrorCode         ierr;
561:   KSPConvergedDefaultCtx *cctx;

564:   PetscNew(&cctx);
565:   *ctx = cctx;
566:   return(0);
567: }

571: /*@
572:    KSPConvergedDefaultSetUIRNorm - makes the default convergence test use || B*(b - A*(initial guess))||
573:       instead of || B*b ||. In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
574:       is used there is no B in the above formula. UIRNorm is short for Use Initial Residual Norm.

576:    Collective on KSP

578:    Input Parameters:
579: .  ksp   - iterative context

581:    Options Database:
582: .   -ksp_converged_use_initial_residual_norm

584:    Notes:
585:    Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.

587:    The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
588:    are defined in petscksp.h.

590:    If the convergence test is not KSPConvergedDefault() then this is ignored.

592:    If right preconditioning is being used then B does not appear in the above formula.


595:    Level: intermediate

597: .keywords: KSP, default, convergence, residual

599: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm()
600: @*/
601: PetscErrorCode  KSPConvergedDefaultSetUIRNorm(KSP ksp)
602: {
603:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;

607:   if (ksp->converged != KSPConvergedDefault) return(0);
608:   if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
609:   ctx->initialrtol = PETSC_TRUE;
610:   return(0);
611: }

615: /*@
616:    KSPConvergedDefaultSetUMIRNorm - makes the default convergence test use min(|| B*(b - A*(initial guess))||,|| B*b ||)
617:       In the case of right preconditioner or if KSPSetNormType(ksp,KSP_NORM_UNPRECONDIITONED)
618:       is used there is no B in the above formula. UMIRNorm is short for Use Minimum Initial Residual Norm.

620:    Collective on KSP

622:    Input Parameters:
623: .  ksp   - iterative context

625:    Options Database:
626: .   -ksp_converged_use_min_initial_residual_norm

628:    Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.

630:    The precise values of reason are macros such as KSP_CONVERGED_RTOL, which
631:    are defined in petscksp.h.

633:    Level: intermediate

635: .keywords: KSP, default, convergence, residual

637: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm()
638: @*/
639: PetscErrorCode  KSPConvergedDefaultSetUMIRNorm(KSP ksp)
640: {
641:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;

645:   if (ksp->converged != KSPConvergedDefault) return(0);
646:   if (ctx->initialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
647:   ctx->mininitialrtol = PETSC_TRUE;
648:   return(0);
649: }

653: /*@C
654:    KSPConvergedDefault - Determines convergence of the linear iterative solvers by default

656:    Collective on KSP

658:    Input Parameters:
659: +  ksp   - iterative context
660: .  n     - iteration number
661: .  rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
662: -  ctx - convergence context which must be created by KSPConvergedDefaultCreate()

664:    Output Parameter:
665: +   positive - if the iteration has converged;
666: .   negative - if residual norm exceeds divergence threshold;
667: -   0 - otherwise.

669:    Notes:
670:    KSPConvergedDefault() reaches convergence when   rnorm < MAX (rtol * rnorm_0, abstol);
671:    Divergence is detected if  rnorm > dtol * rnorm_0,

673:    where:
674: +     rtol = relative tolerance,
675: .     abstol = absolute tolerance.
676: .     dtol = divergence tolerance,
677: -     rnorm_0 is the two norm of the right hand side. When initial guess is non-zero you
678:           can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
679:           as the starting point for relative norm convergence testing, that is as rnorm_0

681:    Use KSPSetTolerances() to alter the defaults for rtol, abstol, dtol.

683:    Use KSPSetNormType() (or -ksp_norm_type <none,preconditioned,unpreconditioned,natural>) to change the norm used for computing rnorm

685:    The precise values of reason are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.

687:    This routine is used by KSP by default so the user generally never needs call it directly.

689:    Use KSPSetConvergenceTest() to provide your own test instead of using this one.

691:    Level: intermediate

693: .keywords: KSP, default, convergence, residual

695: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
696:           KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
697: @*/
698: PetscErrorCode  KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
699: {
700:   PetscErrorCode         ierr;
701:   KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
702:   KSPNormType            normtype;

707:   *reason = KSP_CONVERGED_ITERATING;

709:   KSPGetNormType(ksp,&normtype);
710:   if (normtype == KSP_NORM_NONE) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Use KSPConvergedSkip() with KSPNormType of KSP_NORM_NONE");

712:   if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
713:   if (!n) {
714:     /* if user gives initial guess need to compute norm of b */
715:     if (!ksp->guess_zero && !cctx->initialrtol) {
716:       PetscReal snorm;
717:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
718:         PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
719:         VecNorm(ksp->vec_rhs,NORM_2,&snorm);        /*     <- b'*b */
720:       } else {
721:         Vec z;
722:         /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
723:         VecDuplicate(ksp->vec_rhs,&z);
724:         KSP_PCApply(ksp,ksp->vec_rhs,z);
725:         if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
726:           PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
727:           VecNorm(z,NORM_2,&snorm);                 /*    dp <- b'*B'*B*b */
728:         } else if (ksp->normtype == KSP_NORM_NATURAL) {
729:           PetscScalar norm;
730:           PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
731:           VecDot(ksp->vec_rhs,z,&norm);
732:           snorm = PetscSqrtReal(PetscAbsScalar(norm));                            /*    dp <- b'*B*b */
733:         }
734:         VecDestroy(&z);
735:       }
736:       /* handle special case of zero RHS and nonzero guess */
737:       if (!snorm) {
738:         PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
739:         snorm = rnorm;
740:       }
741:       if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
742:       else ksp->rnorm0 = snorm;
743:     } else {
744:       ksp->rnorm0 = rnorm;
745:     }
746:     ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
747:   }

749:   if (n <= ksp->chknorm) return(0);

751:   if (PetscIsInfOrNanScalar(rnorm)) {
752:     PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
753:     *reason = KSP_DIVERGED_NANORINF;
754:   } else if (rnorm <= ksp->ttol) {
755:     if (rnorm < ksp->abstol) {
756:       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);
757:       *reason = KSP_CONVERGED_ATOL;
758:     } else {
759:       if (cctx->initialrtol) {
760:         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);
761:       } else {
762:         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);
763:       }
764:       *reason = KSP_CONVERGED_RTOL;
765:     }
766:   } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
767:     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);
768:     *reason = KSP_DIVERGED_DTOL;
769:   }
770:   return(0);
771: }

775: /*@C
776:    KSPConvergedDefaultDestroy - Frees the space used by the KSPConvergedDefault() function context

778:    Collective on KSP

780:    Input Parameters:
781: .  ctx - convergence context

783:    Level: intermediate

785: .keywords: KSP, default, convergence, residual

787: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
788:           KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
789: @*/
790: PetscErrorCode  KSPConvergedDefaultDestroy(void *ctx)
791: {
792:   PetscErrorCode         ierr;
793:   KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;

796:   VecDestroy(&cctx->work);
797:   PetscFree(ctx);
798:   return(0);
799: }

803: /*
804:    KSPBuildSolutionDefault - Default code to create/move the solution.

806:    Input Parameters:
807: +  ksp - iterative context
808: -  v   - pointer to the user's vector

810:    Output Parameter:
811: .  V - pointer to a vector containing the solution

813:    Level: advanced

815:    Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations

817: .keywords:  KSP, build, solution, default

819: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
820: */
821: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
822: {

826:   if (ksp->pc_side == PC_RIGHT) {
827:     if (ksp->pc) {
828:       if (v) {
829:         KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
830:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
831:     } else {
832:       if (v) {
833:         VecCopy(ksp->vec_sol,v); *V = v;
834:       } else *V = ksp->vec_sol;
835:     }
836:   } else if (ksp->pc_side == PC_SYMMETRIC) {
837:     if (ksp->pc) {
838:       if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
839:       if (v) {
840:         PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
841:         *V = v;
842:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
843:     } else {
844:       if (v) {
845:         VecCopy(ksp->vec_sol,v); *V = v;
846:       } else *V = ksp->vec_sol;
847:     }
848:   } else {
849:     if (v) {
850:       VecCopy(ksp->vec_sol,v); *V = v;
851:     } else *V = ksp->vec_sol;
852:   }
853:   return(0);
854: }

858: /*
859:    KSPBuildResidualDefault - Default code to compute the residual.

861:    Input Parameters:
862: .  ksp - iterative context
863: .  t   - pointer to temporary vector
864: .  v   - pointer to user vector

866:    Output Parameter:
867: .  V - pointer to a vector containing the residual

869:    Level: advanced

871:    Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations

873: .keywords:  KSP, build, residual, default

875: .seealso: KSPBuildSolutionDefault()
876: */
877: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
878: {
880:   Mat            Amat,Pmat;

883:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
884:   PCGetOperators(ksp->pc,&Amat,&Pmat);
885:   KSPBuildSolution(ksp,t,NULL);
886:   KSP_MatMult(ksp,Amat,t,v);
887:   VecAYPX(v,-1.0,ksp->vec_rhs);
888:   *V   = v;
889:   return(0);
890: }

894: /*@C
895:   KSPGetVecs - Gets a number of work vectors.

897:   Input Parameters:
898: + ksp  - iterative context
899: . rightn  - number of right work vectors
900: - leftn   - number of left work vectors to allocate

902:   Output Parameter:
903: +  right - the array of vectors created
904: -  left - the array of left vectors

906:    Note: The right vector has as many elements as the matrix has columns. The left
907:      vector has as many elements as the matrix has rows.

909:    Level: advanced

911: .seealso:   MatGetVecs()

913: @*/
914: PetscErrorCode KSPGetVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
915: {
917:   Vec            vecr,vecl;

920:   if (rightn) {
921:     if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
922:     if (ksp->vec_sol) vecr = ksp->vec_sol;
923:     else {
924:       0;
925:       if (ksp->dm) {
926:         DMGetGlobalVector(ksp->dm,&vecr); /* don't check for errors -- if any errors, pass down to next block */
927:       }
928:       if (!ksp->dm || ierr) {
929:         Mat mat;
930:         if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
931:         PCGetOperators(ksp->pc,&mat,NULL);
932:         MatGetVecs(mat,&vecr,NULL);
933:       }
934:     }
935:     VecDuplicateVecs(vecr,rightn,right);
936:     if (!ksp->vec_sol) {
937:       if (ksp->dm) {
938:         DMRestoreGlobalVector(ksp->dm,&vecr);
939:       } else {
940:         VecDestroy(&vecr);
941:       }
942:     }
943:   }
944:   if (leftn) {
945:     if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
946:     if (ksp->vec_rhs) vecl = ksp->vec_rhs;
947:     else {
948:       0;
949:       if (ksp->dm) {
950:         DMGetGlobalVector(ksp->dm,&vecl); /* don't check for errors -- if any errors, pass down to next block */
951:       }
952:       if (!ksp->dm || ierr) {
953:         Mat mat;
954:         if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
955:         PCGetOperators(ksp->pc,&mat,NULL);
956:         MatGetVecs(mat,NULL,&vecl);
957:       }
958:     }
959:     VecDuplicateVecs(vecl,leftn,left);
960:     if (!ksp->vec_rhs) {
961:       if (ksp->dm) {
962:         DMRestoreGlobalVector(ksp->dm,&vecl);
963:       } else {
964:         VecDestroy(&vecl);
965:       }
966:     }
967:   }
968:   return(0);
969: }

973: /*
974:   KSPSetWorkVecs - Sets a number of work vectors into a KSP object

976:   Input Parameters:
977: . ksp  - iterative context
978: . nw   - number of work vectors to allocate

980:    Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations

982:  */
983: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
984: {

988:   VecDestroyVecs(ksp->nwork,&ksp->work);
989:   ksp->nwork = nw;
990:   KSPGetVecs(ksp,nw,&ksp->work,0,NULL);
991:   PetscLogObjectParents(ksp,nw,ksp->work);
992:   return(0);
993: }

997: /*
998:   KSPDestroyDefault - Destroys a iterative context variable for methods with
999:   no separate context.  Preferred calling sequence KSPDestroy().

1001:   Input Parameter:
1002: . ksp - the iterative context

1004:    Developers Note: This is PETSC_EXTERN because it may be used by user written plugin KSP implementations

1006: */
1007: PetscErrorCode KSPDestroyDefault(KSP ksp)
1008: {

1013:   PetscFree(ksp->data);
1014:   return(0);
1015: }

1019: /*@
1020:    KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.

1022:    Not Collective

1024:    Input Parameter:
1025: .  ksp - the KSP context

1027:    Output Parameter:
1028: .  reason - negative value indicates diverged, positive value converged, see KSPConvergedReason

1030:    Possible values for reason:
1031: +  KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1032: .  KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1033: .  KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence
1034:            test routine is set.
1035: .  KSP_CONVERGED_CG_NEG_CURVE
1036: .  KSP_CONVERGED_CG_CONSTRAINED
1037: .  KSP_CONVERGED_STEP_LENGTH
1038: .  KSP_DIVERGED_ITS  (required more than its to reach convergence)
1039: .  KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1040: .  KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1041: .  KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1042: -  KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial
1043:                                 residual. Try a different preconditioner, or a different initial Level.)

1045:    See also manual page for each reason.

1047:    guess: beginner

1049:    Notes: Can only be called after the call the KSPSolve() is complete.

1051:    Level: intermediate

1053: .keywords: KSP, nonlinear, set, convergence, test

1055: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason
1056: @*/
1057: PetscErrorCode  KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1058: {
1062:   *reason = ksp->reason;
1063:   return(0);
1064: }

1066: #include <petsc-private/dmimpl.h>
1069: /*@
1070:    KSPSetDM - Sets the DM that may be used by some preconditioners

1072:    Logically Collective on KSP

1074:    Input Parameters:
1075: +  ksp - the preconditioner context
1076: -  dm - the dm

1078:    Notes: If this is used then the KSP will attempt to use the DM to create the matrix and use the routine
1079:           set with DMKSPSetComputeOperators(). Use KSPSetDMActive(ksp,PETSC_FALSE) to instead use the matrix
1080:           you've provided with KSPSetOperators().

1082:    Level: intermediate

1084: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1085: @*/
1086: PetscErrorCode  KSPSetDM(KSP ksp,DM dm)
1087: {
1089:   PC             pc;

1093:   if (dm) {PetscObjectReference((PetscObject)dm);}
1094:   if (ksp->dm) {                /* Move the DMSNES context over to the new DM unless the new DM already has one */
1095:     if (ksp->dm->dmksp && ksp->dmAuto && !dm->dmksp) {
1096:       DMKSP kdm;
1097:       DMCopyDMKSP(ksp->dm,dm);
1098:       DMGetDMKSP(ksp->dm,&kdm);
1099:       if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1100:     }
1101:     DMDestroy(&ksp->dm);
1102:   }
1103:   ksp->dm       = dm;
1104:   ksp->dmAuto   = PETSC_FALSE;
1105:   KSPGetPC(ksp,&pc);
1106:   PCSetDM(pc,dm);
1107:   ksp->dmActive = PETSC_TRUE;
1108:   return(0);
1109: }

1113: /*@
1114:    KSPSetDMActive - Indicates the DM should be used to generate the linear system matrix and right hand side

1116:    Logically Collective on KSP

1118:    Input Parameters:
1119: +  ksp - the preconditioner context
1120: -  flg - use the DM

1122:    Notes:
1123:    By default KSPSetDM() sets the DM as active, call KSPSetDMActive(ksp,PETSC_FALSE); after KSPSetDM(ksp,dm) to not have the KSP object use the DM to generate the matrices.

1125:    Level: intermediate

1127: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1128: @*/
1129: PetscErrorCode  KSPSetDMActive(KSP ksp,PetscBool flg)
1130: {
1134:   ksp->dmActive = flg;
1135:   return(0);
1136: }

1140: /*@
1141:    KSPGetDM - Gets the DM that may be used by some preconditioners

1143:    Not Collective

1145:    Input Parameter:
1146: . ksp - the preconditioner context

1148:    Output Parameter:
1149: .  dm - the dm

1151:    Level: intermediate


1154: .seealso: KSPSetDM(), KSPSetDMActive()
1155: @*/
1156: PetscErrorCode  KSPGetDM(KSP ksp,DM *dm)
1157: {

1162:   if (!ksp->dm) {
1163:     DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1164:     ksp->dmAuto = PETSC_TRUE;
1165:   }
1166:   *dm = ksp->dm;
1167:   return(0);
1168: }

1172: /*@
1173:    KSPSetApplicationContext - Sets the optional user-defined context for the linear solver.

1175:    Logically Collective on KSP

1177:    Input Parameters:
1178: +  ksp - the KSP context
1179: -  usrP - optional user context

1181:    Level: intermediate

1183: .keywords: KSP, set, application, context

1185: .seealso: KSPGetApplicationContext()
1186: @*/
1187: PetscErrorCode  KSPSetApplicationContext(KSP ksp,void *usrP)
1188: {
1190:   PC             pc;

1194:   ksp->user = usrP;
1195:   KSPGetPC(ksp,&pc);
1196:   PCSetApplicationContext(pc,usrP);
1197:   return(0);
1198: }

1202: /*@
1203:    KSPGetApplicationContext - Gets the user-defined context for the linear solver.

1205:    Not Collective

1207:    Input Parameter:
1208: .  ksp - KSP context

1210:    Output Parameter:
1211: .  usrP - user context

1213:    Level: intermediate

1215: .keywords: KSP, get, application, context

1217: .seealso: KSPSetApplicationContext()
1218: @*/
1219: PetscErrorCode  KSPGetApplicationContext(KSP ksp,void *usrP)
1220: {
1223:   *(void**)usrP = ksp->user;
1224:   return(0);
1225: }