Actual source code: iterativ.c

petsc-3.6.1 2015-07-22
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,void *dummy)
132: {
133:   PetscReal      emin,emax,c;
135:   PetscViewer    viewer = (PetscViewer) dummy;

139:   if (!viewer) {
140:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
141:   }
142:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
143:   if (!ksp->calc_sings) {
144:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
145:   } else {
146:     KSPComputeExtremeSingularValues(ksp,&emax,&emin);
147:     c    = emax/emin;
148:     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);
149:   }
150:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
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 - either a viewer or NULL

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,void *dummy)
179: {
181:   Vec            x;
182:   PetscViewer    viewer = (PetscViewer) dummy;

185:   KSPBuildSolution(ksp,NULL,&x);
186:   if (!viewer) {
187:     MPI_Comm comm;
188:     PetscObjectGetComm((PetscObject)ksp,&comm);
189:     viewer = PETSC_VIEWER_DRAW_(comm);
190:   }
191:   VecView(x,viewer);
192:   return(0);
193: }

195: #if defined(PETSC_HAVE_THREADSAFETY)
196: #define KSPMonitorDefault KSPMonitorDefaultUnsafe
197: #endif

201: /*@C
202:    KSPMonitorDefault - Print the residual norm at each iteration of an
203:    iterative solver.

205:    Collective on KSP

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

213:    Level: intermediate

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

217: .seealso: KSPMonitorSet(), KSPMonitorTrueResidualNorm(), KSPMonitorLGResidualNormCreate()
218: @*/
219: PetscErrorCode  KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
220: {
222:   PetscViewer    viewer = (PetscViewer) dummy;

225:   if (!viewer) {
226:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
227:   }
228:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
229:   if (n == 0 && ((PetscObject)ksp)->prefix) {
230:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
231:   }
232:   PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,(double)rnorm);
233:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
234:   return(0);
235: }

237: #if defined(PETSC_HAVE_THREADSAFETY)
238: #undef KSPMonitorDefault
239: PetscErrorCode KSPMonitorDefault(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
240: {
242: #pragma omp critical
243:   KSPMonitorDefaultUnsafe(ksp,n,rnorm,dummy);
244:   return ierr;
245: }
246: #endif

250: /*@C
251:    KSPMonitorTrueResidualNorm - Prints the true residual norm as well as the preconditioned
252:    residual norm at each iteration of an iterative solver.

254:    Collective on KSP

256:    Input Parameters:
257: +  ksp   - iterative context
258: .  n     - iteration number
259: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
260: -  dummy - unused monitor context

262:    Options Database Key:
263: .  -ksp_monitor_true_residual - Activates KSPMonitorTrueResidualNorm()

265:    Notes:
266:    When using right preconditioning, these values are equivalent.

268:    Level: intermediate

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

272: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualMaxNorm()
273: @*/
274: PetscErrorCode  KSPMonitorTrueResidualNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
275: {
277:   Vec            resid;
278:   PetscReal      truenorm,bnorm;
279:   PetscViewer    viewer = (PetscViewer)dummy;
280:   char           normtype[256];

283:   if (!viewer) {
284:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
285:   }
286:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
287:   if (n == 0 && ((PetscObject)ksp)->prefix) {
288:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
289:   }
290:   KSPBuildResidual(ksp,NULL,NULL,&resid);
291:   VecNorm(resid,NORM_2,&truenorm);
292:   VecDestroy(&resid);
293:   VecNorm(ksp->vec_rhs,NORM_2,&bnorm);
294:   PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
295:   PetscStrtolower(normtype);
296:   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));
297:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
298:   return(0);
299: }

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

307:    Collective on KSP

309:    Input Parameters:
310: +  ksp   - iterative context
311: .  n     - iteration number
312: .  rnorm - norm (preconditioned) residual value (may be estimated).
313: -  dummy - unused monitor context

315:    Options Database Key:
316: .  -ksp_monitor_max - Activates KSPMonitorTrueResidualMaxNorm()

318:    Notes:
319:    This could be implemented (better) with a flag in ksp.

321:    Level: intermediate

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

