Actual source code: iterativ.c

      1: #define PETSCKSP_DLL

      3: /*
      4:    This file contains some simple default routines.
      5:    These routines should be SHORT, since they will be included in every
      6:    executable image that uses the iterative routines (note that, through
      7:    the registry system, we provide a way to load only the truely necessary
      8:    files)
      9:  */
     10:  #include src/ksp/ksp/kspimpl.h

     14: /*
     15:   KSPDefaultFreeWork - Free work vectors

     17:   Input Parameters:
     18: . ksp  - iterative context
     19:  */
     20: PetscErrorCode KSPDefaultFreeWork(KSP ksp)
     21: {
     25:   if (ksp->work)  {
     26:     VecDestroyVecs(ksp->work,ksp->nwork);
     27:     ksp->work = PETSC_NULL;
     28:   }
     29:   return(0);
     30: }

     34: /*@C
     35:    KSPGetResidualNorm - Gets the last (approximate preconditioned)
     36:    residual norm that has been computed.
     37:
     38:    Not Collective

     40:    Input Parameters:
     41: .  ksp - the iterative context

     43:    Output Parameters:
     44: .  rnorm - residual norm

     46:    Level: intermediate

     48: .keywords: KSP, get, residual norm

     50: .seealso: KSPBuildResidual()
     51: @*/
     52: PetscErrorCode PETSCKSP_DLLEXPORT KSPGetResidualNorm(KSP ksp,PetscReal *rnorm)
     53: {
     57:   *rnorm = ksp->rnorm;
     58:   return(0);
     59: }

     63: /*@
     64:    KSPGetIterationNumber - Gets the current iteration number; if the
     65:          KSPSolve() is complete, returns the number of iterations
     66:          used.
     67:
     68:    Not Collective

     70:    Input Parameters:
     71: .  ksp - the iterative context

     73:    Output Parameters:
     74: .  its - number of iterations

     76:    Level: intermediate

     78:    Notes:
     79:       During the ith iteration this returns i-1
     80: .keywords: KSP, get, residual norm

     82: .seealso: KSPBuildResidual(), KSPGetResidualNorm()
     83: @*/
     84: PetscErrorCode PETSCKSP_DLLEXPORT KSPGetIterationNumber(KSP ksp,PetscInt *its)
     85: {
     89:   *its = ksp->its;
     90:   return(0);
     91: }

     95: /*@C
     96:     KSPSingularValueMonitor - Prints the two norm of the true residual and
     97:     estimation of the extreme singular values of the preconditioned problem
     98:     at each iteration.
     99:
    100:     Collective on KSP

    102:     Input Parameters:
    103: +   ksp - the iterative context
    104: .   n  - the iteration
    105: -   rnorm - the two norm of the residual

    107:     Options Database Key:
    108: .   -ksp_singmonitor - Activates KSPSingularValueMonitor()

    110:     Notes:
    111:     The CG solver uses the Lanczos technique for eigenvalue computation,
    112:     while GMRES uses the Arnoldi technique; other iterative methods do
    113:     not currently compute singular values.

    115:     Level: intermediate

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

    119: .seealso: KSPComputeExtremeSingularValues()
    120: @*/
    121: PetscErrorCode PETSCKSP_DLLEXPORT KSPSingularValueMonitor(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
    122: {
    123:   PetscReal emin,emax,c;

    128:   if (!ksp->calc_sings) {
    129:     PetscPrintf(ksp->comm,"%3D KSP Residual norm %14.12e \n",n,rnorm);
    130:   } else {
    131:     KSPComputeExtremeSingularValues(ksp,&emax,&emin);
    132:     c = emax/emin;
    133:     PetscPrintf(ksp->comm,"%3D KSP Residual norm %14.12e %% max %g min %g max/min %g\n",n,rnorm,emax,emin,c);
    134:   }
    135:   return(0);
    136: }

    140: /*@C
    141:    KSPVecViewMonitor - Monitors progress of the KSP solvers by calling
    142:    VecView() for the approximate solution at each iteration.

    144:    Collective on KSP

    146:    Input Parameters:
    147: +  ksp - the KSP context
    148: .  its - iteration number
    149: .  fgnorm - 2-norm of residual (or gradient)
    150: -  dummy - either a viewer or PETSC_NULL

    152:    Level: intermediate

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

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

    160: .seealso: KSPSetMonitor(), KSPDefaultMonitor(), VecView()
    161: @*/
    162: PetscErrorCode PETSCKSP_DLLEXPORT KSPVecViewMonitor(KSP ksp,PetscInt its,PetscReal fgnorm,void *dummy)
    163: {
    165:   Vec         x;
    166:   PetscViewer viewer = (PetscViewer) dummy;

    169:   KSPBuildSolution(ksp,PETSC_NULL,&x);
    170:   if (!viewer) {
    171:     MPI_Comm comm;
    172:     PetscObjectGetComm((PetscObject)ksp,&comm);
    173:     viewer = PETSC_VIEWER_DRAW_(comm);
    174:   }
    175:   VecView(x,viewer);

    177:   return(0);
    178: }

    182: /*@C
    183:    KSPDefaultMonitor - Print the residual norm at each iteration of an
    184:    iterative solver.

    186:    Collective on KSP

    188:    Input Parameters:
    189: +  ksp   - iterative context
    190: .  n     - iteration number
    191: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
    192: -  dummy - unused monitor context

    194:    Level: intermediate

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

    198: .seealso: KSPSetMonitor(), KSPTrueMonitor(), KSPLGMonitorCreate()
    199: @*/
    200: PetscErrorCode PETSCKSP_DLLEXPORT KSPDefaultMonitor(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
    201: {
    203:   PetscViewer viewer = (PetscViewer) dummy;

    206:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(ksp->comm);
    207:   PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %14.12e \n",n,rnorm);
    208:   return(0);
    209: }

    213: /*@C
    214:    KSPTrueMonitor - Prints the true residual norm as well as the preconditioned
    215:    residual norm at each iteration of an iterative solver.

    217:    Collective on KSP

    219:    Input Parameters:
    220: +  ksp   - iterative context
    221: .  n     - iteration number
    222: .  rnorm - 2-norm (preconditioned) residual value (may be estimated).
    223: -  dummy - unused monitor context

    225:    Options Database Key:
    226: .  -ksp_truemonitor - Activates KSPTrueMonitor()

    228:    Notes:
    229:    When using right preconditioning, these values are equivalent.

    231:    When using either ICC or ILU preconditioners in BlockSolve95
    232:    (via MATMPIROWBS matrix format), then use this monitor will
    233:    print both the residual norm associated with the original
    234:    (unscaled) matrix.

    236:    Level: intermediate

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

    240: .seealso: KSPSetMonitor(), KSPDefaultMonitor(), KSPLGMonitorCreate()
    241: @*/
    242: PetscErrorCode PETSCKSP_DLLEXPORT KSPTrueMonitor(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
    243: {
    245:   Vec          resid,work;
    246:   PetscReal    scnorm;
    247:   PC           pc;
    248:   Mat          A,B;
    249:   PetscViewer  viewer = (PetscViewer) dummy;
    250:
    252:   VecDuplicate(ksp->vec_rhs,&work);
    253:   KSPBuildResidual(ksp,0,work,&resid);

    255:   /*
    256:      Unscale the residual if the matrix is, for example, a BlockSolve matrix
    257:     but only if both matrices are the same matrix, since only then would
    258:     they be scaled.
    259:   */
    260:   VecCopy(resid,work);
    261:   KSPGetPC(ksp,&pc);
    262:   PCGetOperators(pc,&A,&B,PETSC_NULL);
    263:   if (A == B) {
    264:     MatUnScaleSystem(A,PETSC_NULL,work);
    265:   }
    266:   VecNorm(work,NORM_2,&scnorm);
    267:   VecDestroy(work);
    268:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(ksp->comm);
    269:   PetscViewerASCIIPrintf(viewer,"%3D KSP preconditioned resid norm %14.12e true resid norm %14.12e\n",n,rnorm,scnorm);
    270:   return(0);
    271: }

    275: /*
    276:   Default (short) KSP Monitor, same as KSPDefaultMonitor() except
    277:   it prints fewer digits of the residual as the residual gets smaller.
    278:   This is because the later digits are meaningless and are often
    279:   different on different machines; by using this routine different
    280:   machines will usually generate the same output.
    281: */
    282: PetscErrorCode PETSCKSP_DLLEXPORT KSPDefaultSMonitor(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy)
    283: {
    285:   PetscViewer viewer = (PetscViewer) dummy;

    288:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(ksp->comm);

    290:   if (fnorm > 1.e-9) {
    291:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %g \n",its,fnorm);
    292:   } else if (fnorm > 1.e-11){
    293:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm %5.3e \n",its,fnorm);
    294:   } else {
    295:     PetscViewerASCIIPrintf(viewer,"%3D KSP Residual norm < 1.e-11\n",its);
    296:   }
    297:   return(0);
    298: }

    302: /*@C
    303:    KSPSkipConverged - Convergence test that NEVER returns as converged.

    305:    Collective on KSP

    307:    Input Parameters:
    308: +  ksp   - iterative context
    309: .  n     - iteration number
    310: .  rnorm - 2-norm residual value (may be estimated)
    311: -  dummy - unused convergence context

    313:    Returns:
    314: .  reason - always KSP_CONVERGED_ITERATING

    316:    Notes:
    317:    This is used as the convergence test with the option KSPSetNormType(ksp,KSP_NO_NORM),
    318:    since norms of the residual are not computed. Convergence is then declared
    319:    after a fixed number of iterations have been used. Useful when one is
    320:    using CG or Bi-CG-stab as a smoother.
    321:
    322:    Level: advanced

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

    326: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSetNormType()
    327: @*/
    328: PetscErrorCode PETSCKSP_DLLEXPORT KSPSkipConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
    329: {
    332:   return(0);
    333: }

    337: /*@C
    338:    KSPDefaultConverged - Determines convergence of
    339:    the iterative solvers (default code).

    341:    Collective on KSP

    343:    Input Parameters:
    344: +  ksp   - iterative context
    345: .  n     - iteration number
    346: .  rnorm - 2-norm residual value (may be estimated)
    347: -  dummy - unused convergence context

    349:    Returns:
    350: +   positive - if the iteration has converged;
    351: .   negative - if residual norm exceeds divergence threshold;
    352: -   0 - otherwise.

    354:    Notes:
    355:    KSPDefaultConverged() reaches convergence when
    356: $      rnorm < MAX (rtol * rnorm_0, abstol);
    357:    Divergence is detected if
    358: $      rnorm > dtol * rnorm_0,

    360:    where
    361: +     rtol = relative tolerance,
    362: .     abstol = absolute tolerance.
    363: .     dtol = divergence tolerance,
    364: -     rnorm_0 = initial residual norm

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

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

    371:    Level: intermediate

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

    375: .seealso: KSPSetConvergenceTest(), KSPSetTolerances(), KSPSkipConverged(), KSPConvergedReason, KSPGetConvergedReason()
    376: @*/
    377: PetscErrorCode PETSCKSP_DLLEXPORT KSPDefaultConverged(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *dummy)
    378: {

    384:   *reason = KSP_CONVERGED_ITERATING;

    386:   if (!n) {
    387:     /* if user gives initial guess need to compute norm of b */
    388:     if (!ksp->guess_zero) {
    389:       PetscReal      snorm;
    390:       if (ksp->normtype == KSP_UNPRECONDITIONED_NORM || ksp->pc_side == PC_RIGHT) {
    391:         PetscLogInfo((ksp,"KSPDefaultConverged:user has provided nonzero initial guess, computing 2-norm of RHS\n"));
    392:         VecNorm(ksp->vec_rhs,NORM_2,&snorm);        /*     <- b'*b */
    393:       } else {
    394:         Vec z;
    395:         VecDuplicate(ksp->vec_rhs,&z);
    396:         KSP_PCApply(ksp,ksp->vec_rhs,z);
    397:         if (ksp->normtype == KSP_PRECONDITIONED_NORM) {
    398:           PetscLogInfo((ksp,"KSPDefaultConverged:user has provided nonzero initial guess, computing 2-norm of preconditioned RHS\n"));
    399:           VecNorm(z,NORM_2,&snorm);                 /*    dp <- b'*B'*B*b */
    400:         } else if (ksp->normtype == KSP_NATURAL_NORM) {
    401:           PetscScalar norm;
    402:           PetscLogInfo((ksp,"KSPDefaultConverged:user has provided nonzero initial guess, computing natural norm of RHS\n"));
    403:           VecDot(ksp->vec_rhs,z,&norm);
    404:           snorm = sqrt(PetscAbsScalar(norm));                            /*    dp <- b'*B*b */
    405:         }
    406:         VecDestroy(z);
    407:       }
    408:       /* handle special case of zero RHS and nonzero guess */
    409:       if (!snorm) {
    410:         PetscLogInfo((ksp,"KSPDefaultConverged:Special case, user has provided nonzero initial guess and zero RHS\n"));
    411:         snorm = rnorm;
    412:       }
    413:       ksp->ttol   = PetscMax(ksp->rtol*snorm,ksp->abstol);
    414:       ksp->rnorm0 = snorm;
    415:     } else {
    416:       ksp->ttol   = PetscMax(ksp->rtol*rnorm,ksp->abstol);
    417:       ksp->rnorm0 = rnorm;
    418:     }
    419:   }

    421:   if (rnorm != rnorm) {
    422:     PetscLogInfo((ksp,"KSPDefaultConverged:Linear solver has created a not a number (NaN) as the residual norm, declaring divergence \n"));
    423:     *reason = KSP_DIVERGED_NAN;
    424:   } else if (rnorm <= ksp->ttol) {
    425:     if (rnorm < ksp->abstol) {
    426:       PetscLogInfo((ksp,"KSPDefaultConverged:Linear solver has converged. Residual norm %g is less than absolute tolerance %g at iteration %D\n",rnorm,ksp->abstol,n));
    427:       *reason = KSP_CONVERGED_ATOL;
    428:     } else {
    429:       PetscLogInfo((ksp,"KSPDefaultConverged:Linear solver has converged. Residual norm %g is less than relative tolerance %g times initial right hand side norm %g at iteration %D\n",rnorm,ksp->rtol,ksp->rnorm0,n));
    430:       *reason = KSP_CONVERGED_RTOL;
    431:     }
    432:   } else if (rnorm >= ksp->divtol*ksp->rnorm0) {
    433:     PetscLogInfo((ksp,"KSPDefaultConverged:Linear solver is diverging. Initial right hand size norm %g, current residual norm %g at iteration %D\n",ksp->rnorm0,rnorm,n));
    434:     *reason = KSP_DIVERGED_DTOL;
    435:   }
    436:   return(0);
    437: }

    441: /*
    442:    KSPDefaultBuildSolution - Default code to create/move the solution.

    444:    Input Parameters:
    445: +  ksp - iterative context
    446: -  v   - pointer to the user's vector

    448:    Output Parameter:
    449: .  V - pointer to a vector containing the solution

    451:    Level: advanced

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

    455: .seealso: KSPGetSolution(), KSPDefaultBuildResidual()
    456: */
    457: PetscErrorCode KSPDefaultBuildSolution(KSP ksp,Vec v,Vec *V)
    458: {
    461:   if (ksp->pc_side == PC_RIGHT) {
    462:     if (ksp->pc) {
    463:       if (v) {KSP_PCApply(ksp,ksp->vec_sol,v); *V = v;}
    464:       else {SETERRQ(PETSC_ERR_SUP,"Not working with right preconditioner");}
    465:     } else {
    466:       if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
    467:       else { *V = ksp->vec_sol;}
    468:     }
    469:   } else if (ksp->pc_side == PC_SYMMETRIC) {
    470:     if (ksp->pc) {
    471:       if (ksp->transpose_solve) SETERRQ(PETSC_ERR_SUP,"Not working with symmetric preconditioner and transpose solve");
    472:       if (v) {PCApplySymmetricRight(ksp->pc,ksp->vec_sol,v); *V = v;}
    473:       else {SETERRQ(PETSC_ERR_SUP,"Not working with symmetric preconditioner");}
    474:     } else  {
    475:       if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
    476:       else { *V = ksp->vec_sol;}
    477:     }
    478:   } else {
    479:     if (v) {VecCopy(ksp->vec_sol,v); *V = v;}
    480:     else { *V = ksp->vec_sol; }
    481:   }
    482:   return(0);
    483: }

    487: /*
    488:    KSPDefaultBuildResidual - Default code to compute the residual.

    490:    Input Parameters:
    491: .  ksp - iterative context
    492: .  t   - pointer to temporary vector
    493: .  v   - pointer to user vector

    495:    Output Parameter:
    496: .  V - pointer to a vector containing the residual

    498:    Level: advanced

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

    502: .seealso: KSPDefaultBuildSolution()
    503: */
    504: PetscErrorCode KSPDefaultBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
    505: {
    507:   MatStructure pflag;
    508:   Vec          T;
    509:   PetscScalar  mone = -1.0;
    510:   Mat          Amat,Pmat;

    513:   PCGetOperators(ksp->pc,&Amat,&Pmat,&pflag);
    514:   KSPBuildSolution(ksp,t,&T);
    515:   KSP_MatMult(ksp,Amat,t,v);
    516:   VecAYPX(v,mone,ksp->vec_rhs);
    517:   *V = v;
    518:   return(0);
    519: }

    523: /*
    524:   KSPDefaultGetWork - Gets a number of work vectors.

    526:   Input Parameters:
    527: . ksp  - iterative context
    528: . nw   - number of work vectors to allocate

    530:   Output Parameter:
    531: .  work - the array of vectors created

    533:  */
    534: PetscErrorCode KSPGetVecs(KSP ksp,PetscInt nw,Vec **work)
    535: {
    537:   Vec vec;

    540:   if (ksp->vec_rhs) vec = ksp->vec_rhs;
    541:   else {
    542:     Mat pmat;
    543:     PCGetOperators(ksp->pc,0,&pmat,0);
    544:     MatGetVecs(pmat,&vec,0);
    545:   }
    546:   VecDuplicateVecs(vec,nw,work);
    547:   if (!ksp->vec_rhs) {
    548:     VecDestroy(vec);
    549:   }
    550:   return(0);
    551: }

    555: /*
    556:   KSPDefaultGetWork - Gets a number of work vectors.

    558:   Input Parameters:
    559: . ksp  - iterative context
    560: . nw   - number of work vectors to allocate

    562:   Notes:
    563:   Call this only if no work vectors have been allocated
    564:  */
    565: PetscErrorCode KSPDefaultGetWork(KSP ksp,PetscInt nw)
    566: {

    570:   if (ksp->work) {KSPDefaultFreeWork(ksp);}
    571:   ksp->nwork = nw;
    572:   KSPGetVecs(ksp,nw,&ksp->work);
    573:   PetscLogObjectParents(ksp,nw,ksp->work);
    574:   return(0);
    575: }

    579: /*
    580:   KSPDefaultDestroy - Destroys a iterative context variable for methods with
    581:   no separate context.  Preferred calling sequence KSPDestroy().

    583:   Input Parameter:
    584: . ksp - the iterative context
    585: */
    586: PetscErrorCode KSPDefaultDestroy(KSP ksp)
    587: {

    592:   if (ksp->data) {
    593:     PetscFree(ksp->data);
    594:     ksp->data = 0;
    595:   }

    597:   /* free work vectors */
    598:   KSPDefaultFreeWork(ksp);
    599:   return(0);
    600: }

    604: /*@C
    605:    KSPGetConvergedReason - Gets the reason the KSP iteration was stopped.

    607:    Not Collective

    609:    Input Parameter:
    610: .  ksp - the KSP context

    612:    Output Parameter:
    613: .  reason - negative value indicates diverged, positive value converged, see KSPConvergedReason

    615:    Possible values for reason:
    616: +  KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
    617: .  KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
    618: .  KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration)
    619: .  KSP_CONVERGED_QCG_NEG_CURVE
    620: .  KSP_CONVERGED_QCG_CONSTRAINED
    621: .  KSP_CONVERGED_STEP_LENGTH
    622: .  KSP_DIVERGED_ITS  (required more than its to reach convergence)
    623: .  KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
    624: .  KSP_DIVERGED_NAN (residual norm became Not-a-number likely do to 0/0)
    625: .  KSP_DIVERGED_BREAKDOWN (generic breakdown in method)
    626: -  KSP_DIVERGED_BREAKDOWN_BICG (Initial residual is orthogonal to preconditioned initial
    627:                                 residual. Try a different preconditioner, or a different initial guess.)
    628:

    630:    Level: beginner

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

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

    636: .seealso: KSPSetConvergenceTest(), KSPDefaultConverged(), KSPSetTolerances(), KSPConvergedReason
    637: @*/
    638: PetscErrorCode PETSCKSP_DLLEXPORT KSPGetConvergedReason(KSP ksp,KSPConvergedReason *reason)
    639: {
    643:   *reason = ksp->reason;
    644:   return(0);
    645: }