Actual source code: iterativ.c

petsc-3.7.3 2016-07-24
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(), KSPGetTotalIterations()
 62: @*/
 63: PetscErrorCode  KSPGetIterationNumber(KSP ksp,PetscInt *its)
 64: {
 68:   *its = ksp->its;
 69:   return(0);
 70: }

 74: /*@
 75:    KSPGetTotalIterations - Gets the total number of iterations this KSP object has performed since was created, counted over all linear solves

 77:    Not Collective

 79:    Input Parameters:
 80: .  ksp - the iterative context

 82:    Output Parameters:
 83: .  its - total number of iterations

 85:    Level: intermediate

 87:    Notes: Use KSPGetIterationNumber() to get the count for the most recent solve only
 88:    If this is called within a linear solve (such as in a KSPMonitor routine) then it does not include iterations within that current solve

 90: .keywords: KSP, get, residual norm

 92: .seealso: KSPBuildResidual(), KSPGetResidualNorm(), KSPGetIterationNumber()
 93: @*/
 94: PetscErrorCode  KSPGetTotalIterations(KSP ksp,PetscInt *its)
 95: {
 99:   *its = ksp->totalits;
100:   return(0);
101: }

105: /*@C
106:     KSPMonitorSingularValue - Prints the two norm of the true residual and
107:     estimation of the extreme singular values of the preconditioned problem
108:     at each iteration.

110:     Logically Collective on KSP

112:     Input Parameters:
113: +   ksp - the iterative context
114: .   n  - the iteration
115: -   rnorm - the two norm of the residual

117:     Options Database Key:
118: .   -ksp_monitor_singular_value - Activates KSPMonitorSingularValue()

120:     Notes:
121:     The CG solver uses the Lanczos technique for eigenvalue computation,
122:     while GMRES uses the Arnoldi technique; other iterative methods do
123:     not currently compute singular values.

125:     Level: intermediate

127: .keywords: KSP, CG, default, monitor, extreme, singular values, Lanczos, Arnoldi

129: .seealso: KSPComputeExtremeSingularValues()
130: @*/
131: PetscErrorCode  KSPMonitorSingularValue(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
132: {
133:   PetscReal      emin,emax,c;
135:   PetscViewer    viewer = dummy->viewer;

140:   PetscViewerPushFormat(viewer,dummy->format);
141:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
142:   if (!ksp->calc_sings) {
143:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
144:   } else {
145:     KSPComputeExtremeSingularValues(ksp,&emax,&emin);
146:     c    = emax/emin;
147:     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);
148:   }
149:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
150:   PetscViewerPopFormat(viewer);
151:   return(0);
152: }

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

160:    Collective on KSP

162:    Input Parameters:
163: +  ksp - the KSP context
164: .  its - iteration number
165: .  fgnorm - 2-norm of residual (or gradient)
166: -  dummy - a viewer

168:    Level: intermediate

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

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

176: .seealso: KSPMonitorSet(), KSPMonitorDefault(), VecView()
177: @*/
178: PetscErrorCode  KSPMonitorSolution(KSP ksp,PetscInt its,PetscReal fgnorm,PetscViewerAndFormat *dummy)
179: {
181:   Vec            x;
182:   PetscViewer    viewer = dummy->viewer;

186:   KSPBuildSolution(ksp,NULL,&x);
187:   PetscViewerPushFormat(viewer,dummy->format);
188:   VecView(x,viewer);
189:   PetscViewerPopFormat(viewer);
190:   return(0);
191: }

195: /*@C
196:    KSPMonitorDefault - Print the residual norm at each iteration of an
197:    iterative solver.

199:    Collective on KSP

201:    Input Parameters:
202: +  ksp   - iterative context
203: .  n     - iteration number
204: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
205: -  dummy - an ASCII PetscViewer

207:    Level: intermediate

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

211: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
212: @*/
213: PetscErrorCode  KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
214: {
216:   PetscViewer    viewer =  dummy->viewer;

220:   PetscViewerPushFormat(viewer,dummy->format);
221:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
222:   if (n == 0 && ((PetscObject)ksp)->prefix) {
223:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
224:   }
225:   PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
226:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
227:   PetscViewerPopFormat(viewer);
228:   return(0);
229: }