325: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(),KSPMonitorTrueResidualNorm()
326: @*/
327: PetscErrorCode  KSPMonitorTrueResidualMaxNorm(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
328: {
330:   Vec            resid;
331:   PetscReal      truenorm,bnorm;
332:   PetscViewer    viewer = (PetscViewer)dummy;
333:   char           normtype[256];

336:   if (!viewer) {
337:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
338:   }
339:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
340:   if (n == 0 && ((PetscObject)ksp)->prefix) {
341:     PetscViewerASCIIPrintf(viewer,"  Residual norms (max) for %s solve.\n",((PetscObject)ksp)->prefix);
342:   }
343:   KSPBuildResidual(ksp,NULL,NULL,&resid);
344:   VecNorm(resid,NORM_INFINITY,&truenorm);
345:   VecDestroy(&resid);
346:   VecNorm(ksp->vec_rhs,NORM_INFINITY,&bnorm);
347:   PetscStrncpy(normtype,KSPNormTypes[ksp->normtype],sizeof(normtype));
348:   PetscStrtolower(normtype);
349:   /* 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)); */
350:   PetscViewerASCIIPrintf(viewer,"%3D KSP true resid max norm %14.12e ||r(i)||/||b|| %14.12e\n",n,(double)truenorm,(double)(truenorm/bnorm));
351:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
352:   return(0);
353: }

357: PetscErrorCode  KSPMonitorRange_Private(KSP ksp,PetscInt it,PetscReal *per)
358: {
360:   Vec            resid;
361:   PetscReal      rmax,pwork;
362:   PetscInt       i,n,N;
363:   const PetscScalar *r;

366:   KSPBuildResidual(ksp,NULL,NULL,&resid);
367:   VecNorm(resid,NORM_INFINITY,&rmax);
368:   VecGetLocalSize(resid,&n);
369:   VecGetSize(resid,&N);
370:   VecGetArrayRead(resid,&r);
371:   pwork = 0.0;
372:   for (i=0; i<n; i++) pwork += (PetscAbsScalar(r[i]) > .20*rmax);
373:   VecRestoreArrayRead(resid,&r);
374:   VecDestroy(&resid);
375:   MPI_Allreduce(&pwork,per,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)ksp));
376:   *per = *per/N;
377:   return(0);
378: }

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

385:    Collective on KSP

387:    Input Parameters:
388: +  ksp   - iterative context
389: .  it    - iteration number
390: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
391: -  dummy - unused monitor context

393:    Options Database Key:
394: .  -ksp_monitor_range - Activates KSPMonitorRange()

396:    Level: intermediate

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

400: .seealso: KSPMonitorSet(), KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
401: @*/
402: PetscErrorCode  KSPMonitorRange(KSP ksp,PetscInt it,PetscReal rnorm,void *dummy)
403: {
404:   PetscErrorCode   ierr;
405:   PetscReal        perc,rel;
406:   PetscViewer      viewer = (PetscViewer)dummy;
407:   /* should be in a MonitorRangeContext */
408:   static PetscReal prev;

411:   if (!viewer) {
412:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
413:   }
414:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
415:   if (!it) prev = rnorm;
416:   if (it == 0 && ((PetscObject)ksp)->prefix) {
417:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
418:   }
419:   KSPMonitorRange_Private(ksp,it,&perc);

421:   rel  = (prev - rnorm)/prev;
422:   prev = rnorm;
423:   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));
424:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
425:   return(0);
426: }

