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