233: /*@C
234:    KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
235:    residual norm at each iteration of an iterative solver.

237:    Collective on KSP

239:    Input Parameters:
240: +  ksp   - iterative context
241: .  n     - iteration number
242: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
243: -  dummy - an ASCII PetscViewer

245:    Options Database Key:
246: .  -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()

248:    Notes:
249:    When using right preconditioning, these values are equivalent.

251:    Level: intermediate

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

255: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
256: @*/
257: PetscErrorCode  KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
258: {
260:   Vec            resid;
261:   PetscReal      truenorm,bnorm;
262:   PetscViewer    viewer = dummy->viewer;
263:   char           normtype[256];

267:   PetscViewerPushFormat(viewer,dummy->format);
268:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
269:   if (n == 0 && ((PetscObject)ksp)->prefix) {
270:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
271:   }
272:   KSPBuildResidual(ksp,NULL,NULL,&resid);
273:   VecNorm(resid,NORM_2,&truenorm);
274:   VecDestroy(&resid);
275:   VecNorm(ksp->vec_rhs,NORM_2,&bnorm);
276:   PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
277:   PetscStrtolower(normtype);
278:   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));
279:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
280:   PetscViewerPopFormat(viewer);
281:   return(0);
282: }

286: /*@C
287:    KSPMonitorTrueResidualMaxNorm - Prints the true residual max norm each iteration of an iterative solver.

289:    Collective on KSP

291:    Input Parameters:
292: +  ksp   - iterative context
293: .  n     - iteration number
294: .  rnorm - norm (preconditioned) residual value (may be estimated).
295: -  dummy - an ASCII viewer

297:    Options Database Key:
298: .  -ksp_monitor_max - Activates KSPMonitorTrueResidualMaxNorm()

300:    Notes:
301:    This could be implemented (better) with a flag in ksp.

303:    Level: intermediate

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

307: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualNorm()
308: @*/
309: PetscErrorCode  KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,PetscViewerAndFormat *dummy)
310: {
312:   Vec            resid;
313:   PetscReal      truenorm,bnorm;
314:   PetscViewer    viewer = dummy->viewer;
315:   char           normtype[256];

319:   PetscViewerPushFormat(viewer,dummy->format);
320:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
321:   if (n == 0 && ((PetscObject)ksp)->prefix) {
322:     PetscViewerASCIIPrintf(viewer,"  Residual norms (max) for %s solve.\n",((PetscObject)ksp)->prefix);
323:   }
324:   KSPBuildResidual(ksp,NULL,NULL,&resid);
325:   VecNorm(resid,NORM_INFINITY,&truenorm);
326:   VecDestroy(&resid);
327:   VecNorm(ksp->vec_rhs,NORM_INFINITY,&bnorm);
328:   PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
329:   PetscStrtolower(normtype);
330:   PetscViewerASCIIPrintf(viewer,"%3D KSP true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n",n,(double)truenorm,(double)(truenorm/bnorm));
331:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
332:   PetscViewerPopFormat(viewer);
333:   return(0);
334: }

338: PetscErrorCode  KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
339: {
341:   Vec            resid;
342:   PetscReal      rmax,pwork;
343:   PetscInt       i,n,N;
344:   const PetscScalar *r;

347:   KSPBuildResidual(ksp,NULL,NULL,&resid);
348:   VecNorm(resid,NORM_INFINITY,&rmax);
349:   VecGetLocalSize(resid,&n);
350:   VecGetSize(resid,&N);
351:   VecGetArrayRead(resid,&r);
352:   pwork = 0.0;
353:   for (i=0; i<n; i++) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
354:   VecRestoreArrayRead(resid,&r);
355:   VecDestroy(&resid);
356:   MPIU_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
357:   *per = *per/N;
358:   return(0);
359: }

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

366:    Collective on KSP

368:    Input Parameters:
369: +  ksp   - iterative context
370: .  it    - iteration number
371: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
372: -  dummy - an ASCII viewer

374:    Options Database Key:
375: .  -ksp_monitor_range - Activates KSPMonitorRange()

377:    Level: intermediate

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