430: /*@C
431:    KSPMonitorDynamicTolerance - Recompute the inner tolerance in every
432:    outer iteration in an adaptive way.

434:    Collective on KSP

436:    Input Parameters:
437: +  ksp   - iterative context
438: .  n     - iteration number (not used)
439: .  fnorm - the current residual norm
440: .  dummy - some context as a C struct. fields:
441:              coef: a scaling coefficient. default 1.0. can be passed through
442:                    -sub_ksp_dynamic_tolerance_param
443:              bnrm: norm of the right-hand side. store it to avoid repeated calculation

445:    Notes:
446:    This may be useful for a flexibly preconditioner Krylov method to
447:    control the accuracy of the inner solves needed to gaurantee the
448:    convergence of the outer iterations.

450:    Level: advanced

452: .keywords: KSP, inner tolerance

454: .seealso: KSPMonitorDynamicToleranceDestroy()
455: @*/
456: PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
457: {
459:   PC             pc;
460:   PetscReal      outer_rtol, outer_abstol, outer_dtol, inner_rtol;
461:   PetscInt       outer_maxits,nksp,first,i;
462:   KSPDynTolCtx   *scale   = (KSPDynTolCtx*)dummy;
463:   KSP            kspinner = NULL, *subksp = NULL;

466:   KSPGetPC(ksp, &pc);

468:   /* compute inner_rtol */
469:   if (scale->bnrm < 0.0) {
470:     Vec b;
471:     KSPGetRhs(ksp, &b);
472:     VecNorm(b, NORM_2, &(scale->bnrm));
473:   }
474:   KSPGetTolerances(ksp, &outer_rtol, &outer_abstol, &outer_dtol, &outer_maxits);
475:   inner_rtol = PetscMin(scale->coef * scale->bnrm * outer_rtol / fnorm, 0.999);
476:   /*PetscPrintf(PETSC_COMM_WORLD, "        Inner rtol = %g\n", (double)inner_rtol);*/

478:   /* if pc is ksp */
479:   PCKSPGetKSP(pc, &kspinner);
480:   if (kspinner) {
481:     KSPSetTolerances(kspinner, inner_rtol, outer_abstol, outer_dtol, outer_maxits);
482:     return(0);
483:   }

485:   /* if pc is bjacobi */
486:   PCBJacobiGetSubKSP(pc, &nksp, &first, &subksp);
487:   if (subksp) {
488:     for (i=0; i<nksp; i++) {
489:       KSPSetTolerances(subksp[i], inner_rtol, outer_abstol, outer_dtol, outer_maxits);
490:     }
491:     return(0);
492:   }

494:   /* todo: dynamic tolerance may apply to other types of pc too */
495:   return(0);
496: }

500: /*
501:   Destroy the dummy context used in KSPMonitorDynamicTolerance()
502: */
503: PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy)
504: {

508:   PetscFree(*dummy);
509:   return(0);
510: }

514: /*
515:   Default (short) KSP Monitor, same as KSPMonitorDefault() except
516:   it prints fewer digits of the residual as the residual gets smaller.
517:   This is because the later digits are meaningless and are often
518:   different on different machines; by using this routine different
519:   machines will usually generate the same output.
520: */
521: PetscErrorCode  KSPMonitorDefaultShort(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
522: {
524:   PetscViewer    viewer = (PetscViewer)dummy;

527:   if (!viewer) {
528:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)ksp),&viewer);
529:   }
530:   PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
531:   if (its == 0 && ((PetscObject)ksp)->prefix) {
532:     PetscViewerASCIIPrintf(viewer,"  Residual norms for %s solve.\n",((PetscObject)ksp)->prefix);
533:   }

535:   if (fnorm > 1.e-9) {
536:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %g \n",its,(double)fnorm);
537:   } else if (fnorm > 1.e-11) {
538:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,(double)fnorm);
539:   } else {
540:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);
541:   }
542:   PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
543:   return(0);
544: }

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

552:    Collective on KSP

554:    Input Parameters:
555: +  ksp   - iterative context
556: .  n     - iteration number
557: .  rnorm - 2-norm residual value (may be estimated)
558: -  dummy - unused convergence context

560:    Returns:
561: .  reason - KSP_CONVERGED_ITERATING, KSP_CONVERGED_ITS

563:    Notes:
564:    This should be used as the convergence test with the option
565:    KSPSetNormType(ksp,KSP_NORM_NONE), since norms of the residual are
566:    not computed. Convergence is then declared after the maximum number
567:    of iterations have been reached. Useful when one is using CG or
568:    BiCGStab as a smoother.

570:    Level: advanced

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

