Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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"
      12             : 
      13             :    Method: Polynomial interpolation
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Uses a PEP object to solve the interpolated NEP. Currently supports
      18             :        only Chebyshev interpolation on an interval.
      19             : 
      20             :    References:
      21             : 
      22             :        [1] C. Effenberger and D. Kresser, "Chebyshev interpolation for
      23             :            nonlinear eigenvalue problems", BIT 52:933-951, 2012.
      24             : */
      25             : 
      26             : #include <slepc/private/nepimpl.h>         /*I "slepcnep.h" I*/
      27             : 
      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;
      34             : 
      35          25 : static PetscErrorCode NEPSetUp_Interpol(NEP nep)
      36             : {
      37          25 :   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
      38          25 :   ST             st;
      39          25 :   RG             rg;
      40          25 :   PetscReal      a,b,c,d,s,tol;
      41          25 :   PetscScalar    zero=0.0;
      42          25 :   PetscBool      flg,istrivial,trackall;
      43          25 :   PetscInt       its,in;
      44             : 
      45          25 :   PetscFunctionBegin;
      46          25 :   PetscCall(NEPSetDimensions_Default(nep,nep->nev,&nep->ncv,&nep->mpd));
      47          25 :   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          25 :   if (nep->max_it==PETSC_DEFAULT) nep->max_it = PetscMax(5000,2*nep->n/nep->ncv);
      49          25 :   if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
      50          25 :   PetscCheck(nep->which==NEP_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver supports only target magnitude eigenvalues");
      51          25 :   NEPCheckUnsupported(nep,NEP_FEATURE_CALLBACK | NEP_FEATURE_STOPPING | NEP_FEATURE_TWOSIDED);
      52             : 
      53             :   /* transfer PEP options */
      54          25 :   if (!ctx->pep) PetscCall(NEPInterpolGetPEP(nep,&ctx->pep));
      55          25 :   PetscCall(PEPSetBasis(ctx->pep,PEP_BASIS_CHEBYSHEV1));
      56          25 :   PetscCall(PEPSetWhichEigenpairs(ctx->pep,PEP_TARGET_MAGNITUDE));
      57          25 :   PetscCall(PetscObjectTypeCompare((PetscObject)ctx->pep,PEPJD,&flg));
      58          25 :   if (!flg) {
      59          23 :     PetscCall(PEPGetST(ctx->pep,&st));
      60          23 :     PetscCall(STSetType(st,STSINVERT));
      61             :   }
      62          25 :   PetscCall(PEPSetDimensions(ctx->pep,nep->nev,nep->ncv,nep->mpd));
      63          25 :   PetscCall(PEPGetTolerances(ctx->pep,&tol,&its));
      64          39 :   if (tol==(PetscReal)PETSC_DEFAULT) tol = SlepcDefaultTol(nep->tol);
      65          25 :   if (ctx->tol==(PetscReal)PETSC_DEFAULT) ctx->tol = tol;
      66          25 :   if (its==PETSC_DEFAULT) its = nep->max_it;
      67          25 :   PetscCall(PEPSetTolerances(ctx->pep,tol,its));
      68          25 :   PetscCall(NEPGetTrackAll(nep,&trackall));
      69          25 :   PetscCall(PEPSetTrackAll(ctx->pep,trackall));
      70             : 
      71             :   /* transfer region options */
      72          25 :   PetscCall(RGIsTrivial(nep->rg,&istrivial));
      73          25 :   PetscCheck(!istrivial,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"NEPINTERPOL requires a nontrivial region");
      74          25 :   PetscCall(PetscObjectTypeCompare((PetscObject)nep->rg,RGINTERVAL,&flg));
      75          25 :   PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for interval regions");
      76          25 :   PetscCall(RGIntervalGetEndpoints(nep->rg,&a,&b,&c,&d));
      77          25 :   PetscCheck(a>-PETSC_MAX_REAL && b<PETSC_MAX_REAL,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for bounded intervals");
      78          25 :   PetscCall(PEPGetRG(ctx->pep,&rg));
      79          25 :   PetscCall(RGSetType(rg,RGINTERVAL));
      80          25 :   PetscCheck(a!=b,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Only implemented for intervals on the real axis");
      81          25 :   s = 2.0/(b-a);
      82          25 :   c = c*s;
      83          25 :   d = d*s;
      84          25 :   PetscCall(RGIntervalSetEndpoints(rg,-1.0,1.0,c,d));
      85          25 :   PetscCall(RGCheckInside(nep->rg,1,&nep->target,&zero,&in));
      86          25 :   PetscCheck(in>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The target is not inside the target set");
      87          25 :   PetscCall(PEPSetTarget(ctx->pep,(nep->target-(a+b)/2)*s));
      88             : 
      89          25 :   PetscCall(NEPAllocateSolution(nep,0));
      90          25 :   PetscFunctionReturn(PETSC_SUCCESS);
      91             : }
      92             : 
      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          25 : static PetscErrorCode ChebyshevNodes(PetscInt d,PetscReal a,PetscReal b,PetscScalar *x,PetscReal *dct2)
     102             : {
     103          25 :   PetscInt  j,i;
     104          25 :   PetscReal t;
     105             : 
     106          25 :   PetscFunctionBegin;
     107         263 :   for (j=0;j<d+1;j++) {
     108         238 :     t = ((2*j+1)*PETSC_PI)/(2*(d+1));
     109         238 :     x[j] = (a+b)/2.0+((b-a)/2.0)*PetscCosReal(t);
     110        2870 :     for (i=0;i<d+1;i++) dct2[j*(d+1)+i] = PetscCosReal(i*t);
     111             :   }
     112          25 :   PetscFunctionReturn(PETSC_SUCCESS);
     113             : }
     114             : 
     115          25 : static PetscErrorCode NEPSolve_Interpol(NEP nep)
     116             : {
     117          25 :   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
     118          25 :   Mat            *A,*P;
     119          25 :   PetscScalar    *x,*fx,t;
     120          25 :   PetscReal      *cs,a,b,s,aprox,aprox0=1.0,*matnorm;
     121          25 :   PetscInt       i,j,k,deg=ctx->maxdeg;
     122          25 :   PetscBool      hasmnorm=PETSC_FALSE;
     123          25 :   Vec            vr,vi=NULL;
     124          25 :   ST             st;
     125             : 
     126          25 :   PetscFunctionBegin;
     127          25 :   PetscCall(PetscMalloc4((deg+1)*(deg+1),&cs,deg+1,&x,(deg+1)*nep->nt,&fx,nep->nt,&matnorm));
     128         100 :   for  (j=0;j<nep->nt;j++) {
     129          75 :     PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&hasmnorm));
     130          75 :     if (!hasmnorm) break;
     131          75 :     PetscCall(MatNorm(nep->A[j],NORM_INFINITY,matnorm+j));
     132             :   }
     133          25 :   if (!hasmnorm) for (j=0;j<nep->nt;j++) matnorm[j] = 1.0;
     134          25 :   PetscCall(RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL));
     135          25 :   PetscCall(ChebyshevNodes(deg,a,b,x,cs));
     136         100 :   for (j=0;j<nep->nt;j++) {
     137         789 :     for (i=0;i<=deg;i++) PetscCall(FNEvaluateFunction(nep->f[j],x[i],&fx[i+j*(deg+1)]));
     138             :   }
     139             :   /* Polynomial coefficients */
     140          25 :   PetscCall(PetscMalloc1(deg+1,&A));
     141          25 :   if (nep->P) PetscCall(PetscMalloc1(deg+1,&P));
     142          25 :   ctx->deg = deg;
     143         213 :   for (k=0;k<=deg;k++) {
     144         201 :     PetscCall(MatDuplicate(nep->A[0],MAT_COPY_VALUES,&A[k]));
     145         201 :     if (nep->P) PetscCall(MatDuplicate(nep->P[0],MAT_COPY_VALUES,&P[k]));
     146         201 :     t = 0.0;
     147        2591 :     for (i=0;i<deg+1;i++) t += fx[i]*cs[i*(deg+1)+k];
     148         201 :     t *= 2.0/(deg+1);
     149         201 :     if (k==0) t /= 2.0;
     150         201 :     aprox = matnorm[0]*PetscAbsScalar(t);
     151         201 :     PetscCall(MatScale(A[k],t));
     152         201 :     if (nep->P) PetscCall(MatScale(P[k],t));
     153         603 :     for (j=1;j<nep->nt;j++) {
     154             :       t = 0.0;
     155        5182 :       for (i=0;i<deg+1;i++) t += fx[i+j*(deg+1)]*cs[i*(deg+1)+k];
     156         402 :       t *= 2.0/(deg+1);
     157         402 :       if (k==0) t /= 2.0;
     158         402 :       aprox += matnorm[j]*PetscAbsScalar(t);
     159         402 :       PetscCall(MatAXPY(A[k],t,nep->A[j],nep->mstr));
     160         402 :       if (nep->P) PetscCall(MatAXPY(P[k],t,nep->P[j],nep->mstrp));
     161             :     }
     162         201 :     if (k==0) aprox0 = aprox;
     163         176 :     if (k>1 && aprox/aprox0<ctx->tol) { ctx->deg = k; deg = k; break; }
     164             :   }
     165          25 :   PetscCall(PEPSetOperators(ctx->pep,deg+1,A));
     166          25 :   PetscCall(MatDestroyMatrices(deg+1,&A));
     167          25 :   if (nep->P) {
     168           1 :     PetscCall(PEPGetST(ctx->pep,&st));
     169           1 :     PetscCall(STSetSplitPreconditioner(st,deg+1,P,nep->mstrp));
     170           1 :     PetscCall(MatDestroyMatrices(deg+1,&P));
     171             :   }
     172          25 :   PetscCall(PetscFree4(cs,x,fx,matnorm));
     173             : 
     174             :   /* Solve polynomial eigenproblem */
     175          25 :   PetscCall(PEPSolve(ctx->pep));
     176          25 :   PetscCall(PEPGetConverged(ctx->pep,&nep->nconv));
     177          25 :   PetscCall(PEPGetIterationNumber(ctx->pep,&nep->its));
     178          25 :   PetscCall(PEPGetConvergedReason(ctx->pep,(PEPConvergedReason*)&nep->reason));
     179          25 :   PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
     180          25 :   PetscCall(BVCreateVec(nep->V,&vr));
     181             : #if !defined(PETSC_USE_COMPLEX)
     182             :   PetscCall(VecDuplicate(vr,&vi));
     183             : #endif
     184          25 :   s = 2.0/(b-a);
     185         136 :   for (i=0;i<nep->nconv;i++) {
     186         111 :     PetscCall(PEPGetEigenpair(ctx->pep,i,&nep->eigr[i],&nep->eigi[i],vr,vi));
     187         111 :     nep->eigr[i] /= s;
     188         111 :     nep->eigr[i] += (a+b)/2.0;
     189         111 :     nep->eigi[i] /= s;
     190         111 :     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          25 :   PetscCall(VecDestroy(&vr));
     196          25 :   PetscCall(VecDestroy(&vi));
     197             : 
     198          25 :   nep->state = NEP_STATE_EIGENVECTORS;
     199          25 :   PetscFunctionReturn(PETSC_SUCCESS);
     200             : }
     201             : 
     202         382 : static PetscErrorCode PEPMonitor_Interpol(PEP pep,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
     203             : {
     204         382 :   PetscInt       i,n;
     205         382 :   NEP            nep = (NEP)ctx;
     206         382 :   PetscReal      a,b,s;
     207         382 :   ST             st;
     208             : 
     209         382 :   PetscFunctionBegin;
     210         382 :   n = PetscMin(nest,nep->ncv);
     211        8004 :   for (i=0;i<n;i++) {
     212        7622 :     nep->eigr[i]   = eigr[i];
     213        7622 :     nep->eigi[i]   = eigi[i];
     214        7622 :     nep->errest[i] = errest[i];
     215             :   }
     216         382 :   PetscCall(PEPGetST(pep,&st));
     217         382 :   PetscCall(STBackTransform(st,n,nep->eigr,nep->eigi));
     218         382 :   PetscCall(RGIntervalGetEndpoints(nep->rg,&a,&b,NULL,NULL));
     219         382 :   s = 2.0/(b-a);
     220        8004 :   for (i=0;i<n;i++) {
     221        7622 :     nep->eigr[i] /= s;
     222        7622 :     nep->eigr[i] += (a+b)/2.0;
     223        7622 :     nep->eigi[i] /= s;
     224             :   }
     225         382 :   PetscCall(NEPMonitor(nep,its,nconv,nep->eigr,nep->eigi,nep->errest,nest));
     226         382 :   PetscFunctionReturn(PETSC_SUCCESS);
     227             : }
     228             : 
     229          22 : static PetscErrorCode NEPSetFromOptions_Interpol(NEP nep,PetscOptionItems *PetscOptionsObject)
     230             : {
     231          22 :   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
     232          22 :   PetscInt       i;
     233          22 :   PetscBool      flg1,flg2;
     234          22 :   PetscReal      r;
     235             : 
     236          22 :   PetscFunctionBegin;
     237          22 :   PetscOptionsHeadBegin(PetscOptionsObject,"NEP Interpol Options");
     238             : 
     239          22 :     PetscCall(NEPInterpolGetInterpolation(nep,&r,&i));
     240          22 :     if (!i) i = PETSC_DEFAULT;
     241          22 :     PetscCall(PetscOptionsInt("-nep_interpol_interpolation_degree","Maximum degree of polynomial interpolation","NEPInterpolSetInterpolation",i,&i,&flg1));
     242          22 :     PetscCall(PetscOptionsReal("-nep_interpol_interpolation_tol","Tolerance for interpolation coefficients","NEPInterpolSetInterpolation",r,&r,&flg2));
     243          22 :     if (flg1 || flg2) PetscCall(NEPInterpolSetInterpolation(nep,r,i));
     244             : 
     245          22 :   PetscOptionsHeadEnd();
     246             : 
     247          22 :   if (!ctx->pep) PetscCall(NEPInterpolGetPEP(nep,&ctx->pep));
     248          22 :   PetscCall(PEPSetFromOptions(ctx->pep));
     249          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     250             : }
     251             : 
     252          14 : static PetscErrorCode NEPInterpolSetInterpolation_Interpol(NEP nep,PetscReal tol,PetscInt degree)
     253             : {
     254          14 :   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
     255             : 
     256          14 :   PetscFunctionBegin;
     257          14 :   if (tol == (PetscReal)PETSC_DEFAULT) {
     258          14 :     ctx->tol   = PETSC_DEFAULT;
     259          14 :     nep->state = NEP_STATE_INITIAL;
     260             :   } else {
     261           0 :     PetscCheck(tol>0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of tol. Must be > 0");
     262           0 :     ctx->tol = tol;
     263             :   }
     264          14 :   if (degree == PETSC_DEFAULT || degree == PETSC_DECIDE) {
     265           0 :     ctx->maxdeg = 0;
     266           0 :     if (nep->state) PetscCall(NEPReset(nep));
     267           0 :     nep->state = NEP_STATE_INITIAL;
     268             :   } else {
     269          14 :     PetscCheck(degree>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of degree. Must be > 0");
     270          14 :     if (ctx->maxdeg != degree) {
     271          14 :       ctx->maxdeg = degree;
     272          14 :       if (nep->state) PetscCall(NEPReset(nep));
     273          14 :       nep->state = NEP_STATE_INITIAL;
     274             :     }
     275             :   }
     276          14 :   PetscFunctionReturn(PETSC_SUCCESS);
     277             : }
     278             : 
     279             : /*@
     280             :    NEPInterpolSetInterpolation - Sets the tolerance and maximum degree when building
     281             :    the interpolation polynomial.
     282             : 
     283             :    Collective
     284             : 
     285             :    Input Parameters:
     286             : +  nep - nonlinear eigenvalue solver
     287             : .  tol - tolerance to stop computing polynomial coefficients
     288             : -  deg - maximum degree of interpolation
     289             : 
     290             :    Options Database Key:
     291             : +  -nep_interpol_interpolation_tol <tol> - Sets the tolerance to stop computing polynomial coefficients
     292             : -  -nep_interpol_interpolation_degree <degree> - Sets the maximum degree of interpolation
     293             : 
     294             :    Notes:
     295             :    Use PETSC_DEFAULT for either argument to assign a reasonably good value.
     296             : 
     297             :    Level: advanced
     298             : 
     299             : .seealso: NEPInterpolGetInterpolation()
     300             : @*/
     301          14 : PetscErrorCode NEPInterpolSetInterpolation(NEP nep,PetscReal tol,PetscInt deg)
     302             : {
     303          14 :   PetscFunctionBegin;
     304          14 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     305          56 :   PetscValidLogicalCollectiveReal(nep,tol,2);
     306          56 :   PetscValidLogicalCollectiveInt(nep,deg,3);
     307          14 :   PetscTryMethod(nep,"NEPInterpolSetInterpolation_C",(NEP,PetscReal,PetscInt),(nep,tol,deg));
     308          14 :   PetscFunctionReturn(PETSC_SUCCESS);
     309             : }
     310             : 
     311          24 : static PetscErrorCode NEPInterpolGetInterpolation_Interpol(NEP nep,PetscReal *tol,PetscInt *deg)
     312             : {
     313          24 :   NEP_INTERPOL *ctx = (NEP_INTERPOL*)nep->data;
     314             : 
     315          24 :   PetscFunctionBegin;
     316          24 :   if (tol) *tol = ctx->tol;
     317          24 :   if (deg) *deg = ctx->maxdeg;
     318          24 :   PetscFunctionReturn(PETSC_SUCCESS);
     319             : }
     320             : 
     321             : /*@
     322             :    NEPInterpolGetInterpolation - Gets the tolerance and maximum degree when building
     323             :    the interpolation polynomial.
     324             : 
     325             :    Not Collective
     326             : 
     327             :    Input Parameter:
     328             : .  nep - nonlinear eigenvalue solver
     329             : 
     330             :    Output Parameters:
     331             : +  tol - tolerance to stop computing polynomial coefficients
     332             : -  deg - maximum degree of interpolation
     333             : 
     334             :    Level: advanced
     335             : 
     336             : .seealso: NEPInterpolSetInterpolation()
     337             : @*/
     338          24 : PetscErrorCode NEPInterpolGetInterpolation(NEP nep,PetscReal *tol,PetscInt *deg)
     339             : {
     340          24 :   PetscFunctionBegin;
     341          24 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     342          24 :   PetscUseMethod(nep,"NEPInterpolGetInterpolation_C",(NEP,PetscReal*,PetscInt*),(nep,tol,deg));
     343          24 :   PetscFunctionReturn(PETSC_SUCCESS);
     344             : }
     345             : 
     346           2 : static PetscErrorCode NEPInterpolSetPEP_Interpol(NEP nep,PEP pep)
     347             : {
     348           2 :   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
     349             : 
     350           2 :   PetscFunctionBegin;
     351           2 :   PetscCall(PetscObjectReference((PetscObject)pep));
     352           2 :   PetscCall(PEPDestroy(&ctx->pep));
     353           2 :   ctx->pep = pep;
     354           2 :   nep->state = NEP_STATE_INITIAL;
     355           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     356             : }
     357             : 
     358             : /*@
     359             :    NEPInterpolSetPEP - Associate a polynomial eigensolver object (PEP) to the
     360             :    nonlinear eigenvalue solver.
     361             : 
     362             :    Collective
     363             : 
     364             :    Input Parameters:
     365             : +  nep - nonlinear eigenvalue solver
     366             : -  pep - the polynomial eigensolver object
     367             : 
     368             :    Level: advanced
     369             : 
     370             : .seealso: NEPInterpolGetPEP()
     371             : @*/
     372           2 : PetscErrorCode NEPInterpolSetPEP(NEP nep,PEP pep)
     373             : {
     374           2 :   PetscFunctionBegin;
     375           2 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     376           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,2);
     377           2 :   PetscCheckSameComm(nep,1,pep,2);
     378           2 :   PetscTryMethod(nep,"NEPInterpolSetPEP_C",(NEP,PEP),(nep,pep));
     379           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     380             : }
     381             : 
     382          20 : static PetscErrorCode NEPInterpolGetPEP_Interpol(NEP nep,PEP *pep)
     383             : {
     384          20 :   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
     385             : 
     386          20 :   PetscFunctionBegin;
     387          20 :   if (!ctx->pep) {
     388          20 :     PetscCall(PEPCreate(PetscObjectComm((PetscObject)nep),&ctx->pep));
     389          20 :     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->pep,(PetscObject)nep,1));
     390          20 :     PetscCall(PEPSetOptionsPrefix(ctx->pep,((PetscObject)nep)->prefix));
     391          20 :     PetscCall(PEPAppendOptionsPrefix(ctx->pep,"nep_interpol_"));
     392          20 :     PetscCall(PetscObjectSetOptions((PetscObject)ctx->pep,((PetscObject)nep)->options));
     393          20 :     PetscCall(PEPMonitorSet(ctx->pep,PEPMonitor_Interpol,nep,NULL));
     394             :   }
     395          20 :   *pep = ctx->pep;
     396          20 :   PetscFunctionReturn(PETSC_SUCCESS);
     397             : }
     398             : 
     399             : /*@
     400             :    NEPInterpolGetPEP - Retrieve the polynomial eigensolver object (PEP)
     401             :    associated with the nonlinear eigenvalue solver.
     402             : 
     403             :    Collective
     404             : 
     405             :    Input Parameter:
     406             : .  nep - nonlinear eigenvalue solver
     407             : 
     408             :    Output Parameter:
     409             : .  pep - the polynomial eigensolver object
     410             : 
     411             :    Level: advanced
     412             : 
     413             : .seealso: NEPInterpolSetPEP()
     414             : @*/
     415          20 : PetscErrorCode NEPInterpolGetPEP(NEP nep,PEP *pep)
     416             : {
     417          20 :   PetscFunctionBegin;
     418          20 :   PetscValidHeaderSpecific(nep,NEP_CLASSID,1);
     419          20 :   PetscAssertPointer(pep,2);
     420          20 :   PetscUseMethod(nep,"NEPInterpolGetPEP_C",(NEP,PEP*),(nep,pep));
     421          20 :   PetscFunctionReturn(PETSC_SUCCESS);
     422             : }
     423             : 
     424           0 : static PetscErrorCode NEPView_Interpol(NEP nep,PetscViewer viewer)
     425             : {
     426           0 :   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
     427           0 :   PetscBool      isascii;
     428             : 
     429           0 :   PetscFunctionBegin;
     430           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     431           0 :   if (isascii) {
     432           0 :     if (!ctx->pep) PetscCall(NEPInterpolGetPEP(nep,&ctx->pep));
     433           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  polynomial degree %" PetscInt_FMT ", max=%" PetscInt_FMT "\n",ctx->deg,ctx->maxdeg));
     434           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  tolerance for norm of polynomial coefficients %g\n",(double)ctx->tol));
     435           0 :     PetscCall(PetscViewerASCIIPushTab(viewer));
     436           0 :     PetscCall(PEPView(ctx->pep,viewer));
     437           0 :     PetscCall(PetscViewerASCIIPopTab(viewer));
     438             :   }
     439           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     440             : }
     441             : 
     442          25 : static PetscErrorCode NEPReset_Interpol(NEP nep)
     443             : {
     444          25 :   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
     445             : 
     446          25 :   PetscFunctionBegin;
     447          25 :   PetscCall(PEPReset(ctx->pep));
     448          25 :   PetscFunctionReturn(PETSC_SUCCESS);
     449             : }
     450             : 
     451          22 : static PetscErrorCode NEPDestroy_Interpol(NEP nep)
     452             : {
     453          22 :   NEP_INTERPOL   *ctx = (NEP_INTERPOL*)nep->data;
     454             : 
     455          22 :   PetscFunctionBegin;
     456          22 :   PetscCall(PEPDestroy(&ctx->pep));
     457          22 :   PetscCall(PetscFree(nep->data));
     458          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NULL));
     459          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NULL));
     460          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NULL));
     461          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NULL));
     462          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     463             : }
     464             : 
     465          22 : SLEPC_EXTERN PetscErrorCode NEPCreate_Interpol(NEP nep)
     466             : {
     467          22 :   NEP_INTERPOL   *ctx;
     468             : 
     469          22 :   PetscFunctionBegin;
     470          22 :   PetscCall(PetscNew(&ctx));
     471          22 :   nep->data   = (void*)ctx;
     472          22 :   ctx->maxdeg = 5;
     473          22 :   ctx->tol    = PETSC_DEFAULT;
     474             : 
     475          22 :   nep->ops->solve          = NEPSolve_Interpol;
     476          22 :   nep->ops->setup          = NEPSetUp_Interpol;
     477          22 :   nep->ops->setfromoptions = NEPSetFromOptions_Interpol;
     478          22 :   nep->ops->reset          = NEPReset_Interpol;
     479          22 :   nep->ops->destroy        = NEPDestroy_Interpol;
     480          22 :   nep->ops->view           = NEPView_Interpol;
     481             : 
     482          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetInterpolation_C",NEPInterpolSetInterpolation_Interpol));
     483          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetInterpolation_C",NEPInterpolGetInterpolation_Interpol));
     484          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolSetPEP_C",NEPInterpolSetPEP_Interpol));
     485          22 :   PetscCall(PetscObjectComposeFunction((PetscObject)nep,"NEPInterpolGetPEP_C",NEPInterpolGetPEP_Interpol));
     486          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     487             : }