381: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
382: @*/
383: PetscErrorCode  KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,PetscViewerAndFormat *dummy)
384: {
385:   PetscErrorCode   ierr;
386:   PetscReal        perc,rel;
387:   PetscViewer      viewer = dummy->viewer;
388:   /* should be in a MonitorRangeContext */
389:   static PetscReal prev;

393:   PetscViewerPushFormat(viewer,dummy->format);
394:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
395:   if (!it) prev = rnorm;
396:   if (it == 0 && ((PetscObject)ksp)->prefix) {
397:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
398:   }
399:   KSPMonitorRange_Private(ksp,it,&perc);

401:   rel  = (prev - rnorm)/prev;
402:   prev = rnorm;
403:   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));
404:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
405:   PetscViewerPopFormat(viewer);
406:   return(0);
407: }

411: /*@C
412:    KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
413:    outer iteration in an adaptive way.

415:    Collective on KSP

417:    Input Parameters:
418: +  ksp   - iterative context
419: .  n     - iteration number (not used)
420: .  fnorm - the current residual norm
421: .  dummy - some context as a C struct. fields:
422:              coef: a scaling coefficient. default 1.0. can be passed through
423:                    -sub_ksp_dynamic_tolerance_param
424:              bnrm: norm of the right-hand side. store it to avoid repeated calculation

426:    Notes:
427:    This may be useful for a flexibly preconditioner Krylov method to
428:    control the accuracy of the inner solves needed to gaurantee the
429:    convergence of the outer iterations.

431:    Level: advanced

433: .keywords: KSP, inner tolerance

435: .seealso: KSPMonitorDynamicToleranceDestroy()
436: @*/
437: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
438: {
440:   PC             pc;
441:   PetscReal      outer_rtol, outer_abstol, outer_dtol, inner_rtol;
442:   PetscInt       outer_maxits,nksp,first,i;
443:   KSPDynTolCtx   *scale   = (KSPDynTolCtx*)dummy;
444:   KSP            kspinner = NULL, *subksp = NULL;

447:   KSPGetPC(ksp, &pc);

449:   /* compute inner_rtol */
450:   if (scale->bnrm < 0.0) {
451:     Vec b;
452:     KSPGetRhs(ksp, &b);
453:     VecNorm(b, NORM_2, &(scale->bnrm));
454:   }
455:   KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
456:   inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
457:   /*PetscPrintf(PETSC_COMM_WORLD, "        Inner rtol = %g\n", (double)inner_rtol);*/

459:   /* if pc is ksp */
460:   PCKSPGetKSP(pc, &kspinner);
461:   if (kspinner) {
462:     KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
463:     return(0);
464:   }

466:   /* if pc is bjacobi */
467:   PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
468:   if (subksp) {
469:     for (i=0; i<nksp; i++) {
470:       KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
471:     }
472:     return(0);
473:   }

475:   /* todo: dynamic tolerance may apply to other types of pc too */
476:   return(0);
477: }

481: /*
482:   Destroy the dummy context used in KSPMonitorDynamicTolerance()
483: */
484: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
485: {

489:   PetscFree(*dummy);
490:   return(0);
491: }

495: /*
496:   Default (short) KSP Monitor, same as KSPMonitorDefault() except
497:   it prints fewer digits of the residual as the residual gets smaller.
498:   This is because the later digits are meaningless and are often
499:   different on different machines; by using this routine different
500:   machines will usually generate the same output.
501: */
502: PetscErrorCode  KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,PetscViewerAndFormat *dummy)
503: {
505:   PetscViewer    viewer = dummy->viewer;

509:   PetscViewerPushFormat(viewer,dummy->format);
510:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
511:   if (its == 0 && ((PetscObject)ksp)->prefix) {
512:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
513:   }

515:   if (fnorm > 1.e-9) {
516:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %g \n",its,(double)fnorm);
517:   } else if (fnorm > 1.e-11) {
518:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,(double)fnorm);
519:   } else {
520:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);
521:   }
522:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
523:   PetscViewerPopFormat(viewer);
524:   return(0);
525: }

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

533:    Collective on KSP

535:    Input Parameters:
536: +  ksp   - iterative context
537: .  n     - iteration number
538: .  rnorm - 2-norm residual value (may be estimated)
539: -  dummy - unused convergence context

541:    Returns:
542: .  reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS

544:    Notes:
545:    This should be used as the convergence test with the option
546:    KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
547:    not computed. Convergence is then declared after the maximum number
548:    of iterations have been reached. Useful when one is using CG or
549:    BiCGStab as a smoother.

551:    Level: advanced

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

555: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
556: @*/
557: PetscErrorCode  KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
558: {
562:   *reason = KSP_CONVERGED_ITERATING;
563:   if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
564:   return(0);
565: }


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

