Actual source code: rii.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */
     10: /*
     11:    SLEPc nonlinear eigensolver: "rii"

     13:    Method: Residual inverse iteration

     15:    Algorithm:

     17:        Simple residual inverse iteration with varying shift.

     19:    References:

     21:        [1] A. Neumaier, "Residual inverse iteration for the nonlinear
     22:            eigenvalue problem", SIAM J. Numer. Anal. 22(5):914-923, 1985.
     23: */

     25: #include <slepc/private/nepimpl.h>
     26: #include <../src/nep/impls/nepdefl.h>

     28: typedef struct {
     29:   PetscInt  max_inner_it;     /* maximum number of Newton iterations */
     30:   PetscInt  lag;              /* interval to rebuild preconditioner */
     31:   PetscBool cctol;            /* constant correction tolerance */
     32:   PetscBool herm;             /* whether the Hermitian version of the scalar equation must be used */
     33:   PetscReal deftol;           /* tolerance for the deflation (threshold) */
     34:   KSP       ksp;              /* linear solver object */
     35: } NEP_RII;

     37: static PetscErrorCode NEPSetUp_RII(NEP nep)
     38: {
     39:   PetscFunctionBegin;
     40:   if (nep->ncv!=PETSC_DETERMINE) PetscCall(PetscInfo(nep,"Setting ncv = nev, ignoring user-provided value\n"));
     41:   nep->ncv = nep->nev;
     42:   if (nep->mpd!=PETSC_DETERMINE) PetscCall(PetscInfo(nep,"Setting mpd = nev, ignoring user-provided value\n"));
     43:   nep->mpd = nep->nev;
     44:   if (nep->max_it==PETSC_DETERMINE) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
     45:   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
     46:   PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
     47:   NEPCheckUnsupported(nep,NEP_FEATURE_REGION | NEP_FEATURE_TWOSIDED);
     48:   PetscCall(NEPAllocateSolution(nep,0));
     49:   PetscCall(NEPSetWorkVecs(nep,2));
     50:   PetscFunctionReturn(PETSC_SUCCESS);
     51: }

     53: static PetscErrorCode NEPSolve_RII(NEP nep)
     54: {
     55:   NEP_RII            *ctx = (NEP_RII*)nep->data;
     56:   Mat                T,Tp,H,A;
     57:   Vec                uu,u,r,delta,t;
     58:   PetscScalar        lambda,lambda2,sigma,a1,a2,corr;
     59:   PetscReal          nrm,resnorm=1.0,ktol=0.1,perr,rtol;
     60:   PetscBool          skip=PETSC_FALSE,lock=PETSC_FALSE;
     61:   PetscInt           inner_its,its=0;
     62:   NEP_EXT_OP         extop=NULL;
     63:   KSPConvergedReason kspreason;

     65:   PetscFunctionBegin;
     66:   /* get initial approximation of eigenvalue and eigenvector */
     67:   PetscCall(NEPGetDefaultShift(nep,&sigma));
     68:   lambda = sigma;
     69:   if (!nep->nini) {
     70:     PetscCall(BVSetRandomColumn(nep->V,0));
     71:     PetscCall(BVNormColumn(nep->V,0,NORM_2,&nrm));
     72:     PetscCall(BVScaleColumn(nep->V,0,1.0/nrm));
     73:   }
     74:   if (!ctx->ksp) PetscCall(NEPRIIGetKSP(nep,&ctx->ksp));
     75:   PetscCall(NEPDeflationInitialize(nep,nep->V,ctx->ksp,PETSC_FALSE,nep->nev,&extop));
     76:   PetscCall(NEPDeflationCreateVec(extop,&u));
     77:   PetscCall(VecDuplicate(u,&r));
     78:   PetscCall(VecDuplicate(u,&delta));
     79:   PetscCall(BVGetColumn(nep->V,0,&uu));
     80:   PetscCall(NEPDeflationCopyToExtendedVec(extop,uu,NULL,u,PETSC_FALSE));
     81:   PetscCall(BVRestoreColumn(nep->V,0,&uu));

     83:   /* prepare linear solver */
     84:   PetscCall(NEPDeflationSolveSetUp(extop,sigma));
     85:   PetscCall(KSPGetTolerances(ctx->ksp,&rtol,NULL,NULL,NULL));

     87:   PetscCall(VecCopy(u,r));
     88:   PetscCall(NEPDeflationFunctionSolve(extop,r,u));
     89:   PetscCall(VecNorm(u,NORM_2,&nrm));
     90:   PetscCall(VecScale(u,1.0/nrm));

     92:   /* Restart loop */
     93:   while (nep->reason == NEP_CONVERGED_ITERATING) {
     94:     its++;

     96:     /* Use Newton's method to compute nonlinear Rayleigh functional. Current eigenvalue
     97:        estimate as starting value. */
     98:     inner_its=0;
     99:     lambda2 = lambda;
    100:     do {
    101:       PetscCall(NEPDeflationComputeFunction(extop,lambda,&T));
    102:       PetscCall(MatMult(T,u,r));
    103:       if (!ctx->herm) {
    104:         PetscCall(NEPDeflationFunctionSolve(extop,r,delta));
    105:         PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
    106:         if (kspreason<0) PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed\n",nep->its));
    107:         t = delta;
    108:       } else t = r;
    109:       PetscCall(VecDot(t,u,&a1));
    110:       PetscCall(NEPDeflationComputeJacobian(extop,lambda,&Tp));
    111:       PetscCall(MatMult(Tp,u,r));
    112:       if (!ctx->herm) {
    113:         PetscCall(NEPDeflationFunctionSolve(extop,r,delta));
    114:         PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
    115:         if (kspreason<0) PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed\n",nep->its));
    116:         t = delta;
    117:       } else t = r;
    118:       PetscCall(VecDot(t,u,&a2));
    119:       corr = a1/a2;
    120:       lambda = lambda - corr;
    121:       inner_its++;
    122:     } while (PetscAbsScalar(corr)/PetscAbsScalar(lambda)>PETSC_SQRT_MACHINE_EPSILON && inner_its<ctx->max_inner_it);

    124:     /* form residual,  r = T(lambda)*u */
    125:     PetscCall(NEPDeflationComputeFunction(extop,lambda,&T));
    126:     PetscCall(MatMult(T,u,r));

    128:     /* convergence test */
    129:     perr = nep->errest[nep->nconv];
    130:     PetscCall(VecNorm(r,NORM_2,&resnorm));
    131:     PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
    132:     nep->eigr[nep->nconv] = lambda;
    133:     if (its>1 && (nep->errest[nep->nconv]<=nep->tol || nep->errest[nep->nconv]<=ctx->deftol)) {
    134:       if (nep->errest[nep->nconv]<=ctx->deftol && !extop->ref && nep->nconv) {
    135:         PetscCall(NEPDeflationExtractEigenpair(extop,nep->nconv,u,lambda,nep->ds));
    136:         PetscCall(NEPDeflationSetRefine(extop,PETSC_TRUE));
    137:         PetscCall(MatMult(T,u,r));
    138:         PetscCall(VecNorm(r,NORM_2,&resnorm));
    139:         PetscCall((*nep->converged)(nep,lambda,0,resnorm,&nep->errest[nep->nconv],nep->convergedctx));
    140:         if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
    141:       } else if (nep->errest[nep->nconv]<=nep->tol) lock = PETSC_TRUE;
    142:     }
    143:     if (lock) {
    144:       PetscCall(NEPDeflationSetRefine(extop,PETSC_FALSE));
    145:       nep->nconv = nep->nconv + 1;
    146:       PetscCall(NEPDeflationLocking(extop,u,lambda));
    147:       nep->its += its;
    148:       skip = PETSC_TRUE;
    149:       lock = PETSC_FALSE;
    150:     }
    151:     PetscCall((*nep->stopping)(nep,nep->its+its,nep->max_it,nep->nconv,nep->nev,&nep->reason,nep->stoppingctx));
    152:     if (!skip || nep->reason>0) PetscCall(NEPMonitor(nep,nep->its+its,nep->nconv,nep->eigr,nep->eigi,nep->errest,(nep->reason>0)?nep->nconv:nep->nconv+1));

    154:     if (nep->reason == NEP_CONVERGED_ITERATING) {
    155:       if (!skip) {
    156:         /* update preconditioner and set adaptive tolerance */
    157:         if (ctx->lag && !(its%ctx->lag) && its>=2*ctx->lag && perr && nep->errest[nep->nconv]>.5*perr) PetscCall(NEPDeflationSolveSetUp(extop,lambda2));
    158:         if (!ctx->cctol) {
    159:           ktol = PetscMax(ktol/2.0,rtol);
    160:           PetscCall(KSPSetTolerances(ctx->ksp,ktol,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
    161:         }

    163:         /* eigenvector correction: delta = T(sigma)\r */
    164:         PetscCall(NEPDeflationFunctionSolve(extop,r,delta));
    165:         PetscCall(KSPGetConvergedReason(ctx->ksp,&kspreason));
    166:         if (kspreason<0) {
    167:           PetscCall(PetscInfo(nep,"iter=%" PetscInt_FMT ", linear solve failed, stopping solve\n",nep->its));
    168:           nep->reason = NEP_DIVERGED_LINEAR_SOLVE;
    169:           break;
    170:         }

    172:         /* update eigenvector: u = u - delta */
    173:         PetscCall(VecAXPY(u,-1.0,delta));

    175:         /* normalize eigenvector */
    176:         PetscCall(VecNormalize(u,NULL));
    177:       } else {
    178:         its = -1;
    179:         PetscCall(NEPDeflationSetRandomVec(extop,u));
    180:         PetscCall(NEPDeflationSolveSetUp(extop,sigma));
    181:         PetscCall(VecCopy(u,r));
    182:         PetscCall(NEPDeflationFunctionSolve(extop,r,u));
    183:         PetscCall(VecNorm(u,NORM_2,&nrm));
    184:         PetscCall(VecScale(u,1.0/nrm));
    185:         lambda = sigma;
    186:         skip = PETSC_FALSE;
    187:       }
    188:     }
    189:   }
    190:   PetscCall(NEPDeflationGetInvariantPair(extop,NULL,&H));
    191:   PetscCall(DSSetType(nep->ds,DSNHEP));
    192:   PetscCall(DSAllocate(nep->ds,PetscMax(nep->nconv,1)));
    193:   PetscCall(DSSetDimensions(nep->ds,nep->nconv,0,nep->nconv));
    194:   PetscCall(DSGetMat(nep->ds,DS_MAT_A,&A));
    195:   PetscCall(MatCopy(H,A,SAME_NONZERO_PATTERN));
    196:   PetscCall(DSRestoreMat(nep->ds,DS_MAT_A,&A));
    197:   PetscCall(MatDestroy(&H));
    198:   PetscCall(DSSolve(nep->ds,nep->eigr,nep->eigi));
    199:   PetscCall(NEPDeflationReset(extop));
    200:   PetscCall(VecDestroy(&u));
    201:   PetscCall(VecDestroy(&r));
    202:   PetscCall(VecDestroy(&delta));
    203:   PetscFunctionReturn(PETSC_SUCCESS);
    204: }

    206: static PetscErrorCode NEPSetFromOptions_RII(NEP nep,PetscOptionItems PetscOptionsObject)
    207: {
    208:   NEP_RII        *ctx = (NEP_RII*)nep->data;
    209:   PetscBool      flg;
    210:   PetscInt       i;
    211:   PetscReal      r;

    213:   PetscFunctionBegin;
    214:   PetscOptionsHeadBegin(PetscOptionsObject,"NEP RII Options");

    216:     i = 0;
    217:     PetscCall(PetscOptionsInt("-nep_rii_max_it","Maximum number of Newton iterations for updating Rayleigh functional","NEPRIISetMaximumIterations",ctx->max_inner_it,&i,&flg));
    218:     if (flg) PetscCall(NEPRIISetMaximumIterations(nep,i));

    220:     PetscCall(PetscOptionsBool("-nep_rii_const_correction_tol","Constant correction tolerance for the linear solver","NEPRIISetConstCorrectionTol",ctx->cctol,&ctx->cctol,NULL));

    222:     PetscCall(PetscOptionsBool("-nep_rii_hermitian","Use Hermitian version of the scalar nonlinear equation","NEPRIISetHermitian",ctx->herm,&ctx->herm,NULL));

    224:     i = 0;
    225:     PetscCall(PetscOptionsInt("-nep_rii_lag_preconditioner","Interval to rebuild preconditioner","NEPRIISetLagPreconditioner",ctx->lag,&i,&flg));
    226:     if (flg) PetscCall(NEPRIISetLagPreconditioner(nep,i));

    228:     r = 0.0;
    229:     PetscCall(PetscOptionsReal("-nep_rii_deflation_threshold","Tolerance used as a threshold for including deflated eigenpairs","NEPRIISetDeflationThreshold",ctx->deftol,&r,&flg));
    230:     if (flg) PetscCall(NEPRIISetDeflationThreshold(nep,r));

    232:   PetscOptionsHeadEnd();

    234:   if (!ctx->ksp) PetscCall(NEPRIIGetKSP(nep,&ctx->ksp));
    235:   PetscCall(KSPSetFromOptions(ctx->ksp));
    236:   PetscFunctionReturn(PETSC_SUCCESS);
    237: }

    239: static PetscErrorCode NEPRIISetMaximumIterations_RII(NEP nep,PetscInt its)
    240: {
    241:   NEP_RII *ctx = (NEP_RII*)nep->data;

    243:   PetscFunctionBegin;
    244:   if (its==PETSC_DEFAULT || its==PETSC_DECIDE) ctx->max_inner_it = 10;
    245:   else {
    246:     PetscCheck(its>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations must be >0");
    247:     ctx->max_inner_it = its;
    248:   }
    249:   PetscFunctionReturn(PETSC_SUCCESS);
    250: }

    252: /*@
    253:    NEPRIISetMaximumIterations - Sets the maximum number of inner iterations to be
    254:    used in the RII solver. These are the Newton iterations related to the computation
    255:    of the nonlinear Rayleigh functional.

    257:    Logically Collective

    259:    Input Parameters:
    260: +  nep - the nonlinear eigensolver context
    261: -  its - maximum inner iterations

    263:    Options Database Key:
    264: .  -nep_rii_max_it its - set the maximum number of inner iterations

    266:    Level: advanced

    268: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetMaximumIterations()`
    269: @*/
    270: PetscErrorCode NEPRIISetMaximumIterations(NEP nep,PetscInt its)
    271: {
    272:   PetscFunctionBegin;
    275:   PetscTryMethod(nep,"NEPRIISetMaximumIterations_C",(NEP,PetscInt),(nep,its));
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }

    279: static PetscErrorCode NEPRIIGetMaximumIterations_RII(NEP nep,PetscInt *its)
    280: {
    281:   NEP_RII *ctx = (NEP_RII*)nep->data;

    283:   PetscFunctionBegin;
    284:   *its = ctx->max_inner_it;
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }

    288: /*@
    289:    NEPRIIGetMaximumIterations - Gets the maximum number of inner iterations of RII.

    291:    Not Collective

    293:    Input Parameter:
    294: .  nep - the nonlinear eigensolver context

    296:    Output Parameter:
    297: .  its - maximum inner iterations

    299:    Level: advanced

    301: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetMaximumIterations()`
    302: @*/
    303: PetscErrorCode NEPRIIGetMaximumIterations(NEP nep,PetscInt *its)
    304: {
    305:   PetscFunctionBegin;
    307:   PetscAssertPointer(its,2);
    308:   PetscUseMethod(nep,"NEPRIIGetMaximumIterations_C",(NEP,PetscInt*),(nep,its));
    309:   PetscFunctionReturn(PETSC_SUCCESS);
    310: }

    312: static PetscErrorCode NEPRIISetLagPreconditioner_RII(NEP nep,PetscInt lag)
    313: {
    314:   NEP_RII *ctx = (NEP_RII*)nep->data;

    316:   PetscFunctionBegin;
    317:   PetscCheck(lag>=0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lag must be non-negative");
    318:   ctx->lag = lag;
    319:   PetscFunctionReturn(PETSC_SUCCESS);
    320: }

    322: /*@
    323:    NEPRIISetLagPreconditioner - Determines when the preconditioner is rebuilt in the
    324:    nonlinear solve.

    326:    Logically Collective

    328:    Input Parameters:
    329: +  nep - the nonlinear eigensolver context
    330: -  lag - 0 indicates NEVER rebuild, 1 means rebuild every time the Jacobian is
    331:           computed within the nonlinear iteration, 2 means every second time
    332:           the Jacobian is built, etc.

    334:    Options Database Key:
    335: .  -nep_rii_lag_preconditioner lag - the lag value

    337:    Notes:
    338:    The default is 1.
    339:    The preconditioner is ALWAYS built in the first iteration of a nonlinear solve.

    341:    Level: intermediate

    343: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetLagPreconditioner()`
    344: @*/
    345: PetscErrorCode NEPRIISetLagPreconditioner(NEP nep,PetscInt lag)
    346: {
    347:   PetscFunctionBegin;
    350:   PetscTryMethod(nep,"NEPRIISetLagPreconditioner_C",(NEP,PetscInt),(nep,lag));
    351:   PetscFunctionReturn(PETSC_SUCCESS);
    352: }

    354: static PetscErrorCode NEPRIIGetLagPreconditioner_RII(NEP nep,PetscInt *lag)
    355: {
    356:   NEP_RII *ctx = (NEP_RII*)nep->data;

    358:   PetscFunctionBegin;
    359:   *lag = ctx->lag;
    360:   PetscFunctionReturn(PETSC_SUCCESS);
    361: }

    363: /*@
    364:    NEPRIIGetLagPreconditioner - Indicates how often the preconditioner is rebuilt.

    366:    Not Collective

    368:    Input Parameter:
    369: .  nep - the nonlinear eigensolver context

    371:    Output Parameter:
    372: .  lag - the lag parameter

    374:    Level: intermediate

    376: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetLagPreconditioner()`
    377: @*/
    378: PetscErrorCode NEPRIIGetLagPreconditioner(NEP nep,PetscInt *lag)
    379: {
    380:   PetscFunctionBegin;
    382:   PetscAssertPointer(lag,2);
    383:   PetscUseMethod(nep,"NEPRIIGetLagPreconditioner_C",(NEP,PetscInt*),(nep,lag));
    384:   PetscFunctionReturn(PETSC_SUCCESS);
    385: }

    387: static PetscErrorCode NEPRIISetConstCorrectionTol_RII(NEP nep,PetscBool cct)
    388: {
    389:   NEP_RII *ctx = (NEP_RII*)nep->data;

    391:   PetscFunctionBegin;
    392:   ctx->cctol = cct;
    393:   PetscFunctionReturn(PETSC_SUCCESS);
    394: }

    396: /*@
    397:    NEPRIISetConstCorrectionTol - Sets a flag to keep the tolerance used
    398:    in the linear solver constant.

    400:    Logically Collective

    402:    Input Parameters:
    403: +  nep - the nonlinear eigensolver context
    404: -  cct - if `PETSC_FALSE`, the `KSP` relative tolerance is set to $2^{-i}$

    406:    Options Database Key:
    407: .  -nep_rii_const_correction_tol (true|false) - set a constant or dynamic stopping criterion

    409:    Notes:
    410:    By default, an exponentially decreasing tolerance is set in the `KSP` used
    411:    within the nonlinear iteration, so that each Newton iteration requests
    412:    better accuracy than the previous one. The constant correction tolerance
    413:    flag stops this behavior.

    415:    Level: intermediate

    417: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetConstCorrectionTol()`, `NEPRIIGetKSP()`
    418: @*/
    419: PetscErrorCode NEPRIISetConstCorrectionTol(NEP nep,PetscBool cct)
    420: {
    421:   PetscFunctionBegin;
    424:   PetscTryMethod(nep,"NEPRIISetConstCorrectionTol_C",(NEP,PetscBool),(nep,cct));
    425:   PetscFunctionReturn(PETSC_SUCCESS);
    426: }

    428: static PetscErrorCode NEPRIIGetConstCorrectionTol_RII(NEP nep,PetscBool *cct)
    429: {
    430:   NEP_RII *ctx = (NEP_RII*)nep->data;

    432:   PetscFunctionBegin;
    433:   *cct = ctx->cctol;
    434:   PetscFunctionReturn(PETSC_SUCCESS);
    435: }

    437: /*@
    438:    NEPRIIGetConstCorrectionTol - Returns the constant tolerance flag.

    440:    Not Collective

    442:    Input Parameter:
    443: .  nep - the nonlinear eigensolver context

    445:    Output Parameter:
    446: .  cct - the value of the constant tolerance flag

    448:    Level: intermediate

    450: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetConstCorrectionTol()`
    451: @*/
    452: PetscErrorCode NEPRIIGetConstCorrectionTol(NEP nep,PetscBool *cct)
    453: {
    454:   PetscFunctionBegin;
    456:   PetscAssertPointer(cct,2);
    457:   PetscUseMethod(nep,"NEPRIIGetConstCorrectionTol_C",(NEP,PetscBool*),(nep,cct));
    458:   PetscFunctionReturn(PETSC_SUCCESS);
    459: }

    461: static PetscErrorCode NEPRIISetHermitian_RII(NEP nep,PetscBool herm)
    462: {
    463:   NEP_RII *ctx = (NEP_RII*)nep->data;

    465:   PetscFunctionBegin;
    466:   ctx->herm = herm;
    467:   PetscFunctionReturn(PETSC_SUCCESS);
    468: }

    470: /*@
    471:    NEPRIISetHermitian - Sets a flag to indicate if the Hermitian version of the
    472:    scalar nonlinear equation must be used by the solver.

    474:    Logically Collective

    476:    Input Parameters:
    477: +  nep  - the nonlinear eigensolver context
    478: -  herm - `PETSC_TRUE` if the Hermitian version is preferred

    480:    Options Database Key:
    481: .  -nep_rii_hermitian (true|false) - toggle the Hermitian version

    483:    Notes:
    484:    By default, the scalar nonlinear equation $x^*T(\sigma)^{-1}T(z)x=0$ is solved
    485:    at each step of the nonlinear iteration. When this flag is set the simpler
    486:    form $x^*T(z)x=0$ is used, which is supposed to be valid only for Hermitian
    487:    problems.

    489:    Level: intermediate

    491: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetHermitian()`
    492: @*/
    493: PetscErrorCode NEPRIISetHermitian(NEP nep,PetscBool herm)
    494: {
    495:   PetscFunctionBegin;
    498:   PetscTryMethod(nep,"NEPRIISetHermitian_C",(NEP,PetscBool),(nep,herm));
    499:   PetscFunctionReturn(PETSC_SUCCESS);
    500: }

    502: static PetscErrorCode NEPRIIGetHermitian_RII(NEP nep,PetscBool *herm)
    503: {
    504:   NEP_RII *ctx = (NEP_RII*)nep->data;

    506:   PetscFunctionBegin;
    507:   *herm = ctx->herm;
    508:   PetscFunctionReturn(PETSC_SUCCESS);
    509: }

    511: /*@
    512:    NEPRIIGetHermitian - Returns the flag about using the Hermitian version of
    513:    the scalar nonlinear equation.

    515:    Not Collective

    517:    Input Parameter:
    518: .  nep - the nonlinear eigensolver context

    520:    Output Parameter:
    521: .  herm - the value of the Hermitian flag

    523:    Level: intermediate

    525: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetHermitian()`
    526: @*/
    527: PetscErrorCode NEPRIIGetHermitian(NEP nep,PetscBool *herm)
    528: {
    529:   PetscFunctionBegin;
    531:   PetscAssertPointer(herm,2);
    532:   PetscUseMethod(nep,"NEPRIIGetHermitian_C",(NEP,PetscBool*),(nep,herm));
    533:   PetscFunctionReturn(PETSC_SUCCESS);
    534: }

    536: static PetscErrorCode NEPRIISetDeflationThreshold_RII(NEP nep,PetscReal deftol)
    537: {
    538:   NEP_RII *ctx = (NEP_RII*)nep->data;

    540:   PetscFunctionBegin;
    541:   ctx->deftol = deftol;
    542:   PetscFunctionReturn(PETSC_SUCCESS);
    543: }

    545: /*@
    546:    NEPRIISetDeflationThreshold - Sets the threshold value used to switch between
    547:    deflated and non-deflated iteration.

    549:    Logically Collective

    551:    Input Parameters:
    552: +  nep    - the nonlinear eigensolver context
    553: -  deftol - the threshold value

    555:    Options Database Key:
    556: .  -nep_rii_deflation_threshold deftol - set the threshold

    558:    Note:
    559:    Normally, the solver iterates on the extended problem in order to deflate
    560:    previously converged eigenpairs. If this threshold is set to a nonzero value,
    561:    then once the residual error is below this threshold the solver will
    562:    continue the iteration without deflation. The intention is to be able to
    563:    improve the current eigenpair further, despite having previous eigenpairs
    564:    with somewhat bad precision.

    566:    Level: advanced

    568: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetDeflationThreshold()`
    569: @*/
    570: PetscErrorCode NEPRIISetDeflationThreshold(NEP nep,PetscReal deftol)
    571: {
    572:   PetscFunctionBegin;
    575:   PetscTryMethod(nep,"NEPRIISetDeflationThreshold_C",(NEP,PetscReal),(nep,deftol));
    576:   PetscFunctionReturn(PETSC_SUCCESS);
    577: }

    579: static PetscErrorCode NEPRIIGetDeflationThreshold_RII(NEP nep,PetscReal *deftol)
    580: {
    581:   NEP_RII *ctx = (NEP_RII*)nep->data;

    583:   PetscFunctionBegin;
    584:   *deftol = ctx->deftol;
    585:   PetscFunctionReturn(PETSC_SUCCESS);
    586: }

    588: /*@
    589:    NEPRIIGetDeflationThreshold - Returns the threshold value that controls deflation.

    591:    Not Collective

    593:    Input Parameter:
    594: .  nep - the nonlinear eigensolver context

    596:    Output Parameter:
    597: .  deftol - the threshold

    599:    Level: advanced

    601: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetDeflationThreshold()`
    602: @*/
    603: PetscErrorCode NEPRIIGetDeflationThreshold(NEP nep,PetscReal *deftol)
    604: {
    605:   PetscFunctionBegin;
    607:   PetscAssertPointer(deftol,2);
    608:   PetscUseMethod(nep,"NEPRIIGetDeflationThreshold_C",(NEP,PetscReal*),(nep,deftol));
    609:   PetscFunctionReturn(PETSC_SUCCESS);
    610: }

    612: static PetscErrorCode NEPRIISetKSP_RII(NEP nep,KSP ksp)
    613: {
    614:   NEP_RII        *ctx = (NEP_RII*)nep->data;

    616:   PetscFunctionBegin;
    617:   PetscCall(PetscObjectReference((PetscObject)ksp));
    618:   PetscCall(KSPDestroy(&ctx->ksp));
    619:   ctx->ksp = ksp;
    620:   nep->state = NEP_STATE_INITIAL;
    621:   PetscFunctionReturn(PETSC_SUCCESS);
    622: }

    624: /*@
    625:    NEPRIISetKSP - Associate a linear solver object (`KSP`) to the nonlinear
    626:    eigenvalue solver.

    628:    Collective

    630:    Input Parameters:
    631: +  nep - the nonlinear eigensolver context
    632: -  ksp - the linear solver object

    634:    Level: advanced

    636: .seealso: [](ch:nep), `NEPRII`, `NEPRIIGetKSP()`
    637: @*/
    638: PetscErrorCode NEPRIISetKSP(NEP nep,KSP ksp)
    639: {
    640:   PetscFunctionBegin;
    643:   PetscCheckSameComm(nep,1,ksp,2);
    644:   PetscTryMethod(nep,"NEPRIISetKSP_C",(NEP,KSP),(nep,ksp));
    645:   PetscFunctionReturn(PETSC_SUCCESS);
    646: }

    648: static PetscErrorCode NEPRIIGetKSP_RII(NEP nep,KSP *ksp)
    649: {
    650:   NEP_RII        *ctx = (NEP_RII*)nep->data;

    652:   PetscFunctionBegin;
    653:   if (!ctx->ksp) {
    654:     PetscCall(KSPCreate(PetscObjectComm((PetscObject)nep),&ctx->ksp));
    655:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)nep,1));
    656:     PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)nep)->prefix));
    657:     PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"nep_rii_"));
    658:     PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)nep)->options));
    659:     PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
    660:     PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(nep->tol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
    661:   }
    662:   *ksp = ctx->ksp;
    663:   PetscFunctionReturn(PETSC_SUCCESS);
    664: }

    666: /*@
    667:    NEPRIIGetKSP - Retrieve the linear solver object (`KSP`) associated with
    668:    the nonlinear eigenvalue solver.

    670:    Collective

    672:    Input Parameter:
    673: .  nep - the nonlinear eigensolver context

    675:    Output Parameter:
    676: .  ksp - the linear solver object

    678:    Level: advanced

    680: .seealso: [](ch:nep), `NEPRII`, `NEPRIISetKSP()`
    681: @*/
    682: PetscErrorCode NEPRIIGetKSP(NEP nep,KSP *ksp)
    683: {
    684:   PetscFunctionBegin;
    686:   PetscAssertPointer(ksp,2);
    687:   PetscUseMethod(nep,"NEPRIIGetKSP_C",(NEP,KSP*),(nep,ksp));
    688:   PetscFunctionReturn(PETSC_SUCCESS);
    689: }

    691: static PetscErrorCode NEPView_RII(NEP nep,PetscViewer viewer)
    692: {
    693:   NEP_RII        *ctx = (NEP_RII*)nep->data;
    694:   PetscBool      isascii;

    696:   PetscFunctionBegin;
    697:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    698:   if (isascii) {
    699:     PetscCall(PetscViewerASCIIPrintf(viewer,"  maximum number of inner iterations: %" PetscInt_FMT "\n",ctx->max_inner_it));
    700:     if (ctx->cctol) PetscCall(PetscViewerASCIIPrintf(viewer,"  using a constant tolerance for the linear solver\n"));
    701:     if (ctx->herm) PetscCall(PetscViewerASCIIPrintf(viewer,"  using the Hermitian version of the scalar nonlinear equation\n"));
    702:     if (ctx->lag) PetscCall(PetscViewerASCIIPrintf(viewer,"  updating the preconditioner every %" PetscInt_FMT " iterations\n",ctx->lag));
    703:     if (ctx->deftol) PetscCall(PetscViewerASCIIPrintf(viewer,"  deflation threshold: %g\n",(double)ctx->deftol));
    704:     if (!ctx->ksp) PetscCall(NEPRIIGetKSP(nep,&ctx->ksp));
    705:     PetscCall(PetscViewerASCIIPushTab(viewer));
    706:     PetscCall(KSPView(ctx->ksp,viewer));
    707:     PetscCall(PetscViewerASCIIPopTab(viewer));
    708:   }
    709:   PetscFunctionReturn(PETSC_SUCCESS);
    710: }

    712: static PetscErrorCode NEPReset_RII(NEP nep)
    713: {
    714:   NEP_RII        *ctx = (NEP_RII*)nep->data;

    716:   PetscFunctionBegin;
    717:   PetscCall(KSPReset(ctx->ksp));
    718:   PetscFunctionReturn(PETSC_SUCCESS);
    719: }

    721: static PetscErrorCode NEPDestroy_RII(NEP nep)
    722: {
    723:   NEP_RII        *ctx = (NEP_RII*)nep->data;

    725:   PetscFunctionBegin;
    726:   PetscCall(KSPDestroy(&ctx->ksp));
    727:   PetscCall(PetscFree(nep->data));
    728:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NULL));
    729:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NULL));
    730:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NULL));
    731:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NULL));
    732:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NULL));
    733:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NULL));
    734:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NULL));
    735:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NULL));
    736:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NULL));
    737:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NULL));
    738:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NULL));
    739:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NULL));
    740:   PetscFunctionReturn(PETSC_SUCCESS);
    741: }

    743: /*MC
    744:    NEPRII - NEPRII = "rii" - Simple residual inverse iteration with
    745:    varying shift.

    747:    Notes:
    748:    This is the default solver, although it is very basic. Users are
    749:    advised to try other solvers.

    751:    This solver is based on the modification proposed by {cite:t}`Neu85`
    752:    of the classical residual inverse iteration method. At each step,
    753:    this method has to solve a system of linear equations. Call
    754:    `NEPRIIGetKSP()` to configure the `KSP` object used for this.
    755:    When using iterative linear solvers, the preconditioner will be updated
    756:    in each step or not, depending on `NEPRIISetLagPreconditioner()`, and
    757:    the tolerance will be constant or not depending on
    758:    `NEPRIISetConstCorrectionTol()`.

    760:    The solver incorporates deflation, so that several eigenpairs con be
    761:    computed. Details of the implementation in SLEPc can be found in
    762:    {cite:p}`Cam21`.

    764:    Level: beginner

    766: .seealso: [](ch:nep), `NEP`, `NEPType`, `NEPSetType()`, `NEPRIIGetKSP()`, `NEPRIISetLagPreconditioner()`, `NEPRIISetConstCorrectionTol()`
    767: M*/
    768: SLEPC_EXTERN PetscErrorCode NEPCreate_RII(NEP nep)
    769: {
    770:   NEP_RII        *ctx;

    772:   PetscFunctionBegin;
    773:   PetscCall(PetscNew(&ctx));
    774:   nep->data = (void*)ctx;
    775:   ctx->max_inner_it = 10;
    776:   ctx->lag          = 1;
    777:   ctx->cctol        = PETSC_FALSE;
    778:   ctx->herm         = PETSC_FALSE;
    779:   ctx->deftol       = 0.0;

    781:   nep->useds = PETSC_TRUE;

    783:   nep->ops->solve          = NEPSolve_RII;
    784:   nep->ops->setup          = NEPSetUp_RII;
    785:   nep->ops->setfromoptions = NEPSetFromOptions_RII;
    786:   nep->ops->reset          = NEPReset_RII;
    787:   nep->ops->destroy        = NEPDestroy_RII;
    788:   nep->ops->view           = NEPView_RII;
    789:   nep->ops->computevectors = NEPComputeVectors_Schur;

    791:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetMaximumIterations_C",NEPRIISetMaximumIterations_RII));
    792:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetMaximumIterations_C",NEPRIIGetMaximumIterations_RII));
    793:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetLagPreconditioner_C",NEPRIISetLagPreconditioner_RII));
    794:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetLagPreconditioner_C",NEPRIIGetLagPreconditioner_RII));
    795:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetConstCorrectionTol_C",NEPRIISetConstCorrectionTol_RII));
    796:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetConstCorrectionTol_C",NEPRIIGetConstCorrectionTol_RII));
    797:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetHermitian_C",NEPRIISetHermitian_RII));
    798:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetHermitian_C",NEPRIIGetHermitian_RII));
    799:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetDeflationThreshold_C",NEPRIISetDeflationThreshold_RII));
    800:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetDeflationThreshold_C",NEPRIIGetDeflationThreshold_RII));
    801:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIISetKSP_C",NEPRIISetKSP_RII));
    802:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPRIIGetKSP_C",NEPRIIGetKSP_RII));
    803:   PetscFunctionReturn(PETSC_SUCCESS);
    804: }