574: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
575: @*/
576: PetscErrorCode  KSPConvergedSkip(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
577: {
581:   *reason = KSP_CONVERGED_ITERATING;
582:   if (n >= ksp->max_it) *reason = KSP_CONVERGED_ITS;
583:   return(0);
584: }


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

592:    Collective on KSP

594:    Output Parameter:
595: .  ctx - convergence context

597:    Level: intermediate

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

601: .seealso: KSPConvergedDefault(), KSPConvergedDefaultDestroy(), KSPSetConvergenceTest(), KSPSetTolerances(),
602:           KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
603: @*/
604: PetscErrorCode  KSPConvergedDefaultCreate(void **ctx)
605: {
606:   PetscErrorCode         ierr;
607:   KSPConvergedDefaultCtx *cctx;

610:   PetscNew(&cctx);
611:   *ctx = cctx;
612:   return(0);
613: }

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

622:    Collective on KSP

624:    Input Parameters:
625: .  ksp   - iterative context

627:    Options Database:
628: .   -ksp_converged_use_initial_residual_norm

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

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

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

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


641:    Level: intermediate

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

645: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUMIRNorm()
646: @*/
647: PetscErrorCode  KSPConvergedDefaultSetUIRNorm(KSP ksp)
648: {
649:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;

653:   if (ksp->converged != KSPConvergedDefault) return(0);
654:   if (ctx->mininitialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
655:   ctx->initialrtol = PETSC_TRUE;
656:   return(0);
657: }

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

666:    Collective on KSP

668:    Input Parameters:
669: .  ksp   - iterative context

671:    Options Database:
672: .   -ksp_converged_use_min_initial_residual_norm

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

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

679:    Level: intermediate

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

683: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm()
684: @*/
685: PetscErrorCode  KSPConvergedDefaultSetUMIRNorm(KSP ksp)
686: {
687:   KSPConvergedDefaultCtx *ctx = (KSPConvergedDefaultCtx*) ksp->cnvP;

691:   if (ksp->converged != KSPConvergedDefault) return(0);
692:   if (ctx->initialrtol) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"Cannot use KSPConvergedDefaultSetUIRNorm() and KSPConvergedDefaultSetUMIRNorm() together");
693:   ctx->mininitialrtol = PETSC_TRUE;
694:   return(0);
695: }

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

702:    Collective on KSP

704:    Input Parameters:
705: +  ksp   - iterative context
706: .  n     - iteration number
707: .  rnorm - residual norm (may be estimated, depending on the method may be the preconditioned residual norm)
708: -  ctx - convergence context which must be created by KSPConvergedDefaultCreate()

710:    Output Parameter:
711: +   positive - if the iteration has converged;
712: .   negative - if residual norm exceeds divergence threshold;
713: -   0 - otherwise.

715:    Notes:
716:    KSPConvergedDefault() reaches convergence when   rnorm < MAX (rtol * rnorm_0, abstol);
717:    Divergence is detected if  rnorm > dtol * rnorm_0,

719:    where:
720: +     rtol = relative tolerance,
721: .     abstol = absolute tolerance.
722: .     dtol = divergence tolerance,
723: -     rnorm_0 is the two norm of the right hand side. When initial guess is non-zero you
724:           can call KSPConvergedDefaultSetUIRNorm() to use the norm of (b - A*(initial guess))
725:           as the starting point for relative norm convergence testing, that is as rnorm_0

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

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

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

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

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

737:    Level: intermediate

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

741: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(), KSPConvergedReason, KSPGetConvergedReason(),
742:           KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm(), KSPConvergedDefaultCreate(), KSPConvergedDefaultDestroy()
743: @*/
744: PetscErrorCode  KSPConvergedDefault(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx)
745: {
746:   PetscErrorCode         ierr;
747:   KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;
748:   KSPNormType            normtype;

753:   *reason = KSP_CONVERGED_ITERATING;

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

758:   if (!cctx) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_NULL,"Convergence context must have been created with KSPConvergedDefaultCreate()");
759:   if (!n) {
760:     /* if user gives initial guess need to compute norm of b */
761:     if (!ksp->guess_zero && !cctx->initialrtol) {
762:       PetscReal snorm;
763:       if (ksp->normtype == KSP_NORM_UNPRECONDITIONED || ksp->pc_side == PC_RIGHT) {
764:         PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of RHS\n");
765:         VecNorm(ksp->vec_rhs,NORM_2,&snorm);        /*     <- b'*b */
766:       } else {
767:         Vec z;
768:         /* Should avoid allocating the z vector each time but cannot stash it in cctx because if KSPReset() is called the vector size might change */
769:         VecDuplicate(ksp->vec_rhs,&z);
770:         KSP_PCApply(ksp,ksp->vec_rhs,z);
771:         if (ksp->normtype == KSP_NORM_PRECONDITIONED) {
772:           PetscInfo(ksp,"user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n");
773:           VecNorm(z,NORM_2,&snorm);                 /*    dp <- b'*B'*B*b */
774:         } else if (ksp->normtype == KSP_NORM_NATURAL) {
775:           PetscScalar norm;
776:           PetscInfo(ksp,"user has provided nonzero initial guess, computing natural norm of RHS\n");
777:           VecDot(ksp->vec_rhs,z,&norm);
778:           snorm = PetscSqrtReal(PetscAbsScalar(norm));                            /*    dp <- b'*B*b */
779:         }
780:         VecDestroy(&z);
781:       }
782:       /* handle special case of zero RHS and nonzero guess */
783:       if (!snorm) {
784:         PetscInfo(ksp,"Special case, user has provided nonzero initial guess and zero RHS\n");
785:         snorm = rnorm;
786:       }
787:       if (cctx->mininitialrtol) ksp->rnorm0 = PetscMin(snorm,rnorm);
788:       else ksp->rnorm0 = snorm;
789:     } else {
790:       ksp->rnorm0 = rnorm;
791:     }
792:     ksp->ttol = PetscMax(ksp->rtol*ksp->rnorm0,ksp->abstol);
793:   }

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

797:   if (PetscIsInfOrNanReal(rnorm)) {
798:     PetscInfo(ksp,"Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n");
799:     *reason = KSP_DIVERGED_NANORINF;
800:   } else if (rnorm <= ksp->ttol) {
801:     if (rnorm < ksp->abstol) {
802:       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);
803:       *reason = KSP_CONVERGED_ATOL;
804:     } else {
805:       if (cctx->initialrtol) {
806:         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);
807:       } else {
808:         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);
809:       }
810:       *reason = KSP_CONVERGED_RTOL;
811:     }
812:   } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
813:     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);
814:     *reason = KSP_DIVERGED_DTOL;
815:   }
816:   return(0);
817: }

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