573:    Collective on KSP

575:    Output Parameter:
576: .  ctx - convergence context

578:    Level: intermediate

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

582: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
583:           KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
584: @*/
585: PetscErrorCode  KSPConvergedDefaultCreate(void **ctx)
586: {
587:   PetscErrorCode         ierr;
588:   KSPConvergedDefaultCtx *cctx;

591:   PetscNew(&cctx);
592:   *ctx = cctx;
593:   return(0);
594: }

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

603:    Collective on KSP

605:    Input Parameters:
606: .  ksp   - iterative context

608:    Options Database:
609: .   -ksp_converged_use_initial_residual_norm

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

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

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

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


622:    Level: intermediate

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

626: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm()
627: @*/
628: PetscErrorCode  KSPConvergedDefaultSetUIRNorm(KSP ksp)
629: {
630:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;

634:   if (ksp->converged != KSPConvergedDefault) return(0);
635:   if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
636:   ctx->initialrtol = PETSC_TRUE;
637:   return(0);
638: }

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

647:    Collective on KSP

649:    Input Parameters:
650: .  ksp   - iterative context

652:    Options Database:
653: .   -ksp_converged_use_min_initial_residual_norm

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

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

660:    Level: intermediate

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

664: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm()
665: @*/
666: PetscErrorCode  KSPConvergedDefaultSetUMIRNorm(KSP ksp)
667: {
668:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;

672:   if (ksp->converged != KSPConvergedDefault) return(0);
673:   if (ctx->initialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
674:   ctx->mininitialrtol = PETSC_TRUE;
675:   return(0);
676: }

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

683:    Collective on KSP

685:    Input Parameters:
686: +  ksp   - iterative context
687: .  n     - iteration number
688: .  rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
689: -  ctx - convergence context which must be created by KSPConvergedDefaultCreate()

691:    Output Parameter:
692: +   positive - if the iteration has converged;
693: .   negative - if residual norm exceeds divergence threshold;
694: -   0 - otherwise.

696:    Notes:
697:    KSPConvergedDefault() reaches convergence when   rnorm < MAX (rtol * rnorm_0, abstol);
698:    Divergence is detected if  rnorm > dtol * rnorm_0,

700:    where:
701: +     rtol = relative tolerance,
702: .     abstol = absolute tolerance.
703: .     dtol = divergence tolerance,
704: -     rnorm_0 is the two norm of the right hand side. When initial guess is non-zero you
705:           can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
706:           as the starting point for relative norm convergence testing, that is as rnorm_0

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

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

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

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

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

718:    Level: intermediate

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

722: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
723:           KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
724: @*/
725: PetscErrorCode  KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
726: {
727:   PetscErrorCode         ierr;
728:   KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
729:   KSPNormType            normtype;

734:   *reason = KSP_CONVERGED_ITERATING;

736:   KSPGetNormType(ksp,&normtype);
737:   if (normtype == KSP_NORM_NONE) return(0);

739:   if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
740:   if (!n) {
741:     /* if user gives initial guess need to compute norm of b */
742:     if (!ksp->guess_zero && !cctx->initialrtol) {
743:       PetscReal snorm = 0.0;
744:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
745:         PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
746:         VecNorm(ksp->vec_rhs,NORM_2,&snorm);        /*     <- b'*b */
747:       } else {
748:         Vec z;
749:         /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
750:         VecDuplicate(ksp->vec_rhs,&z);
751:         KSP_PCApply(ksp,ksp->vec_rhs,z);
752:         if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
753:           PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
754:           VecNorm(z,NORM_2,&snorm);                 /*    dp <- b'*B'*B*b */
755:         } else if (ksp->normtype == KSP_NORM_NATURAL) {
756:           PetscScalar norm;
757:           PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
758:           VecDot(ksp->vec_rhs,z,&norm);
759:           snorm = PetscSqrtReal(PetscAbsScalar(norm));                            /*    dp <- b'*B*b */
760:         }
761:         VecDestroy(&z);
762:       }
763:       /* handle special case of zero RHS and nonzero guess */
764:       if (!snorm) {
765:         PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
766:         snorm = rnorm;
767:       }
768:       if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
769:       else ksp->rnorm0 = snorm;
770:     } else {
771:       ksp->rnorm0 = rnorm;
772:     }
773:     ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
774:   }

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

