Actual source code: interpol.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: "interpol"

     13:    Method: Polynomial interpolation

     15:    Algorithm:

     17:        Uses a PEP object to solve the interpolated NEP. Currently supports
     18:        only Chebyshev interpolation on an interval.

     20:    References:

     22:        [1] C. Effenberger and D. Kresser, "Chebyshev interpolation for
     23:            nonlinear eigenvalue problems", BIT 52:933-951, 2012.
     24: */

     26: #include <slepc/private/nepimpl.h>

     28: typedef struct {
     29:   PEP       pep;
     30:   PetscReal tol;       /* tolerance for norm of polynomial coefficients */
     31:   PetscInt  maxdeg;    /* maximum degree of interpolation polynomial */
     32:   PetscInt  deg;       /* actual degree of interpolation polynomial */
     33: } NEP_INTERPOL;

     35: static PetscErrorCode NEPSetUp_Interpol(NEP nep)
     36: {
     37:   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
     38:   ST             st;
     39:   RG             rg;
     40:   PetscReal      a,b,c,d,s,tol;
     41:   PetscScalar    zero=0.0;
     42:   PetscBool      flg,istrivial,trackall;
     43:   PetscInt       its,in;

     45:   PetscFunctionBegin;
     46:   PetscCall(NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd));
     47:   PetscCheck(nep->ncv<=nep->nev+nep->mpd,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
     48:   if (nep->max_it==PETSC_DETERMINE) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
     49:   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
     50:   PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
     51:   NEPCheckUnsupported(nep,NEP_FEATURE_CALLBACK | NEP_FEATURE_STOPPING | NEP_FEATURE_TWOSIDED);

     53:   /* transfer PEP options */
     54:   if (!ctx->pep) PetscCall(NEPInterpolGetPEP(nep,&ctx->pep));
     55:   PetscCall(PEPSetBasis(ctx->pep,PEP_BASIS_CHEBYSHEV1));
     56:   PetscCall(PEPSetWhichEigenpairs(ctx->pep,PEP_TARGET_MAGNITUDE));
     57:   PetscCall(PetscObjectTypeCompare((PetscObject)ctx->pep,PEPJD,&flg));
     58:   if (!flg) {
     59:     PetscCall(PEPGetST(ctx->pep,&st));
     60:     PetscCall(STSetType(st,STSINVERT));
     61:   }
     62:   PetscCall(PEPSetDimensions(ctx->pep,nep->nev,nep->ncv,nep->mpd));
     63:   PetscCall(PEPGetTolerances(ctx->pep,&tol,&its));
     64:   if (tol==(PetscReal)PETSC_DETERMINE) tol = SlepcDefaultTol(nep->tol);
     65:   if (ctx->tol==(PetscReal)PETSC_DETERMINE) ctx->tol = tol;
     66:   if (its==PETSC_DETERMINE) its = nep->max_it;
     67:   PetscCall(PEPSetTolerances(ctx->pep,tol,its));
     68:   PetscCall(NEPGetTrackAll(nep,&trackall));
     69:   PetscCall(PEPSetTrackAll(ctx->pep,trackall));

     71:   /* transfer region options */
     72:   PetscCall(RGIsTrivial(nep->rg,&istrivial));
     73:   PetscCheck(!istrivial,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPINTERPOL requires a nontrivial region");
     74:   PetscCall(PetscObjectTypeCompare((PetscObject)nep->rg,RGINTERVAL,&flg));
     75:   PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for interval regions");
     76:   PetscCall(RGIntervalGetEndpoints(nep->rg,&a,&b,&c,&d));
     77:   PetscCheck(a>-PETSC_MAX_REAL && b<PETSC_MAX_REAL,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for bounded intervals");
     78:   PetscCall(PEPGetRG(ctx->pep,&rg));
     79:   PetscCall(RGSetType(rg,RGINTERVAL));
     80:   PetscCheck(a!=b,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for intervals on the real axis");
     81:   s = 2.0/(b-a);
     82:   c = c*s;
     83:   d = d*s;
     84:   PetscCall(RGIntervalSetEndpoints(rg,-1.0,1.0,c,d));
     85:   PetscCall(RGCheckInside(nep->rg,1,&nep->target,&zero,&in));
     86:   PetscCheck(in>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
     87:   PetscCall(PEPSetTarget(ctx->pep,(nep->target-(a+b)/2)*s));

     89:   PetscCall(NEPAllocateSolution(nep,0));
     90:   PetscFunctionReturn(PETSC_SUCCESS);
     91: }

     93: /*
     94:   Input:
     95:     d, number of nodes to compute
     96:     a,b, interval extremes
     97:   Output:
     98:     *x, array containing the d Chebyshev nodes of the interval [a,b]
     99:     *dct2, coefficients to compute a discrete cosine transformation (DCT-II)
    100: */
    101: static PetscErrorCode ChebyshevNodes(PetscInt d,PetscReal a,PetscReal b,PetscScalar *x,PetscReal *dct2)
    102: {
    103:   PetscInt  j,i;
    104:   PetscReal t;

    106:   PetscFunctionBegin;
    107:   for (j=0;j<d+1;j++) {
    108:     t = ((2*j+1)*PETSC_PI)/(2*(d+1));
    109:     x[j] = (a+b)/2.0+((b-a)/2.0)*PetscCosReal(t);
    110:     for (i=0;i<d+1;i++) dct2[j*(d+1)+i] = PetscCosReal(i*t);
    111:   }
    112:   PetscFunctionReturn(PETSC_SUCCESS);
    113: }

    115: static PetscErrorCode NEPSolve_Interpol(NEP nep)
    116: {
    117:   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
    118:   Mat            *A,*P;
    119:   PetscScalar    *x,*fx,t;
    120:   PetscReal      *cs,a,b,s,aprox,aprox0=1.0,*matnorm;
    121:   PetscInt       i,j,k,deg=ctx->maxdeg;
    122:   PetscBool      hasmnorm=PETSC_FALSE;
    123:   Vec            vr,vi=NULL;
    124:   ST             st;

    126:   PetscFunctionBegin;
    127:   PetscCall(PetscMalloc4((deg+1)*(deg+1),&cs,deg+1,&x,(deg+1)*nep->nt,&fx,nep->nt,&matnorm));
    128:   for  (j=0;j<nep->nt;j++) {
    129:     PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm));
    130:     if (!hasmnorm) break;
    131:     PetscCall(MatNorm(nep->A[j],NORM_INFINITY,matnorm+j));
    132:   }
    133:   if (!hasmnorm) for (j=0;j<nep->nt;j++) matnorm[j] = 1.0;
    134:   PetscCall(RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL));
    135:   PetscCall(ChebyshevNodes(deg,a,b,x,cs));
    136:   for (j=0;j<nep->nt;j++) {
    137:     for (i=0;i<=deg;i++) PetscCall(FNEvaluateFunction(nep->f[j],x[i],&fx[i+j*(deg+1)]));
    138:   }
    139:   /* Polynomial coefficients */
    140:   PetscCall(PetscMalloc1(deg+1,&A));
    141:   if (nep->P) PetscCall(PetscMalloc1(deg+1,&P));
    142:   ctx->deg = deg;
    143:   for (k=0;k<=deg;k++) {
    144:     PetscCall(MatDuplicate(nep->A[0],MAT_COPY_VALUES,&A[k]));
    145:     if (nep->P) PetscCall(MatDuplicate(nep->P[0],MAT_COPY_VALUES,&P[k]));
    146:     t = 0.0;
    147:     for (i=0;i<deg+1;i++) t += fx[i]*cs[i*(deg+1)+k];
    148:     t *= 2.0/(deg+1);
    149:     if (k==0) t /= 2.0;
    150:     aprox = matnorm[0]*PetscAbsScalar(t);
    151:     PetscCall(MatScale(A[k],t));
    152:     if (nep->P) PetscCall(MatScale(P[k],t));
    153:     for (j=1;j<nep->nt;j++) {
    154:       t = 0.0;
    155:       for (i=0;i<deg+1;i++) t += fx[i+j*(deg+1)]*cs[i*(deg+1)+k];
    156:       t *= 2.0/(deg+1);
    157:       if (k==0) t /= 2.0;
    158:       aprox += matnorm[j]*PetscAbsScalar(t);
    159:       PetscCall(MatAXPY(A[k],t,nep->A[j],nep->mstr));
    160:       if (nep->P) PetscCall(MatAXPY(P[k],t,nep->P[j],nep->mstrp));
    161:     }
    162:     if (k==0) aprox0 = aprox;
    163:     if (k>1 && aprox/aprox0<ctx->tol) { ctx->deg = k; deg = k; break; }
    164:   }
    165:   PetscCall(PEPSetOperators(ctx->pep,deg+1,A));
    166:   PetscCall(MatDestroyMatrices(deg+1,&A));
    167:   if (nep->P) {
    168:     PetscCall(PEPGetST(ctx->pep,&st));
    169:     PetscCall(STSetSplitPreconditioner(st,deg+1,P,nep->mstrp));
    170:     PetscCall(MatDestroyMatrices(deg+1,&P));
    171:   }
    172:   PetscCall(PetscFree4(cs,x,fx,matnorm));

    174:   /* Solve polynomial eigenproblem */
    175:   PetscCall(PEPSolve(ctx->pep));
    176:   PetscCall(PEPGetConverged(ctx->pep,&nep->nconv));
    177:   PetscCall(PEPGetIterationNumber(ctx->pep,&nep->its));
    178:   PetscCall(PEPGetConvergedReason(ctx->pep,(PEPConvergedReason*)&nep->reason));
    179:   PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
    180:   PetscCall(BVCreateVec(nep->V,&vr));
    181: #if !defined(PETSC_USE_COMPLEX)
    182:   PetscCall(VecDuplicate(vr,&vi));
    183: #endif
    184:   s = 2.0/(b-a);
    185:   for (i=0;i<nep->nconv;i++) {
    186:     PetscCall(PEPGetEigenpair(ctx->pep,i,&nep->eigr[i],&nep->eigi[i],vr,vi));
    187:     nep->eigr[i] /= s;
    188:     nep->eigr[i] += (a+b)/2.0;
    189:     nep->eigi[i] /= s;
    190:     PetscCall(BVInsertVec(nep->V,i,vr));
    191: #if !defined(PETSC_USE_COMPLEX)
    192:     if (nep->eigi[i]!=0.0) PetscCall(BVInsertVec(nep->V,++i,vi));
    193: #endif
    194:   }
    195:   PetscCall(VecDestroy(&vr));
    196:   PetscCall(VecDestroy(&vi));

    198:   nep->state = NEP_STATE_EIGENVECTORS;
    199:   PetscFunctionReturn(PETSC_SUCCESS);
    200: }

    202: static PetscErrorCode PEPMonitor_Interpol(PEP pep,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
    203: {
    204:   PetscInt       i,n;
    205:   NEP            nep = (NEP)ctx;
    206:   PetscReal      a,b,s;
    207:   ST             st;

    209:   PetscFunctionBegin;
    210:   n = PetscMin(nest,nep->ncv);
    211:   for (i=0;i<n;i++) {
    212:     nep->eigr[i]   = eigr[i];
    213:     nep->eigi[i]   = eigi[i];
    214:     nep->errest[i] = errest[i];
    215:   }
    216:   PetscCall(PEPGetST(pep,&st));
    217:   PetscCall(STBackTransform(st,n,nep->eigr,nep->eigi));
    218:   PetscCall(RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL));
    219:   s = 2.0/(b-a);
    220:   for (i=0;i<n;i++) {
    221:     nep->eigr[i] /= s;
    222:     nep->eigr[i] += (a+b)/2.0;
    223:     nep->eigi[i] /= s;
    224:   }
    225:   PetscCall(NEPMonitor(nep,its,nconv,nep->eigr,nep->eigi,nep->errest,nest));
    226:   PetscFunctionReturn(PETSC_SUCCESS);
    227: }

    229: static PetscErrorCode NEPSetFromOptions_Interpol(NEP nep,PetscOptionItems PetscOptionsObject)
    230: {
    231:   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
    232:   PetscInt       i;
    233:   PetscBool      flg1,flg2;
    234:   PetscReal      r;

    236:   PetscFunctionBegin;
    237:   PetscOptionsHeadBegin(PetscOptionsObject,"NEP Interpol Options");

    239:     PetscCall(NEPInterpolGetInterpolation(nep,&r,&i));
    240:     if (!i) i = PETSC_DETERMINE;
    241:     PetscCall(PetscOptionsInt("-nep_interpol_interpolation_degree","Maximum degree of polynomial interpolation","NEPInterpolSetInterpolation",i,&i,&flg1));
    242:     PetscCall(PetscOptionsReal("-nep_interpol_interpolation_tol","Tolerance for interpolation coefficients","NEPInterpolSetInterpolation",r,&r,&flg2));
    243:     if (flg1 || flg2) PetscCall(NEPInterpolSetInterpolation(nep,r,i));

    245:   PetscOptionsHeadEnd();

    247:   if (!ctx->pep) PetscCall(NEPInterpolGetPEP(nep,&ctx->pep));
    248:   PetscCall(PEPSetFromOptions(ctx->pep));
    249:   PetscFunctionReturn(PETSC_SUCCESS);
    250: }

    252: static PetscErrorCode NEPInterpolSetInterpolation_Interpol(NEP nep,PetscReal tol,PetscInt degree)
    253: {
    254:   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;

    256:   PetscFunctionBegin;
    257:   if (tol == (PetscReal)PETSC_DETERMINE) {
    258:     ctx->tol   = PETSC_DETERMINE;
    259:     nep->state = NEP_STATE_INITIAL;
    260:   } else if (tol != (PetscReal)PETSC_CURRENT) {
    261:     PetscCheck(tol>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
    262:     ctx->tol = tol;
    263:   }
    264:   if (degree == PETSC_DETERMINE) {
    265:     ctx->maxdeg = 0;
    266:     if (nep->state) PetscCall(NEPReset(nep));
    267:     nep->state = NEP_STATE_INITIAL;
    268:   } else if (degree != PETSC_CURRENT) {
    269:     PetscCheck(degree>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
    270:     if (ctx->maxdeg != degree) {
    271:       ctx->maxdeg = degree;
    272:       if (nep->state) PetscCall(NEPReset(nep));
    273:       nep->state = NEP_STATE_INITIAL;
    274:     }
    275:   }
    276:   PetscFunctionReturn(PETSC_SUCCESS);
    277: }

    279: /*@
    280:    NEPInterpolSetInterpolation - Sets the tolerance and maximum degree when building
    281:    the interpolation polynomial.

    283:    Collective

    285:    Input Parameters:
    286: +  nep - the nonlinear eigensolver context
    287: .  tol - tolerance to stop computing polynomial coefficients
    288: -  deg - maximum degree of interpolation

    290:    Options Database Keys:
    291: +  -nep_interpol_interpolation_tol tol    - sets the tolerance to stop computing polynomial coefficients
    292: -  -nep_interpol_interpolation_degree deg - sets the maximum degree of interpolation

    294:    Note:
    295:    `PETSC_CURRENT` can be used to preserve the current value of any of the
    296:    arguments, and `PETSC_DETERMINE` to set them to a default value.

    298:    Level: advanced

    300: .seealso: [](ch:nep), `NEPINTERPOL`, `NEPInterpolGetInterpolation()`
    301: @*/
    302: PetscErrorCode NEPInterpolSetInterpolation(NEP nep,PetscReal tol,PetscInt deg)
    303: {
    304:   PetscFunctionBegin;
    308:   PetscTryMethod(nep,"NEPInterpolSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,deg));
    309:   PetscFunctionReturn(PETSC_SUCCESS);
    310: }

    312: static PetscErrorCode NEPInterpolGetInterpolation_Interpol(NEP nep,PetscReal *tol,PetscInt *deg)
    313: {
    314:   NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;

    316:   PetscFunctionBegin;
    317:   if (tol) *tol = ctx->tol;
    318:   if (deg) *deg = ctx->maxdeg;
    319:   PetscFunctionReturn(PETSC_SUCCESS);
    320: }

    322: /*@
    323:    NEPInterpolGetInterpolation - Gets the tolerance and maximum degree when building
    324:    the interpolation polynomial.

    326:    Not Collective

    328:    Input Parameter:
    329: .  nep - the nonlinear eigensolver context

    331:    Output Parameters:
    332: +  tol - tolerance to stop computing polynomial coefficients
    333: -  deg - maximum degree of interpolation

    335:    Level: advanced

    337: .seealso: [](ch:nep), `NEPINTERPOL`, `NEPInterpolSetInterpolation()`
    338: @*/
    339: PetscErrorCode NEPInterpolGetInterpolation(NEP nep,PetscReal *tol,PetscInt *deg)
    340: {
    341:   PetscFunctionBegin;
    343:   PetscUseMethod(nep,"NEPInterpolGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,deg));
    344:   PetscFunctionReturn(PETSC_SUCCESS);
    345: }

    347: static PetscErrorCode NEPInterpolSetPEP_Interpol(NEP nep,PEP pep)
    348: {
    349:   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;

    351:   PetscFunctionBegin;
    352:   PetscCall(PetscObjectReference((PetscObject)pep));
    353:   PetscCall(PEPDestroy(&ctx->pep));
    354:   ctx->pep = pep;
    355:   nep->state = NEP_STATE_INITIAL;
    356:   PetscFunctionReturn(PETSC_SUCCESS);
    357: }

    359: /*@
    360:    NEPInterpolSetPEP - Associate a polynomial eigensolver object (`PEP`) to the
    361:    nonlinear eigenvalue solver.

    363:    Collective

    365:    Input Parameters:
    366: +  nep - the nonlinear eigensolver context
    367: -  pep - the polynomial eigensolver context

    369:    Level: advanced

    371: .seealso: [](ch:nep), `NEPINTERPOL`, `NEPInterpolGetPEP()`
    372: @*/
    373: PetscErrorCode NEPInterpolSetPEP(NEP nep,PEP pep)
    374: {
    375:   PetscFunctionBegin;
    378:   PetscCheckSameComm(nep,1,pep,2);
    379:   PetscTryMethod(nep,"NEPInterpolSetPEP_C",(NEP,PEP),(nep,pep));
    380:   PetscFunctionReturn(PETSC_SUCCESS);
    381: }

    383: static PetscErrorCode NEPInterpolGetPEP_Interpol(NEP nep,PEP *pep)
    384: {
    385:   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;

    387:   PetscFunctionBegin;
    388:   if (!ctx->pep) {
    389:     PetscCall(PEPCreate(PetscObjectComm((PetscObject)nep),&ctx->pep));
    390:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->pep,(PetscObject)nep,1));
    391:     PetscCall(PEPSetOptionsPrefix(ctx->pep,((PetscObject)nep)->prefix));
    392:     PetscCall(PEPAppendOptionsPrefix(ctx->pep,"nep_interpol_"));
    393:     PetscCall(PetscObjectSetOptions((PetscObject)ctx->pep,((PetscObject)nep)->options));
    394:     PetscCall(PEPMonitorSet(ctx->pep,PEPMonitor_Interpol,nep,NULL));
    395:   }
    396:   *pep = ctx->pep;
    397:   PetscFunctionReturn(PETSC_SUCCESS);
    398: }

    400: /*@
    401:    NEPInterpolGetPEP - Retrieve the polynomial eigensolver object (`PEP`)
    402:    associated with the nonlinear eigenvalue solver.

    404:    Collective

    406:    Input Parameter:
    407: .  nep - the nonlinear eigensolver context

    409:    Output Parameter:
    410: .  pep - the polynomial eigensolver context

    412:    Level: advanced

    414: .seealso: [](ch:nep), `NEPINTERPOL`, `NEPInterpolSetPEP()`
    415: @*/
    416: PetscErrorCode NEPInterpolGetPEP(NEP nep,PEP *pep)
    417: {
    418:   PetscFunctionBegin;
    420:   PetscAssertPointer(pep,2);
    421:   PetscUseMethod(nep,"NEPInterpolGetPEP_C",(NEP,PEP*),(nep,pep));
    422:   PetscFunctionReturn(PETSC_SUCCESS);
    423: }

    425: static PetscErrorCode NEPView_Interpol(NEP nep,PetscViewer viewer)
    426: {
    427:   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
    428:   PetscBool      isascii;

    430:   PetscFunctionBegin;
    431:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    432:   if (isascii) {
    433:     if (!ctx->pep) PetscCall(NEPInterpolGetPEP(nep,&ctx->pep));
    434:     PetscCall(PetscViewerASCIIPrintf(viewer,"  polynomial degree %" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->deg,ctx->maxdeg));
    435:     PetscCall(PetscViewerASCIIPrintf(viewer,"  tolerance for norm of polynomial coefficients %g\n",(double)ctx->tol));
    436:     PetscCall(PetscViewerASCIIPushTab(viewer));
    437:     PetscCall(PEPView(ctx->pep,viewer));
    438:     PetscCall(PetscViewerASCIIPopTab(viewer));
    439:   }
    440:   PetscFunctionReturn(PETSC_SUCCESS);
    441: }

    443: static PetscErrorCode NEPReset_Interpol(NEP nep)
    444: {
    445:   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;

    447:   PetscFunctionBegin;
    448:   PetscCall(PEPReset(ctx->pep));
    449:   PetscFunctionReturn(PETSC_SUCCESS);
    450: }

    452: static PetscErrorCode NEPDestroy_Interpol(NEP nep)
    453: {
    454:   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;

    456:   PetscFunctionBegin;
    457:   PetscCall(PEPDestroy(&ctx->pep));
    458:   PetscCall(PetscFree(nep->data));
    459:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NULL));
    460:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NULL));
    461:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NULL));
    462:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NULL));
    463:   PetscFunctionReturn(PETSC_SUCCESS);
    464: }

    466: /*MC
    467:    NEPINTERPOL - NEPINTERPOL = "interpol" - Solve the nonlinear
    468:    eigenvalue problem via polynomial interpolation.

    470:    Note:
    471:    This solver will approximate the nonlinear eigenvalue problem
    472:    with a polynomial eigenvalue problem and use `PEP` to solve it.
    473:    This is appropriate in nonlinear problems without singularities
    474:    and whose eigenvalues are all real. The degree of the polynomial
    475:    should generally be large, in which case a non-monomial representation
    476:    such as Chebyshev should be used {cite:p}`Kre14`.

    478:    Level: beginner

    480: .seealso: [](ch:nep), `NEP`, `NEPType`, `NEPSetType()`, `NEPSLPGetEPS()`, `NEPSetTwoSided()`, `NEPSLPGetEPSLeft()`
    481: M*/
    482: SLEPC_EXTERN PetscErrorCode NEPCreate_Interpol(NEP nep)
    483: {
    484:   NEP_INTERPOL   *ctx;

    486:   PetscFunctionBegin;
    487:   PetscCall(PetscNew(&ctx));
    488:   nep->data   = (void*)ctx;
    489:   ctx->maxdeg = 5;
    490:   ctx->tol    = PETSC_DETERMINE;

    492:   nep->ops->solve          = NEPSolve_Interpol;
    493:   nep->ops->setup          = NEPSetUp_Interpol;
    494:   nep->ops->setfromoptions = NEPSetFromOptions_Interpol;
    495:   nep->ops->reset          = NEPReset_Interpol;
    496:   nep->ops->destroy        = NEPDestroy_Interpol;
    497:   nep->ops->view           = NEPView_Interpol;

    499:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NEPInterpolSetInterpolation_Interpol));
    500:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NEPInterpolGetInterpolation_Interpol));
    501:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NEPInterpolSetPEP_Interpol));
    502:   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NEPInterpolGetPEP_Interpol));
    503:   PetscFunctionReturn(PETSC_SUCCESS);
    504: }