824:    Collective on KSP

826:    Input Parameters:
827: .  ctx - convergence context

829:    Level: intermediate

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

833: .seealso: KSPConvergedDefault(), KSPConvergedDefaultCreate(), KSPSetConvergenceTest(), KSPSetTolerances(), KSPConvergedSkip(),
834:           KSPConvergedReason, KSPGetConvergedReason(), KSPConvergedDefaultSetUIRNorm(), KSPConvergedDefaultSetUMIRNorm()
835: @*/
836: PetscErrorCode  KSPConvergedDefaultDestroy(void *ctx)
837: {
838:   PetscErrorCode         ierr;
839:   KSPConvergedDefaultCtx *cctx = (KSPConvergedDefaultCtx*) ctx;

842:   VecDestroy(&cctx->work);
843:   PetscFree(ctx);
844:   return(0);
845: }

849: /*
850:    KSPBuildSolutionDefault - Default code to create/move the solution.

852:    Input Parameters:
853: +  ksp - iterative context
854: -  v   - pointer to the user's vector

856:    Output Parameter:
857: .  V - pointer to a vector containing the solution

859:    Level: advanced

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

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

865: .seealso: KSPGetSolution(), KSPBuildResidualDefault()
866: */
867: PetscErrorCode KSPBuildSolutionDefault(KSP ksp,Vec v,Vec *V)
868: {

872:   if (ksp->pc_side == PC_RIGHT) {
873:     if (ksp->pc) {
874:       if (v) {
875:         KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;
876:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with right preconditioner");
877:     } else {
878:       if (v) {
879:         VecCopy(ksp->vec_sol,v); *V = v;
880:       } else *V = ksp->vec_sol;
881:     }
882:   } else if (ksp->pc_side == PC_SYMMETRIC) {
883:     if (ksp->pc) {
884:       if (ksp->transpose_solve) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
885:       if (v) {
886:         PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v);
887:         *V = v;
888:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Not working with symmetric preconditioner");
889:     } else {
890:       if (v) {
891:         VecCopy(ksp->vec_sol,v); *V = v;
892:       } else *V = ksp->vec_sol;
893:     }
894:   } else {
895:     if (v) {
896:       VecCopy(ksp->vec_sol,v); *V = v;
897:     } else *V = ksp->vec_sol;
898:   }
899:   return(0);
900: }

904: /*
905:    KSPBuildResidualDefault - Default code to compute the residual.

907:    Input Parameters:
908: .  ksp - iterative context
909: .  t   - pointer to temporary vector
910: .  v   - pointer to user vector

912:    Output Parameter:
913: .  V - pointer to a vector containing the residual

915:    Level: advanced

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

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

921: .seealso: KSPBuildSolutionDefault()
922: */
923: PetscErrorCode KSPBuildResidualDefault(KSP ksp,Vec t,Vec v,Vec *V)
924: {
926:   Mat            Amat,Pmat;

929:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
930:   PCGetOperators(ksp->pc,&Amat,&Pmat);
931:   KSPBuildSolution(ksp,t,NULL);
932:   KSP_MatMult(ksp,Amat,t,v);
933:   VecAYPX(v,-1.0,ksp->vec_rhs);
934:   *V   = v;
935:   return(0);
936: }

940: /*@C
941:   KSPCreateVecs - Gets a number of work vectors.

943:   Input Parameters:
944: + ksp  - iterative context
945: . rightn  - number of right work vectors
946: - leftn   - number of left work vectors to allocate

948:   Output Parameter:
949: +  right - the array of vectors created
950: -  left - the array of left vectors

952:    Note: The right vector has as many elements as the matrix has columns. The left
953:      vector has as many elements as the matrix has rows.

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

957:    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
958:                     that does not exist tries to get them from the DM (if it is provided).

960:    Level: advanced

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

964: @*/
965: PetscErrorCode KSPCreateVecs(KSP ksp,PetscInt rightn, Vec **right,PetscInt leftn,Vec **left)
966: {
968:   Vec            vecr = NULL,vecl = NULL;
969:   PetscBool      matset,pmatset;
970:   Mat            mat = NULL;

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

1036: /*
1037:   KSPSetWorkVecs - Sets a number of work vectors into a KSP object

1039:   Input Parameters:
1040: . ksp  - iterative context
1041: . nw   - number of work vectors to allocate

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

1045:  */
1046: PetscErrorCode KSPSetWorkVecs(KSP ksp,PetscInt nw)
1047: {

1051:   VecDestroyVecs(ksp->nwork,&ksp->work);
1052:   ksp->nwork = nw;
1053:   KSPCreateVecs(ksp,nw,&ksp->work,0,NULL);
1054:   PetscLogObjectParents(ksp,nw,ksp->work);
1055:   return(0);
1056: }

1060: /*
1061:   KSPDestroyDefault - Destroys a iterative context variable for methods with
1062:   no separate context.  Preferred calling sequence KSPDestroy().

1064:   Input Parameter:
1065: . ksp - the iterative context

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

1069: */
1070: PetscErrorCode KSPDestroyDefault(KSP ksp)
1071: {

1076:   PetscFree(ksp->data);
1077:   return(0);
1078: }

1082: /*@
1083:    KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.

1085:    Not Collective

1087:    Input Parameter:
1088: .  ksp - the KSP context

1090:    Output Parameter:
1091: .  reason - negative value indicates diverged, positive value converged, see KSPConvergedReason

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

1108:    See also manual page for each reason.

1110:    guess: beginner

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

1114:    Level: intermediate

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

1118: .seealso: KSPSetConvergenceTest(), KSPConvergedDefault(), KSPSetTolerances(), KSPConvergedReason
1119: @*/
1120: PetscErrorCode  KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
1121: {
1125:   *reason = ksp->reason;
1126:   return(0);
1127: }

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

1135:    Logically Collective on KSP

1137:    Input Parameters:
1138: +  ksp - the preconditioner context
1139: -  dm - the dm

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

1145:    Level: intermediate

1147: .seealso: KSPGetDM(), KSPSetDMActive(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess(), DMKSPSetComputeOperators(), DMKSPSetComputeRHS(), DMKSPSetComputeInitialGuess()
1148: @*/
1149: PetscErrorCode  KSPSetDM(KSP ksp,DM dm)
1150: {
1152:   PC             pc;

1156:   if (dm) {PetscObjectReference((PetscObject)dm);}
1157:   if (ksp->dm) {                /* Move the DMSNES context over to the new DM unless the new DM already has one */
1158:     if (ksp->dm->dmksp && ksp->dmAuto && !dm->dmksp) {
1159:       DMKSP kdm;
1160:       DMCopyDMKSP(ksp->dm,dm);
1161:       DMGetDMKSP(ksp->dm,&kdm);
1162:       if (kdm->originaldm == ksp->dm) kdm->originaldm = dm; /* Grant write privileges to the replacement DM */
1163:     }
1164:     DMDestroy(&ksp->dm);
1165:   }
1166:   ksp->dm       = dm;
1167:   ksp->dmAuto   = PETSC_FALSE;
1168:   KSPGetPC(ksp,&pc);
1169:   PCSetDM(pc,dm);
1170:   ksp->dmActive = PETSC_TRUE;
1171:   return(0);
1172: }

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

1179:    Logically Collective on KSP

1181:    Input Parameters:
1182: +  ksp - the preconditioner context
1183: -  flg - use the DM

1185:    Notes:
1186:    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.

1188:    Level: intermediate

1190: .seealso: KSPGetDM(), KSPSetDM(), SNESSetDM(), KSPSetComputeOperators(), KSPSetComputeRHS(), KSPSetComputeInitialGuess()
1191: @*/
1192: PetscErrorCode  KSPSetDMActive(KSP ksp,PetscBool flg)
1193: {
1197:   ksp->dmActive = flg;
1198:   return(0);
1199: }

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

1206:    Not Collective

1208:    Input Parameter:
1209: . ksp - the preconditioner context

1211:    Output Parameter:
1212: .  dm - the dm

1214:    Level: intermediate


1217: .seealso: KSPSetDM(), KSPSetDMActive()
1218: @*/
1219: PetscErrorCode  KSPGetDM(KSP ksp,DM *dm)
1220: {

1225:   if (!ksp->dm) {
1226:     DMShellCreate(PetscObjectComm((PetscObject)ksp),&ksp->dm);
1227:     ksp->dmAuto = PETSC_TRUE;
1228:   }
1229:   *dm = ksp->dm;
1230:   return(0);
1231: }

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

1238:    Logically Collective on KSP

1240:    Input Parameters:
1241: +  ksp - the KSP context
1242: -  usrP - optional user context

1244:    Level: intermediate

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

1248: .seealso: KSPGetApplicationContext()
1249: @*/
1250: PetscErrorCode  KSPSetApplicationContext(KSP ksp,void *usrP)
1251: {
1253:   PC             pc;

1257:   ksp->user = usrP;
1258:   KSPGetPC(ksp,&pc);
1259:   PCSetApplicationContext(pc,usrP);
1260:   return(0);
1261: }

1265: /*@
1266:    KSPGetApplicationContext - Gets the user-defined context for the linear solver.

1268:    Not Collective

1270:    Input Parameter:
1271: .  ksp - KSP context

1273:    Output Parameter:
1274: .  usrP - user context

1276:    Level: intermediate

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

1280: .seealso: KSPSetApplicationContext()
1281: @*/
1282: PetscErrorCode  KSPGetApplicationContext(KSP ksp,void *usrP)
1283: {
1286:   *(void**)usrP = ksp->user;
1287:   return(0);
1288: }