778:   if (PetscIsInfOrNanReal(rnorm)) {
779:     PCFailedReason pcreason;
780:     PetscInt       sendbuf,pcreason_max;
781:     PCGetSetUpFailedReason(ksp->pc,&pcreason);
782:     sendbuf = (PetscInt)pcreason;
783:     MPI_Allreduce(&sendbuf,&pcreason_max,1,MPIU_INT,MPIU_MAX,PetscObjectComm((PetscObject)ksp));
784:     if (pcreason_max) {
785:       *reason = KSP_DIVERGED_PCSETUP_FAILED;
786:       VecSetInf(ksp->vec_sol);
787:       PetscInfo(ksp,"Linear solver pcsetup fails, declaring divergence \n");
788:     } else {
789:       *reason = KSP_DIVERGED_NANORINF;
790:       PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
791:     }
792:   } else if (rnorm <= ksp->ttol) {
793:     if (rnorm < ksp->abstol) {
794:       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);
795:       *reason = KSP_CONVERGED_ATOL;
796:     } else {
797:       if (cctx->initialrtol) {
798:         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);
799:       } else {
800:         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);
801:       }
802:       *reason = KSP_CONVERGED_RTOL;
803:     }
804:   } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
805:     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);
806:     *reason = KSP_DIVERGED_DTOL;
807:   }
808:   return(0);
809: }

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

816:    Collective on KSP

818:    Input Parameters:
819: .  ctx - convergence context

821:    Level: intermediate

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

825: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
826:           KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
827: @*/
828: PetscErrorCode  KSPConvergedDefaultDestroy(void *ctx)
829: {
830:   PetscErrorCode         ierr;
831:   KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;

834:   VecDestroy(&cctx->work);
835:   PetscFree(ctx);
836:   return(0);
837: }

841: /*
842:    KSPBuildSolutionDefault - Default code to create/move the solution.

844:    Input Parameters:
845: +  ksp - iterative context
846: -  v   - pointer to the user's vector

848:    Output Parameter:
849: .  V - pointer to a vector containing the solution

851:    Level: advanced

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

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

857: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
858: */
859: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
860: {

864:   if (ksp->pc_side == PC_RIGHT) {
865:     if (ksp->pc) {
866:       if (v) {
867:         KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
868:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
869:     } else {
870:       if (v) {
871:         VecCopy(ksp->vec_sol,v); *V = v;
872:       } else *V = ksp->vec_sol;
873:     }
874:   } else if (ksp->pc_side == PC_SYMMETRIC) {
875:     if (ksp->pc) {
876:       if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
877:       if (v) {
878:         PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
879:         *V = v;
880:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
881:     } else {
882:       if (v) {
883:         VecCopy(ksp->vec_sol,v); *V = v;
884:       } else *V = ksp->vec_sol;
885:     }
886:   } else {
887:     if (v) {
888:       VecCopy(ksp->vec_sol,v); *V = v;
889:     } else *V = ksp->vec_sol;
890:   }
891:   return(0);
892: }

896: /*
897:    KSPBuildResidualDefault - Default code to compute the residual.

899:    Input Parameters:
900: .  ksp - iterative context
901: .  t   - pointer to temporary vector
902: .  v   - pointer to user vector

904:    Output Parameter:
905: .  V - pointer to a vector containing the residual

907:    Level: advanced

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

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

913: .seealso: KSPBuildSolutionDefault()
914: */
915: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
916: {
918:   Mat            Amat,Pmat;

921:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
922:   PCGetOperators(ksp->pc,&Amat,&Pmat);
923:   KSPBuildSolution(ksp,t,NULL);
924:   KSP_MatMult(ksp,Amat,t,v);
925:   VecAYPX(v,-1.0,ksp->vec_rhs);
926:   *V   = v;
927:   return(0);
928: }

932: /*@C
933:   KSPCreateVecs - Gets a number of work vectors.

935:   Input Parameters:
936: + ksp  - iterative context
937: . rightn  - number of right work vectors
938: - leftn   - number of left work vectors to allocate

940:   Output Parameter:
941: +  right - the array of vectors created
942: -  left - the array of left vectors

944:    Note: The right vector has as many elements as the matrix has columns. The left
945:      vector has as many elements as the matrix has rows.

947:    The vectors are new vectors that are not owned by the KSP, they should be destroyed with calls to VecDestroyVecs() when no longer needed.

949:    Developers Note: First tries to duplicate the rhs and solution vectors of the KSP, if they do not exist tries to get them from the matrix, if
950:                     that does not exist tries to get them from the DM (if it is provided).

952:    Level: advanced

954: .seealso:   MatCreateVecs(), VecDestroyVecs()

956: @*/
957: PetscErrorCode KSPCreateVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
958: {
960:   Vec            vecr = NULL,vecl = NULL;
961:   PetscBool      matset,pmatset;
962:   Mat            mat = NULL;

965:   if (rightn) {
966:     if (!right) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for right vectors but did not pass a pointer to hold them");
967:     if (ksp->vec_sol) vecr = ksp->vec_sol;
968:     else {
969:       if (ksp->pc) {
970:         PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
971:         /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
972:         if (matset) {
973:           PCGetOperators(ksp->pc,&mat,NULL);
974:           MatCreateVecs(mat,&vecr,NULL);
975:         } else if (pmatset) {
976:           PCGetOperators(ksp->pc,NULL,&mat);
977:           MatCreateVecs(mat,&vecr,NULL);
978:         }
979:       }
980:       if (!vecr) {
981:         if (ksp->dm) {
982:           DMGetGlobalVector(ksp->dm,&vecr);
983:         } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
984:       }
985:     }
986:     VecDuplicateVecs(vecr,rightn,right);
987:     if (!ksp->vec_sol) {
988:       if (mat) {
989:         VecDestroy(&vecr);
990:       } else if (ksp->dm) {
991:         DMRestoreGlobalVector(ksp->dm,&vecr);
992:       }
993:     }
994:   }
995:   if (leftn) {
996:     if (!left) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_INCOMP,"You asked for left vectors but did not pass a pointer to hold them");
997:     if (ksp->vec_rhs) vecl = ksp->vec_rhs;
998:     else {
999:       if (ksp->pc) {
1000:         PCGetOperatorsSet(ksp->pc,&matset,&pmatset);
1001:         /* check for mat before pmat because for KSPLSQR pmat may be a different size than mat since pmat maybe mat'*mat */
1002:         if (matset) {
1003:           PCGetOperators(ksp->pc,&mat,NULL);
1004:           MatCreateVecs(mat,NULL,&vecl);
1005:         } else if (pmatset) {
1006:           PCGetOperators(ksp->pc,NULL,&mat);
1007:           MatCreateVecs(mat,NULL,&vecl);
1008:         }
1009:       }
1010:       if (!vecl) {
1011:         if (ksp->dm) {
1012:           DMGetGlobalVector(ksp->dm,&vecl);
1013:         } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You requested a vector from a KSP that cannot provide one");
1014:       }
1015:     }
1016:     VecDuplicateVecs(vecl,leftn,left);
1017:     if (!ksp->vec_rhs) {
1018:       if (mat) {
1019:         VecDestroy(&vecl);
1020:       } else if (ksp->dm) {
1021:         DMRestoreGlobalVector(ksp->dm,&vecl);
1022:       }
1023:     }
1024:   }
1025:   return(0);
1026: }

1030: /*
1031:   KSPSetWorkVecs - Sets a number of work vectors into a KSP object

1033:   Input Parameters:
1034: . ksp  - iterative context
1035: . nw   - number of work vectors to allocate

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

1039:  */
1040: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1041: {

1045:   VecDestroyVecs(ksp->nwork,&ksp->work);
1046:   ksp->nwork = nw;
1047:   KSPCreateVecs(ksp,nw,&ksp->work,0,NULL);
1048:   PetscLogObjectParents(ksp,nw,ksp->work);
1049:   return(0);
1050: }

1054: /*
1055:   KSPDestroyDefault - Destroys a iterative context variable for methods with
1056:   no separate context.  Preferred calling sequence KSPDestroy().

1058:   Input Parameter:
1059: . ksp - the iterative context

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

1063: */
1064: PetscErrorCode KSPDestroyDefault(KSP ksp)
1065: {

1070:   PetscFree(ksp->data);
1071:   return(0);
1072: }

1076: /*@
1077:    KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.

1079:    Not Collective

1081:    Input Parameter:
1082: .  ksp - the KSP context

1084:    Output Parameter:
1085: .  reason - negative value indicates diverged, positive value converged, see KSPConvergedReason

1087:    Possible values for reason: See also manual page for each reason
1088: $  KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
1089: $  KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
1090: $  KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration, or when the KSPConvergedSkip() convergence test routine is set.
1091: $  KSP_CONVERGED_CG_NEG_CURVE
1092: $  KSP_CONVERGED_CG_CONSTRAINED
1093: $  KSP_CONVERGED_STEP_LENGTH
1094: $  KSP_DIVERGED_ITS  (required more than its to reach convergence)
1095: $  KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
1096: $  KSP_DIVERGED_NANORINF (residual norm became Not-a-number or Inf likely due to 0/0)
1097: $  KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
1098: $  KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial residual. Try a different preconditioner, or a different initial Level.)

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

1102:    Level: intermediate

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

1106: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason
1107: @*/
1108: PetscErrorCode  KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1109: {
1113:   *reason = ksp->reason;
1114:   return(0);
1115: }

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

1123:    Logically Collective on KSP

1125:    Input Parameters:
1126: +  ksp - the preconditioner context
1127: -  dm - the dm

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

1133:    Level: intermediate

1135: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1136: @*/
1137: PetscErrorCode  KSPSetDM(KSP ksp,DM dm)
1138: {
1140:   PC             pc;

1144:   if (dm) {PetscObjectReference((PetscObject)dm);}
1145:   if (ksp->dm) {                /* Move the DMSNES context over to the new DM unless the new DM already has one */
1146:     if (ksp->dm->dmksp && ksp->dmAuto && !dm->dmksp) {
1147:       DMKSP kdm;
1148:       DMCopyDMKSP(ksp->dm,dm);
1149:       DMGetDMKSP(ksp->dm,&kdm);
1150:       if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1151:     }
1152:     DMDestroy(&ksp->dm);
1153:   }
1154:   ksp->dm       = dm;
1155:   ksp->dmAuto   = PETSC_FALSE;
1156:   KSPGetPC(ksp,&pc);
1157:   PCSetDM(pc,dm);
1158:   ksp->dmActive = PETSC_TRUE;
1159:   return(0);
1160: }

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

1167:    Logically Collective on KSP

1169:    Input Parameters:
1170: +  ksp - the preconditioner context
1171: -  flg - use the DM

1173:    Notes:
1174:    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.

1176:    Level: intermediate

1178: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1179: @*/
1180: PetscErrorCode  KSPSetDMActive(KSP ksp,PetscBool flg)
1181: {
1185:   ksp->dmActive = flg;
1186:   return(0);
1187: }

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

1194:    Not Collective

1196:    Input Parameter:
1197: . ksp - the preconditioner context

1199:    Output Parameter:
1200: .  dm - the dm

1202:    Level: intermediate


1205: .seealso: KSPSetDM(), KSPSetDMActive()
1206: @*/
1207: PetscErrorCode  KSPGetDM(KSP ksp,DM *dm)
1208: {

1213:   if (!ksp->dm) {
1214:     DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1215:     ksp->dmAuto = PETSC_TRUE;
1216:   }
1217:   *dm = ksp->dm;
1218:   return(0);
1219: }

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

1226:    Logically Collective on KSP

1228:    Input Parameters:
1229: +  ksp - the KSP context
1230: -  usrP - optional user context

1232:    Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
1233:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

1235:    Level: intermediate

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

1239: .seealso: KSPGetApplicationContext()
1240: @*/
1241: PetscErrorCode  KSPSetApplicationContext(KSP ksp,void *usrP)
1242: {
1244:   PC             pc;

1248:   ksp->user = usrP;
1249:   KSPGetPC(ksp,&pc);
1250:   PCSetApplicationContext(pc,usrP);
1251:   return(0);
1252: }

1256: /*@
1257:    KSPGetApplicationContext - Gets the user-defined context for the linear solver.

1259:    Not Collective

1261:    Input Parameter:
1262: .  ksp - KSP context

1264:    Output Parameter:
1265: .  usrP - user context

1267:    Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
1268:     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.

1270:    Level: intermediate

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

1274: .seealso: KSPSetApplicationContext()
1275: @*/
1276: PetscErrorCode  KSPGetApplicationContext(KSP ksp,void *usrP)
1277: {
1280:   *(void**)usrP = ksp->user;
1281:   return(0);
1282: }