Actual source code: dsnep.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: */

     11: #include <slepc/private/dsimpl.h>
     12: #include <slepcblaslapack.h>

     14: typedef struct {
     15:   PetscInt              nf;                     /* number of functions in f[] */
     16:   FN                    f[DS_NUM_EXTRA];        /* functions defining the nonlinear operator */
     17:   PetscInt              max_mid;                /* maximum minimality index */
     18:   PetscInt              nnod;                   /* number of nodes for quadrature rules */
     19:   PetscInt              spls;                   /* number of sampling columns for quadrature rules */
     20:   PetscInt              Nit;                    /* number of refinement iterations */
     21:   PetscReal             rtol;                   /* tolerance of Newton refinement */
     22:   RG                    rg;                     /* region for contour integral */
     23:   PetscLayout           map;                    /* used to distribute work among MPI processes */
     24:   DSNEPMatrixFunctionFn *computematrix;         /* user-provided compute matrix function */
     25:   void                  *computematrixctx;      /* context for the compute matrix function */
     26:   PetscCtxDestroyFn     *computematrixdestroy;  /* context destroy function */
     27: } DS_NEP;

     29: /*
     30:    DSNEPComputeMatrix - Build the matrix associated with a nonlinear operator
     31:    T(lambda) or its derivative T'(lambda), given the parameter lambda, where
     32:    T(lambda) = sum_i E_i*f_i(lambda). The result is written in mat.
     33: */
     34: static PetscErrorCode DSNEPComputeMatrix(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat)
     35: {
     36:   DS_NEP            *ctx = (DS_NEP*)ds->data;
     37:   PetscScalar       *T,alpha;
     38:   const PetscScalar *E;
     39:   PetscInt          i,ld,n;
     40:   PetscBLASInt      k,inc=1;

     42:   PetscFunctionBegin;
     43:   PetscCall(PetscLogEventBegin(DS_Other,ds,0,0,0));
     44:   if (ctx->computematrix) PetscCall((*ctx->computematrix)(ds,lambda,deriv,mat,ctx->computematrixctx));
     45:   else {
     46:     PetscCall(DSGetDimensions(ds,&n,NULL,NULL,NULL));
     47:     PetscCall(DSGetLeadingDimension(ds,&ld));
     48:     PetscCall(PetscBLASIntCast(ld*n,&k));
     49:     PetscCall(MatDenseGetArray(ds->omat[mat],&T));
     50:     PetscCall(PetscArrayzero(T,k));
     51:     for (i=0;i<ctx->nf;i++) {
     52:       if (deriv) PetscCall(FNEvaluateDerivative(ctx->f[i],lambda,&alpha));
     53:       else PetscCall(FNEvaluateFunction(ctx->f[i],lambda,&alpha));
     54:       PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&E));
     55:       PetscCallBLAS("BLASaxpy",BLASaxpy_(&k,&alpha,E,&inc,T,&inc));
     56:       PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&E));
     57:     }
     58:     PetscCall(MatDenseRestoreArray(ds->omat[mat],&T));
     59:   }
     60:   PetscCall(PetscLogEventEnd(DS_Other,ds,0,0,0));
     61:   PetscFunctionReturn(PETSC_SUCCESS);
     62: }

     64: static PetscErrorCode DSAllocate_NEP(DS ds,PetscInt ld)
     65: {
     66:   DS_NEP         *ctx = (DS_NEP*)ds->data;
     67:   PetscInt       i;

     69:   PetscFunctionBegin;
     70:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
     71:   for (i=0;i<ctx->nf;i++) PetscCall(DSAllocateMat_Private(ds,DSMatExtra[i]));
     72:   PetscCall(PetscFree(ds->perm));
     73:   PetscCall(PetscMalloc1(ld*ctx->max_mid,&ds->perm));
     74:   PetscFunctionReturn(PETSC_SUCCESS);
     75: }

     77: static PetscErrorCode DSView_NEP(DS ds,PetscViewer viewer)
     78: {
     79:   DS_NEP            *ctx = (DS_NEP*)ds->data;
     80:   PetscViewerFormat format;
     81:   PetscInt          i;
     82:   const char        *methodname[] = {
     83:                      "Successive Linear Problems",
     84:                      "Contour Integral"
     85:   };
     86:   const int         nmeth=PETSC_STATIC_ARRAY_LENGTH(methodname);

     88:   PetscFunctionBegin;
     89:   PetscCall(PetscViewerGetFormat(viewer,&format));
     90:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
     91:     if (ds->method<nmeth) PetscCall(PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]));
     92: #if defined(PETSC_USE_COMPLEX)
     93:     if (ds->method==1) {  /* contour integral method */
     94:       PetscCall(PetscViewerASCIIPrintf(viewer,"number of integration points: %" PetscInt_FMT "\n",ctx->nnod));
     95:       PetscCall(PetscViewerASCIIPrintf(viewer,"maximum minimality index: %" PetscInt_FMT "\n",ctx->max_mid));
     96:       if (ctx->spls) PetscCall(PetscViewerASCIIPrintf(viewer,"number of sampling columns for quadrature: %" PetscInt_FMT "\n",ctx->spls));
     97:       if (ctx->Nit) PetscCall(PetscViewerASCIIPrintf(viewer,"doing iterative refinement (%" PetscInt_FMT " its, tolerance %g)\n",ctx->Nit,(double)ctx->rtol));
     98:       PetscCall(RGView(ctx->rg,viewer));
     99:     }
    100: #endif
    101:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscCall(PetscViewerASCIIPrintf(viewer,"number of functions: %" PetscInt_FMT "\n",ctx->nf));
    102:     PetscFunctionReturn(PETSC_SUCCESS);
    103:   }
    104:   for (i=0;i<ctx->nf;i++) {
    105:     PetscCall(FNView(ctx->f[i],viewer));
    106:     PetscCall(DSViewMat(ds,viewer,DSMatExtra[i]));
    107:   }
    108:   if (ds->state>DS_STATE_INTERMEDIATE) PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
    109:   PetscFunctionReturn(PETSC_SUCCESS);
    110: }

    112: static PetscErrorCode DSVectors_NEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
    113: {
    114:   PetscFunctionBegin;
    115:   PetscCheck(!rnorm,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
    116:   switch (mat) {
    117:     case DS_MAT_X:
    118:       break;
    119:     case DS_MAT_Y:
    120:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
    121:     default:
    122:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
    123:   }
    124:   PetscFunctionReturn(PETSC_SUCCESS);
    125: }

    127: static PetscErrorCode DSSort_NEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *dummy)
    128: {
    129:   DS_NEP         *ctx = (DS_NEP*)ds->data;
    130:   PetscInt       n,l,i,*perm,lds;
    131:   PetscScalar    *Q;

    133:   PetscFunctionBegin;
    134:   if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
    135:   if (!ds->method) PetscFunctionReturn(PETSC_SUCCESS);  /* SLP computes just one eigenvalue */
    136:   n = ds->n*ctx->max_mid;
    137:   lds = ds->ld*ctx->max_mid;
    138:   l = ds->l;
    139:   perm = ds->perm;
    140:   for (i=0;i<n;i++) perm[i] = i;
    141:   if (rr) PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
    142:   else PetscCall(DSSortEigenvalues_Private(ds,wr,NULL,perm,PETSC_FALSE));
    143:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
    144:   for (i=l;i<ds->t;i++) Q[i+i*lds] = wr[perm[i]];
    145:   for (i=l;i<ds->t;i++) wr[i] = Q[i+i*lds];
    146:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
    147:   /* n != ds->n */
    148:   PetscCall(DSPermuteColumns_Private(ds,0,ds->t,ds->n,DS_MAT_X,perm));
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: #if defined(SLEPC_MISSING_LAPACK_GGEV3)
    153: #define LAPGEEV "ggev"
    154: #else
    155: #define LAPGEEV "ggev3"
    156: #endif

    158: static PetscErrorCode DSSolve_NEP_SLP(DS ds,PetscScalar *wr,PetscScalar *wi)
    159: {
    160:   PetscScalar    *A,*B,*W,*X,*work,*alpha,*beta,a;
    161:   PetscScalar    sigma,lambda,mu,re,re2,sone=1.0,szero=0.0;
    162:   PetscBLASInt   info,n,ld,lwork,one=1,zero=0;
    163:   PetscInt       it,pos,j,maxit=100,result;
    164:   PetscReal      norm,tol,done=1.0;
    165: #if !defined(PETSC_USE_COMPLEX)
    166:   PetscReal      *alphai,im,im2;
    167: #endif

    169:   PetscFunctionBegin;
    170:   PetscCall(PetscBLASIntCast(ds->n,&n));
    171:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
    172:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
    173:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
    174:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
    175:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
    176:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
    177:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
    178:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));

    180:   /* workspace query and memory allocation */
    181:   lwork = -1;
    182: #if defined(PETSC_USE_COMPLEX)
    183:   PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&n,A,&ld,B,&ld,NULL,NULL,NULL,&ld,W,&ld,&a,&lwork,NULL,&info));
    184:   PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
    185:   PetscCall(DSAllocateWork_Private(ds,lwork+2*ds->n,8*ds->n,0));
    186:   alpha = ds->work;
    187:   beta  = ds->work + ds->n;
    188:   work  = ds->work + 2*ds->n;
    189: #else
    190:   PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&n,A,&ld,B,&ld,NULL,NULL,NULL,NULL,&ld,W,&ld,&a,&lwork,&info));
    191:   PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
    192:   PetscCall(DSAllocateWork_Private(ds,lwork+3*ds->n,0,0));
    193:   alpha  = ds->work;
    194:   beta   = ds->work + ds->n;
    195:   alphai = ds->work + 2*ds->n;
    196:   work   = ds->work + 3*ds->n;
    197: #endif

    199:   sigma = 0.0;
    200:   if (ds->sc->comparison==SlepcCompareTargetMagnitude || ds->sc->comparison==SlepcCompareTargetReal) sigma = *(PetscScalar*)ds->sc->comparisonctx;
    201:   lambda = sigma;
    202:   tol = n*PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);

    204:   for (it=0;it<maxit;it++) {

    206:     /* evaluate T and T' */
    207:     PetscCall(DSNEPComputeMatrix(ds,lambda,PETSC_FALSE,DS_MAT_A));
    208:     if (it) {
    209:       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,A,&ld,X,&one,&szero,X+ld,&one));
    210:       norm = BLASnrm2_(&n,X+ld,&one);
    211:       if (norm/PetscAbsScalar(lambda)<=tol) break;
    212:     }
    213:     PetscCall(DSNEPComputeMatrix(ds,lambda,PETSC_TRUE,DS_MAT_B));

    215:     /* compute eigenvalue correction mu and eigenvector u */
    216: #if defined(PETSC_USE_COMPLEX)
    217:     PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&n,A,&ld,B,&ld,alpha,beta,NULL,&ld,W,&ld,work,&lwork,ds->rwork,&info));
    218: #else
    219:     PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&n,A,&ld,B,&ld,alpha,alphai,beta,NULL,&ld,W,&ld,work,&lwork,&info));
    220: #endif
    221:     SlepcCheckLapackInfo(LAPGEEV,info);

    223:     /* find smallest eigenvalue */
    224:     j = 0;
    225:     if (beta[j]==0.0) re = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
    226:     else re = alpha[j]/beta[j];
    227: #if !defined(PETSC_USE_COMPLEX)
    228:     if (beta[j]==0.0) im = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
    229:     else im = alphai[j]/beta[j];
    230: #endif
    231:     pos = 0;
    232:     for (j=1;j<n;j++) {
    233:       if (beta[j]==0.0) re2 = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
    234:       else re2 = alpha[j]/beta[j];
    235: #if !defined(PETSC_USE_COMPLEX)
    236:       if (beta[j]==0.0) im2 = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
    237:       else im2 = alphai[j]/beta[j];
    238:       PetscCall(SlepcCompareSmallestMagnitude(re,im,re2,im2,&result,NULL));
    239: #else
    240:       PetscCall(SlepcCompareSmallestMagnitude(re,0.0,re2,0.0,&result,NULL));
    241: #endif
    242:       if (result > 0) {
    243:         re = re2;
    244: #if !defined(PETSC_USE_COMPLEX)
    245:         im = im2;
    246: #endif
    247:         pos = j;
    248:       }
    249:     }

    251: #if !defined(PETSC_USE_COMPLEX)
    252:     PetscCheck(im==0.0,PETSC_COMM_SELF,PETSC_ERR_SUP,"DSNEP found a complex eigenvalue; try rerunning with complex scalars");
    253: #endif
    254:     mu = alpha[pos]/beta[pos];
    255:     PetscCall(PetscArraycpy(X,W+pos*ld,n));
    256:     norm = BLASnrm2_(&n,X,&one);
    257:     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&one,X,&n,&info));
    258:     SlepcCheckLapackInfo("lascl",info);

    260:     /* correct eigenvalue approximation */
    261:     lambda = lambda - mu;
    262:   }
    263:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
    264:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
    265:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
    266:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));

    268:   PetscCheck(it<maxit,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"DSNEP did not converge");
    269:   ds->t = 1;
    270:   wr[0] = lambda;
    271:   if (wi) wi[0] = 0.0;
    272:   PetscFunctionReturn(PETSC_SUCCESS);
    273: }

    275: #if defined(PETSC_USE_COMPLEX)
    276: /*
    277:   Newton refinement for eigenpairs computed with contour integral.
    278:   k  - number of eigenpairs to refine
    279:   wr - eigenvalues (eigenvectors are stored in DS_MAT_X)
    280: */
    281: static PetscErrorCode DSNEPNewtonRefine(DS ds,PetscInt k,PetscScalar *wr)
    282: {
    283:   DS_NEP         *ctx = (DS_NEP*)ds->data;
    284:   PetscScalar    *X,*W,*U,*R,sone=1.0,szero=0.0;
    285:   PetscReal      norm;
    286:   PetscInt       i,j,ii,nwu=0,*p,jstart=0,jend=k;
    287:   const PetscInt *range;
    288:   PetscBLASInt   n,*perm,info,ld,one=1,n1;
    289:   PetscMPIInt    len,size,root;
    290:   PetscLayout    map;

    292:   PetscFunctionBegin;
    293:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
    294:   PetscCall(PetscBLASIntCast(ds->n,&n));
    295:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
    296:   n1 = n+1;
    297:   p  = ds->perm;
    298:   PetscCall(PetscArrayzero(p,k));
    299:   PetscCall(DSAllocateWork_Private(ds,(n+2)*(n+1),0,n+1));
    300:   U    = ds->work+nwu;    nwu += (n+1)*(n+1);
    301:   R    = ds->work+nwu;    /*nwu += n+1;*/
    302:   perm = ds->iwork;
    303:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {
    304:     PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)ds),PETSC_DECIDE,k,1,&map));
    305:     PetscCall(PetscLayoutGetRange(map,&jstart,&jend));
    306:   }
    307:   for (ii=0;ii<ctx->Nit;ii++) {
    308:     for (j=jstart;j<jend;j++) {
    309:       if (p[j]<2) {
    310:         PetscCall(DSNEPComputeMatrix(ds,wr[j],PETSC_FALSE,DS_MAT_W));
    311:         PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
    312:         PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,W,&ld,X+ld*j,&one,&szero,R,&one));
    313:         PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
    314:         norm = BLASnrm2_(&n,R,&one);
    315:         if (norm/PetscAbsScalar(wr[j]) > ctx->rtol) {
    316:           PetscCall(PetscInfo(NULL,"Refining eigenpair %" PetscInt_FMT ", residual=%g\n",j,(double)(norm/PetscAbsScalar(wr[j]))));
    317:           p[j] = 1;
    318:           R[n] = 0.0;
    319:           PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
    320:           for (i=0;i<n;i++) {
    321:             PetscCall(PetscArraycpy(U+i*n1,W+i*ld,n));
    322:             U[n+i*n1] = PetscConj(X[j*ld+i]);
    323:           }
    324:           PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
    325:           U[n+n*n1] = 0.0;
    326:           PetscCall(DSNEPComputeMatrix(ds,wr[j],PETSC_TRUE,DS_MAT_W));
    327:           PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
    328:           PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,W,&ld,X+ld*j,&one,&szero,U+n*(n+1),&one));
    329:           PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
    330:           /* solve system  */
    331:           PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n1,&n1,U,&n1,perm,&info));
    332:           SlepcCheckLapackInfo("getrf",info);
    333:           PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n1,&one,U,&n1,perm,R,&n1,&info));
    334:           SlepcCheckLapackInfo("getrs",info);
    335:           wr[j] -= R[n];
    336:           for (i=0;i<n;i++) X[j*ld+i] -= R[i];
    337:           /* normalization */
    338:           norm = BLASnrm2_(&n,X+ld*j,&one);
    339:           for (i=0;i<n;i++) X[ld*j+i] /= norm;
    340:         } else p[j] = 2;
    341:       }
    342:     }
    343:   }
    344:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {  /* communicate results */
    345:     PetscCall(PetscMPIIntCast(k,&len));
    346:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,p,len,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)ds)));
    347:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size));
    348:     PetscCall(PetscLayoutGetRanges(map,&range));
    349:     for (j=0;j<k;j++) {
    350:       if (p[j]) {  /* j-th eigenpair has been refined */
    351:         for (root=0;root<size;root++) if (range[root+1]>j) break;
    352:         PetscCall(PetscMPIIntCast(1,&len));
    353:         PetscCallMPI(MPI_Bcast(wr+j,len,MPIU_SCALAR,root,PetscObjectComm((PetscObject)ds)));
    354:         PetscCall(PetscMPIIntCast(n,&len));
    355:         PetscCallMPI(MPI_Bcast(X+ld*j,len,MPIU_SCALAR,root,PetscObjectComm((PetscObject)ds)));
    356:       }
    357:     }
    358:     PetscCall(PetscLayoutDestroy(&map));
    359:   }
    360:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
    361:   PetscFunctionReturn(PETSC_SUCCESS);
    362: }

    364: PetscErrorCode DSSolve_NEP_Contour(DS ds,PetscScalar *wr,PetscScalar *wi)
    365: {
    366:   DS_NEP         *ctx = (DS_NEP*)ds->data;
    367:   PetscScalar    *alpha,*beta,*Q,*Z,*X,*U,*V,*W,*work,*Rc,*R,*w,*z,*zn,*S;
    368:   PetscScalar    sone=1.0,szero=0.0,center,a;
    369:   PetscReal      *rwork,norm,radius,vscale,rgscale,*sigma;
    370:   PetscBLASInt   info,n,*perm,p,pp,ld,lwork,k_,rk_,colA,rowA,one=1;
    371:   PetscInt       mid,lds,nnod=ctx->nnod,k,i,ii,jj,j,s,off,rk,nwu=0,nw,lrwork,*inside,kstart=0,kend=nnod;
    372:   PetscMPIInt    len;
    373:   PetscBool      isellipse;
    374:   PetscRandom    rand;

    376:   PetscFunctionBegin;
    377:   PetscCheck(ctx->rg,PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"The contour solver requires a region passed with DSNEPSetRG()");
    378:   /* Contour parameters */
    379:   PetscCall(PetscObjectTypeCompare((PetscObject)ctx->rg,RGELLIPSE,&isellipse));
    380:   PetscCheck(isellipse,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Region must be Ellipse");
    381:   PetscCall(RGEllipseGetParameters(ctx->rg,&center,&radius,&vscale));
    382:   PetscCall(RGGetScale(ctx->rg,&rgscale));
    383:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {
    384:     if (!ctx->map) PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)ds),PETSC_DECIDE,ctx->nnod,1,&ctx->map));
    385:     PetscCall(PetscLayoutGetRange(ctx->map,&kstart,&kend));
    386:   }

    388:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W)); /* size n */
    389:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_Q)); /* size mid*n */
    390:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_Z)); /* size mid*n */
    391:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_U)); /* size mid*n */
    392:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_V)); /* size mid*n */
    393:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Q],&Q));
    394:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Z],&Z));
    395:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
    396:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_V],&V));
    397:   mid  = ctx->max_mid;
    398:   PetscCall(PetscBLASIntCast(ds->n,&n));
    399:   p    = n;   /* maximum number of columns for the probing matrix */
    400:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
    401:   PetscCall(PetscBLASIntCast(mid*n,&rowA));
    402:   nw     = 2*n*(p+mid)+3*nnod+2*mid*n*p;
    403:   lrwork = 9*mid*n;

    405:   /* workspace query and memory allocation */
    406:   lwork = -1;
    407:   PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&rowA,Q,&rowA,Z,&rowA,NULL,NULL,NULL,&ld,V,&rowA,&a,&lwork,NULL,&info));
    408:   PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
    409:   PetscCall(DSAllocateWork_Private(ds,lwork+nw,lrwork,n+1));

    411:   sigma = ds->rwork;
    412:   rwork = ds->rwork+mid*n;
    413:   perm  = ds->iwork;
    414:   z     = ds->work+nwu;    nwu += nnod;         /* quadrature points */
    415:   zn    = ds->work+nwu;    nwu += nnod;         /* normalized quadrature points */
    416:   w     = ds->work+nwu;    nwu += nnod;         /* quadrature weights */
    417:   Rc    = ds->work+nwu;    nwu += n*p;
    418:   R     = ds->work+nwu;    nwu += n*p;
    419:   alpha = ds->work+nwu;    nwu += mid*n;
    420:   beta  = ds->work+nwu;    nwu += mid*n;
    421:   S     = ds->work+nwu;    nwu += 2*mid*n*p;
    422:   work  = ds->work+nwu;

    424:   /* Compute quadrature parameters */
    425:   PetscCall(RGComputeQuadrature(ctx->rg,RG_QUADRULE_TRAPEZOIDAL,nnod,z,zn,w));

    427:   /* Set random matrix */
    428:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)ds),&rand));
    429:   PetscCall(PetscRandomSetSeed(rand,0x12345678));
    430:   PetscCall(PetscRandomSeed(rand));
    431:   for (j=0;j<p;j++)
    432:     for (i=0;i<n;i++) PetscCall(PetscRandomGetValue(rand,Rc+i+j*n));
    433:   PetscCall(PetscArrayzero(S,2*mid*n*p));
    434:   /* Loop of integration points */
    435:   for (k=kstart;k<kend;k++) {
    436:     PetscCall(PetscInfo(NULL,"Solving integration point %" PetscInt_FMT "\n",k));
    437:     PetscCall(PetscArraycpy(R,Rc,p*n));
    438:     PetscCall(DSNEPComputeMatrix(ds,z[k],PETSC_FALSE,DS_MAT_W));

    440:     /* LU factorization */
    441:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
    442:     PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n,&n,W,&ld,perm,&info));
    443:     SlepcCheckLapackInfo("getrf",info);
    444:     PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&p,W,&ld,perm,R,&n,&info));
    445:     SlepcCheckLapackInfo("getrs",info);
    446:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));

    448:     /* Moments computation */
    449:     for (s=0;s<2*ctx->max_mid;s++) {
    450:       off = s*n*p;
    451:       for (j=0;j<p;j++)
    452:         for (i=0;i<n;i++) S[off+i+j*n] += w[k]*R[j*n+i];
    453:       w[k] *= zn[k];
    454:     }
    455:   }

    457:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {  /* compute final S via reduction */
    458:     PetscCall(PetscMPIIntCast(2*mid*n*p,&len));
    459:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE,S,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ds)));
    460:   }
    461:   PetscCall(PetscBLASIntCast(ctx->spls?PetscMin(ctx->spls,n):n,&p));
    462:   pp = p;
    463:   do {
    464:     p = pp;
    465:     PetscCall(PetscBLASIntCast(mid*p,&colA));

    467:     PetscCall(PetscInfo(ds,"Computing SVD of size %" PetscBLASInt_FMT "x%" PetscBLASInt_FMT "\n",rowA,colA));
    468:     for (jj=0;jj<mid;jj++) {
    469:       for (ii=0;ii<mid;ii++) {
    470:         off = jj*p*rowA+ii*n;
    471:         for (j=0;j<p;j++)
    472:           for (i=0;i<n;i++) Q[off+j*rowA+i] = S[((jj+ii)*n+j)*n+i];
    473:       }
    474:     }
    475:     PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&rowA,&colA,Q,&rowA,sigma,U,&rowA,V,&colA,work,&lwork,rwork,&info));
    476:     SlepcCheckLapackInfo("gesvd",info);

    478:     rk = colA;
    479:     for (i=1;i<colA;i++) if (sigma[i]/sigma[0]<PETSC_MACHINE_EPSILON*1e4) {rk = i; break;}
    480:     if (rk<colA || p==n) break;
    481:     pp *= 2;
    482:   } while (pp<=n);
    483:   PetscCall(PetscInfo(ds,"Solving generalized eigenproblem of size %" PetscInt_FMT "\n",rk));
    484:   for (jj=0;jj<mid;jj++) {
    485:     for (ii=0;ii<mid;ii++) {
    486:       off = jj*p*rowA+ii*n;
    487:       for (j=0;j<p;j++)
    488:         for (i=0;i<n;i++) Q[off+j*rowA+i] = S[((jj+ii+1)*n+j)*n+i];
    489:     }
    490:   }
    491:   PetscCall(PetscBLASIntCast(rk,&rk_));
    492:   PetscCallBLAS("BLASgemm",BLASgemm_("N","C",&rowA,&rk_,&colA,&sone,Q,&rowA,V,&colA,&szero,Z,&rowA));
    493:   PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&rk_,&rk_,&rowA,&sone,U,&rowA,Z,&rowA,&szero,Q,&rk_));
    494:   PetscCall(PetscArrayzero(Z,n*mid*n*mid));
    495:   for (j=0;j<rk;j++) Z[j+j*rk_] = sigma[j];
    496:   PetscCallBLAS("LAPACK" LAPGEEV,LAPACKggevalt_("N","V",&rk_,Q,&rk_,Z,&rk_,alpha,beta,NULL,&ld,V,&rk_,work,&lwork,rwork,&info));
    497:   SlepcCheckLapackInfo(LAPGEEV,info);
    498:   for (i=0;i<rk;i++) wr[i] = (center+alpha[i]*radius/beta[i])*rgscale;
    499:   PetscCall(PetscMalloc1(rk,&inside));
    500:   PetscCall(RGCheckInside(ctx->rg,rk,wr,wi,inside));
    501:   k=0;
    502:   for (i=0;i<rk;i++)
    503:     if (inside[i]==1) inside[k++] = i;
    504:   /* Discard values outside region */
    505:   lds = ld*mid;
    506:   PetscCall(PetscArrayzero(Q,lds*lds));
    507:   PetscCall(PetscArrayzero(Z,lds*lds));
    508:   for (i=0;i<k;i++) Q[i+i*lds] = (center*beta[inside[i]]+radius*alpha[inside[i]])*rgscale;
    509:   for (i=0;i<k;i++) Z[i+i*lds] = beta[inside[i]];
    510:   for (i=0;i<k;i++) wr[i] = Q[i+i*lds]/Z[i+i*lds];
    511:   for (j=0;j<k;j++) for (i=0;i<rk;i++) V[j*rk+i] = sigma[i]*V[inside[j]*rk+i];
    512:   PetscCall(PetscBLASIntCast(k,&k_));
    513:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
    514:   PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k_,&rk_,&sone,U,&rowA,V,&rk_,&szero,X,&ld));
    515:   /* Normalize */
    516:   for (j=0;j<k;j++) {
    517:     norm = BLASnrm2_(&n,X+ld*j,&one);
    518:     for (i=0;i<n;i++) X[ld*j+i] /= norm;
    519:   }
    520:   PetscCall(PetscFree(inside));
    521:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
    522:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Q],&Q));
    523:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Z],&Z));
    524:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
    525:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_V],&V));

    527:   /* Newton refinement */
    528:   if (ctx->Nit) PetscCall(DSNEPNewtonRefine(ds,k,wr));
    529:   ds->t = k;
    530:   PetscCall(PetscRandomDestroy(&rand));
    531:   PetscFunctionReturn(PETSC_SUCCESS);
    532: }
    533: #endif

    535: #if !defined(PETSC_HAVE_MPIUNI)
    536: static PetscErrorCode DSSynchronize_NEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
    537: {
    538:   DS_NEP         *ctx = (DS_NEP*)ds->data;
    539:   PetscInt       ld=ds->ld,k=0;
    540:   PetscMPIInt    n,n2,rank,size,off=0;
    541:   PetscScalar    *X;

    543:   PetscFunctionBegin;
    544:   if (!ds->method) { /* SLP */
    545:     if (ds->state>=DS_STATE_CONDENSED) k += ds->n;
    546:     if (eigr) k += 1;
    547:     if (eigi) k += 1;
    548:     PetscCall(PetscMPIIntCast(1,&n));
    549:     PetscCall(PetscMPIIntCast(ds->n,&n2));
    550:   } else { /* Contour */
    551:     if (ds->state>=DS_STATE_CONDENSED) k += ctx->max_mid*ds->n*ld;
    552:     if (eigr) k += ctx->max_mid*ds->n;
    553:     if (eigi) k += ctx->max_mid*ds->n;
    554:     PetscCall(PetscMPIIntCast(ctx->max_mid*ds->n,&n));
    555:     PetscCall(PetscMPIIntCast(ctx->max_mid*ds->n*ld,&n2));
    556:   }
    557:   PetscCall(DSAllocateWork_Private(ds,k,0,0));
    558:   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
    559:   if (ds->state>=DS_STATE_CONDENSED) PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
    560:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
    561:   if (!rank) {
    562:     if (ds->state>=DS_STATE_CONDENSED) PetscCallMPI(MPI_Pack(X,n2,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
    563:     if (eigr) PetscCallMPI(MPI_Pack(eigr,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
    564: #if !defined(PETSC_USE_COMPLEX)
    565:     if (eigi) PetscCallMPI(MPI_Pack(eigi,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
    566: #endif
    567:   }
    568:   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
    569:   if (rank) {
    570:     if (ds->state>=DS_STATE_CONDENSED) PetscCallMPI(MPI_Unpack(ds->work,size,&off,X,n2,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
    571:     if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
    572: #if !defined(PETSC_USE_COMPLEX)
    573:     if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
    574: #endif
    575:   }
    576:   if (ds->state>=DS_STATE_CONDENSED) PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
    577:   PetscFunctionReturn(PETSC_SUCCESS);
    578: }
    579: #endif

    581: static PetscErrorCode DSNEPSetFN_NEP(DS ds,PetscInt n,FN fn[])
    582: {
    583:   DS_NEP         *ctx = (DS_NEP*)ds->data;
    584:   PetscInt       i;

    586:   PetscFunctionBegin;
    587:   PetscCheck(n>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more functions, you have %" PetscInt_FMT,n);
    588:   PetscCheck(n<=DS_NUM_EXTRA,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Too many functions, you specified %" PetscInt_FMT " but the limit is %d",n,DS_NUM_EXTRA);
    589:   if (ds->ld) PetscCall(PetscInfo(ds,"DSNEPSetFN() called after DSAllocate()\n"));
    590:   for (i=0;i<n;i++) PetscCall(PetscObjectReference((PetscObject)fn[i]));
    591:   for (i=0;i<ctx->nf;i++) PetscCall(FNDestroy(&ctx->f[i]));
    592:   for (i=0;i<n;i++) ctx->f[i] = fn[i];
    593:   ctx->nf = n;
    594:   PetscFunctionReturn(PETSC_SUCCESS);
    595: }

    597: /*@
    598:    DSNEPSetFN - Sets a number of functions that define the nonlinear
    599:    eigenproblem.

    601:    Collective

    603:    Input Parameters:
    604: +  ds - the direct solver context
    605: .  n  - number of functions
    606: -  fn - array of functions

    608:    Notes:
    609:    The nonlinear eigenproblem is defined in terms of the split nonlinear
    610:    operator $T(\lambda) = \sum_i E_i f_i(\lambda)$.

    612:    This function must be called before `DSAllocate()`. Then `DSAllocate()`
    613:    will allocate an extra matrix $E_i$ per each function, that can be
    614:    filled in the usual way.

    616:    Level: advanced

    618: .seealso: [](sec:ds), `DSNEP`, `DSNEPGetFN()`, `DSAllocate()`
    619: @*/
    620: PetscErrorCode DSNEPSetFN(DS ds,PetscInt n,FN fn[])
    621: {
    622:   PetscInt       i;

    624:   PetscFunctionBegin;
    627:   PetscAssertPointer(fn,3);
    628:   for (i=0;i<n;i++) {
    630:     PetscCheckSameComm(ds,1,fn[i],3);
    631:   }
    632:   PetscTryMethod(ds,"DSNEPSetFN_C",(DS,PetscInt,FN[]),(ds,n,fn));
    633:   PetscFunctionReturn(PETSC_SUCCESS);
    634: }

    636: static PetscErrorCode DSNEPGetFN_NEP(DS ds,PetscInt k,FN *fn)
    637: {
    638:   DS_NEP *ctx = (DS_NEP*)ds->data;

    640:   PetscFunctionBegin;
    641:   PetscCheck(k>=0 && k<ctx->nf,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,ctx->nf-1);
    642:   *fn = ctx->f[k];
    643:   PetscFunctionReturn(PETSC_SUCCESS);
    644: }

    646: /*@
    647:    DSNEPGetFN - Gets the functions associated with the nonlinear `DS`.

    649:    Not Collective

    651:    Input Parameters:
    652: +  ds - the direct solver context
    653: -  k  - the index of the requested function (starting in 0)

    655:    Output Parameter:
    656: .  fn - the function

    658:    Level: advanced

    660: .seealso: [](sec:ds), `DSNEP`, `DSNEPSetFN()`
    661: @*/
    662: PetscErrorCode DSNEPGetFN(DS ds,PetscInt k,FN *fn)
    663: {
    664:   PetscFunctionBegin;
    666:   PetscAssertPointer(fn,3);
    667:   PetscUseMethod(ds,"DSNEPGetFN_C",(DS,PetscInt,FN*),(ds,k,fn));
    668:   PetscFunctionReturn(PETSC_SUCCESS);
    669: }

    671: static PetscErrorCode DSNEPGetNumFN_NEP(DS ds,PetscInt *n)
    672: {
    673:   DS_NEP *ctx = (DS_NEP*)ds->data;

    675:   PetscFunctionBegin;
    676:   *n = ctx->nf;
    677:   PetscFunctionReturn(PETSC_SUCCESS);
    678: }

    680: /*@
    681:    DSNEPGetNumFN - Returns the number of functions stored internally by
    682:    the `DS`.

    684:    Not Collective

    686:    Input Parameter:
    687: .  ds - the direct solver context

    689:    Output Parameter:
    690: .  n - the number of functions passed in `DSNEPSetFN()`

    692:    Level: advanced

    694: .seealso: [](sec:ds), `DSNEP`, `DSNEPSetFN()`
    695: @*/
    696: PetscErrorCode DSNEPGetNumFN(DS ds,PetscInt *n)
    697: {
    698:   PetscFunctionBegin;
    700:   PetscAssertPointer(n,2);
    701:   PetscUseMethod(ds,"DSNEPGetNumFN_C",(DS,PetscInt*),(ds,n));
    702:   PetscFunctionReturn(PETSC_SUCCESS);
    703: }

    705: static PetscErrorCode DSNEPSetMinimality_NEP(DS ds,PetscInt n)
    706: {
    707:   DS_NEP *ctx = (DS_NEP*)ds->data;

    709:   PetscFunctionBegin;
    710:   if (n == PETSC_DECIDE || n == PETSC_DEFAULT) ctx->max_mid = 4;
    711:   else {
    712:     PetscCheck(n>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The minimality value must be > 0");
    713:     ctx->max_mid = n;
    714:   }
    715:   PetscFunctionReturn(PETSC_SUCCESS);
    716: }

    718: /*@
    719:    DSNEPSetMinimality - Sets the maximum minimality index used internally by
    720:    the `DSNEP`.

    722:    Logically Collective

    724:    Input Parameters:
    725: +  ds - the direct solver context
    726: -  n  - the maximum minimality index

    728:    Options Database Key:
    729: .  -ds_nep_minimality n - sets the maximum minimality index

    731:    Notes:
    732:    The maximum minimality index is used only in the contour integral method,
    733:    and is related to the highest moments used in the method. The default
    734:    value is 1, a larger value might give better accuracy in some cases, but
    735:    at a higher cost.

    737:    Level: advanced

    739: .seealso: [](sec:ds), `DSNEP`, `DSNEPGetMinimality()`
    740: @*/
    741: PetscErrorCode DSNEPSetMinimality(DS ds,PetscInt n)
    742: {
    743:   PetscFunctionBegin;
    746:   PetscTryMethod(ds,"DSNEPSetMinimality_C",(DS,PetscInt),(ds,n));
    747:   PetscFunctionReturn(PETSC_SUCCESS);
    748: }

    750: static PetscErrorCode DSNEPGetMinimality_NEP(DS ds,PetscInt *n)
    751: {
    752:   DS_NEP *ctx = (DS_NEP*)ds->data;

    754:   PetscFunctionBegin;
    755:   *n = ctx->max_mid;
    756:   PetscFunctionReturn(PETSC_SUCCESS);
    757: }

    759: /*@
    760:    DSNEPGetMinimality - Returns the maximum minimality index used internally by
    761:    the `DSNEP`.

    763:    Not Collective

    765:    Input Parameter:
    766: .  ds - the direct solver context

    768:    Output Parameter:
    769: .  n - the maximum minimality index passed in `DSNEPSetMinimality()`

    771:    Level: advanced

    773: .seealso: [](sec:ds), `DSNEP`, `DSNEPSetMinimality()`
    774: @*/
    775: PetscErrorCode DSNEPGetMinimality(DS ds,PetscInt *n)
    776: {
    777:   PetscFunctionBegin;
    779:   PetscAssertPointer(n,2);
    780:   PetscUseMethod(ds,"DSNEPGetMinimality_C",(DS,PetscInt*),(ds,n));
    781:   PetscFunctionReturn(PETSC_SUCCESS);
    782: }

    784: static PetscErrorCode DSNEPSetRefine_NEP(DS ds,PetscReal tol,PetscInt its)
    785: {
    786:   DS_NEP *ctx = (DS_NEP*)ds->data;

    788:   PetscFunctionBegin;
    789:   if (tol == (PetscReal)PETSC_DETERMINE) {
    790:     ctx->rtol = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);
    791:   } else if (tol != (PetscReal)PETSC_CURRENT) {
    792:     PetscCheck(tol>0.0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The tolerance must be > 0");
    793:     ctx->rtol = tol;
    794:   }
    795:   if (its == PETSC_DETERMINE) {
    796:     ctx->Nit = 3;
    797:   } else if (its != PETSC_CURRENT) {
    798:     PetscCheck(its>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The number of iterations must be >= 0");
    799:     ctx->Nit = its;
    800:   }
    801:   PetscFunctionReturn(PETSC_SUCCESS);
    802: }

    804: /*@
    805:    DSNEPSetRefine - Sets the tolerance and the number of iterations of Newton iterative
    806:    refinement for eigenpairs when solving a `DSNEP`.

    808:    Logically Collective

    810:    Input Parameters:
    811: +  ds  - the direct solver context
    812: .  tol - the tolerance
    813: -  its - the number of iterations

    815:    Options Database Keys:
    816: +  -ds_nep_refine_tol tol - sets the tolerance
    817: -  -ds_nep_refine_its its - sets the number of Newton iterations

    819:    Notes:
    820:    Iterative refinement of eigenpairs is currently used only in the contour
    821:    integral method.

    823:    Use `PETSC_CURRENT` to retain the current value of any of the parameters.
    824:    Use `PETSC_DETERMINE` for either argument to assign a default value computed
    825:    internally.

    827:    Level: advanced

    829: .seealso: [](sec:ds), `DSNEP`, `DSNEPGetRefine()`
    830: @*/
    831: PetscErrorCode DSNEPSetRefine(DS ds,PetscReal tol,PetscInt its)
    832: {
    833:   PetscFunctionBegin;
    837:   PetscTryMethod(ds,"DSNEPSetRefine_C",(DS,PetscReal,PetscInt),(ds,tol,its));
    838:   PetscFunctionReturn(PETSC_SUCCESS);
    839: }

    841: static PetscErrorCode DSNEPGetRefine_NEP(DS ds,PetscReal *tol,PetscInt *its)
    842: {
    843:   DS_NEP *ctx = (DS_NEP*)ds->data;

    845:   PetscFunctionBegin;
    846:   if (tol) *tol = ctx->rtol;
    847:   if (its) *its = ctx->Nit;
    848:   PetscFunctionReturn(PETSC_SUCCESS);
    849: }

    851: /*@
    852:    DSNEPGetRefine - Returns the tolerance and the number of iterations of Newton iterative
    853:    refinement for eigenpairs.

    855:    Not Collective

    857:    Input Parameter:
    858: .  ds - the direct solver context

    860:    Output Parameters:
    861: +  tol - the tolerance
    862: -  its - the number of iterations

    864:    Level: advanced

    866: .seealso: [](sec:ds), `DSNEP`, `DSNEPSetRefine()`
    867: @*/
    868: PetscErrorCode DSNEPGetRefine(DS ds,PetscReal *tol,PetscInt *its)
    869: {
    870:   PetscFunctionBegin;
    872:   PetscUseMethod(ds,"DSNEPGetRefine_C",(DS,PetscReal*,PetscInt*),(ds,tol,its));
    873:   PetscFunctionReturn(PETSC_SUCCESS);
    874: }

    876: static PetscErrorCode DSNEPSetIntegrationPoints_NEP(DS ds,PetscInt ip)
    877: {
    878:   DS_NEP         *ctx = (DS_NEP*)ds->data;

    880:   PetscFunctionBegin;
    881:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) ctx->nnod = 64;
    882:   else {
    883:     PetscCheck(ip>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The number of integration points must be > 0");
    884:     ctx->nnod = ip;
    885:   }
    886:   PetscCall(PetscLayoutDestroy(&ctx->map));  /* need to redistribute at next solve */
    887:   PetscFunctionReturn(PETSC_SUCCESS);
    888: }

    890: /*@
    891:    DSNEPSetIntegrationPoints - Sets the number of integration points to be
    892:    used in the contour integral method when solving a `DSNEP`.

    894:    Logically Collective

    896:    Input Parameters:
    897: +  ds - the direct solver context
    898: -  ip - the number of integration points

    900:    Options Database Key:
    901: .  -ds_nep_integration_points ip - sets the number of integration points

    903:    Notes:
    904:    This parameter is relevant only in the contour integral method.

    906:    Level: advanced

    908: .seealso: [](sec:ds), `DSNEP`, `DSNEPGetIntegrationPoints()`
    909: @*/
    910: PetscErrorCode DSNEPSetIntegrationPoints(DS ds,PetscInt ip)
    911: {
    912:   PetscFunctionBegin;
    915:   PetscTryMethod(ds,"DSNEPSetIntegrationPoints_C",(DS,PetscInt),(ds,ip));
    916:   PetscFunctionReturn(PETSC_SUCCESS);
    917: }

    919: static PetscErrorCode DSNEPGetIntegrationPoints_NEP(DS ds,PetscInt *ip)
    920: {
    921:   DS_NEP *ctx = (DS_NEP*)ds->data;

    923:   PetscFunctionBegin;
    924:   *ip = ctx->nnod;
    925:   PetscFunctionReturn(PETSC_SUCCESS);
    926: }

    928: /*@
    929:    DSNEPGetIntegrationPoints - Returns the number of integration points used
    930:    in the contour integral method.

    932:    Not Collective

    934:    Input Parameter:
    935: .  ds - the direct solver context

    937:    Output Parameter:
    938: .  ip - the number of integration points

    940:    Level: advanced

    942: .seealso: [](sec:ds), `DSNEP`, `DSNEPSetIntegrationPoints()`
    943: @*/
    944: PetscErrorCode DSNEPGetIntegrationPoints(DS ds,PetscInt *ip)
    945: {
    946:   PetscFunctionBegin;
    948:   PetscAssertPointer(ip,2);
    949:   PetscUseMethod(ds,"DSNEPGetIntegrationPoints_C",(DS,PetscInt*),(ds,ip));
    950:   PetscFunctionReturn(PETSC_SUCCESS);
    951: }

    953: static PetscErrorCode DSNEPSetSamplingSize_NEP(DS ds,PetscInt p)
    954: {
    955:   DS_NEP *ctx = (DS_NEP*)ds->data;

    957:   PetscFunctionBegin;
    958:   if (p == PETSC_DECIDE || p == PETSC_DEFAULT) ctx->spls = 0;
    959:   else {
    960:     PetscCheck(p>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The sample size must be > 0");
    961:     PetscCheck(p>=20,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The sample size cannot be smaller than 20");
    962:     ctx->spls = p;
    963:   }
    964:   PetscFunctionReturn(PETSC_SUCCESS);
    965: }

    967: /*@
    968:    DSNEPSetSamplingSize - Sets the number of sampling columns to be
    969:    used in the contour integral method when solving a `DSNEP`.

    971:    Logically Collective

    973:    Input Parameters:
    974: +  ds - the direct solver context
    975: -  p  - the number of columns for the sampling matrix

    977:    Options Database Key:
    978: .  -ds_nep_sampling_size p - set the number of sampling columns

    980:    Note:
    981:    This parameter is relevant only in the contour integral method.

    983:    Level: advanced

    985: .seealso: [](sec:ds), `DSNEP`, `DSNEPGetSamplingSize()`
    986: @*/
    987: PetscErrorCode DSNEPSetSamplingSize(DS ds,PetscInt p)
    988: {
    989:   PetscFunctionBegin;
    992:   PetscTryMethod(ds,"DSNEPSetSamplingSize_C",(DS,PetscInt),(ds,p));
    993:   PetscFunctionReturn(PETSC_SUCCESS);
    994: }

    996: static PetscErrorCode DSNEPGetSamplingSize_NEP(DS ds,PetscInt *p)
    997: {
    998:   DS_NEP *ctx = (DS_NEP*)ds->data;

    1000:   PetscFunctionBegin;
    1001:   *p = ctx->spls;
    1002:   PetscFunctionReturn(PETSC_SUCCESS);
    1003: }

    1005: /*@
    1006:    DSNEPGetSamplingSize - Returns the number of sampling columns used
    1007:    in the contour integral method.

    1009:    Not Collective

    1011:    Input Parameter:
    1012: .  ds - the direct solver context

    1014:    Output Parameter:
    1015: .  p -  the number of columns for the sampling matrix

    1017:    Level: advanced

    1019: .seealso: [](sec:ds), `DSNEP`, `DSNEPSetSamplingSize()`
    1020: @*/
    1021: PetscErrorCode DSNEPGetSamplingSize(DS ds,PetscInt *p)
    1022: {
    1023:   PetscFunctionBegin;
    1025:   PetscAssertPointer(p,2);
    1026:   PetscUseMethod(ds,"DSNEPGetSamplingSize_C",(DS,PetscInt*),(ds,p));
    1027:   PetscFunctionReturn(PETSC_SUCCESS);
    1028: }

    1030: static PetscErrorCode DSNEPSetComputeMatrixFunction_NEP(DS ds,DSNEPMatrixFunctionFn *fun,PetscCtx ctx,PetscCtxDestroyFn *destroy)
    1031: {
    1032:   DS_NEP *dsctx = (DS_NEP*)ds->data;

    1034:   PetscFunctionBegin;
    1035:   if (dsctx->computematrixdestroy) PetscCall((*dsctx->computematrixdestroy)(&dsctx->computematrixctx));
    1036:   dsctx->computematrix        = fun;
    1037:   dsctx->computematrixctx     = ctx;
    1038:   dsctx->computematrixdestroy = destroy;
    1039:   PetscFunctionReturn(PETSC_SUCCESS);
    1040: }

    1042: /*@C
    1043:    DSNEPSetComputeMatrixFunction - Sets a user-provided subroutine to compute
    1044:    the matrices $T(\lambda)$ or $T'(\lambda)$.

    1046:    Logically Collective

    1048:    Input Parameters:
    1049: +  ds      - the direct solver context
    1050: .  fun     - matrix function evaluation routine, see `DSNEPMatrixFunctionFn` for the calling sequence
    1051: .  ctx     - a context pointer (the last parameter to the user function)
    1052: -  destroy - a routine for destroying the context (may be `NULL`), see `PetscCtxDestroyFn`
    1053:              for the calling sequence

    1055:    Note:
    1056:    The result is computed as $T(\lambda) = \sum_i E_i f_i(\lambda)$, and similarly
    1057:    for the derivative, where $E_i$ are the extra matrices, see `DSMatType`.

    1059:    Level: developer

    1061: .seealso: [](sec:ds), `DSNEP`, `DSNEPGetComputeMatrixFunction()`
    1062: @*/
    1063: PetscErrorCode DSNEPSetComputeMatrixFunction(DS ds,DSNEPMatrixFunctionFn *fun,PetscCtx ctx,PetscCtxDestroyFn *destroy)
    1064: {
    1065:   PetscFunctionBegin;
    1067:   PetscTryMethod(ds,"DSNEPSetComputeMatrixFunction_C",(DS,DSNEPMatrixFunctionFn*,PetscCtx,PetscCtxDestroyFn*),(ds,fun,ctx,destroy));
    1068:   PetscFunctionReturn(PETSC_SUCCESS);
    1069: }

    1071: static PetscErrorCode DSNEPGetComputeMatrixFunction_NEP(DS ds,DSNEPMatrixFunctionFn **fun,PetscCtxRt ctx,PetscCtxDestroyFn **destroy)
    1072: {
    1073:   DS_NEP *dsctx = (DS_NEP*)ds->data;

    1075:   PetscFunctionBegin;
    1076:   if (fun) *fun = dsctx->computematrix;
    1077:   if (ctx) *(void**)ctx = dsctx->computematrixctx;
    1078:   if (destroy) *destroy = dsctx->computematrixdestroy;
    1079:   PetscFunctionReturn(PETSC_SUCCESS);
    1080: }

    1082: /*@C
    1083:    DSNEPGetComputeMatrixFunction - Returns the user-provided callback function
    1084:    set in `DSNEPSetComputeMatrixFunction()`.

    1086:    Not Collective

    1088:    Input Parameter:
    1089: .  ds  - the direct solver context

    1091:    Output Parameters:
    1092: +  fun     - the pointer to the user function
    1093: .  ctx     - the context pointer
    1094: -  destroy - a routine for destroying the context (may be `NULL`)

    1096:    Level: developer

    1098: .seealso: [](sec:ds), `DSNEP`, `DSNEPSetComputeMatrixFunction()`
    1099: @*/
    1100: PetscErrorCode DSNEPGetComputeMatrixFunction(DS ds,DSNEPMatrixFunctionFn **fun,PetscCtxRt ctx,PetscCtxDestroyFn **destroy)
    1101: {
    1102:   PetscFunctionBegin;
    1104:   PetscUseMethod(ds,"DSNEPGetComputeMatrixFunction_C",(DS,DSNEPMatrixFunctionFn**,PetscCtxRt,PetscCtxDestroyFn**),(ds,fun,ctx,destroy));
    1105:   PetscFunctionReturn(PETSC_SUCCESS);
    1106: }

    1108: static PetscErrorCode DSNEPSetRG_NEP(DS ds,RG rg)
    1109: {
    1110:   DS_NEP         *dsctx = (DS_NEP*)ds->data;

    1112:   PetscFunctionBegin;
    1113:   PetscCall(PetscObjectReference((PetscObject)rg));
    1114:   PetscCall(RGDestroy(&dsctx->rg));
    1115:   dsctx->rg = rg;
    1116:   PetscFunctionReturn(PETSC_SUCCESS);
    1117: }

    1119: /*@
    1120:    DSNEPSetRG - Associates a region object to the `DSNEP` solver.

    1122:    Collective

    1124:    Input Parameters:
    1125: +  ds  - the direct solver context
    1126: -  rg  - the region context

    1128:    Notes:
    1129:    The region is used only in the contour integral method, and
    1130:    should enclose the wanted eigenvalues.

    1132:    Level: developer

    1134: .seealso: [](sec:ds), `DSNEP`, `DSNEPGetRG()`
    1135: @*/
    1136: PetscErrorCode DSNEPSetRG(DS ds,RG rg)
    1137: {
    1138:   PetscFunctionBegin;
    1140:   if (rg) {
    1142:     PetscCheckSameComm(ds,1,rg,2);
    1143:   }
    1144:   PetscTryMethod(ds,"DSNEPSetRG_C",(DS,RG),(ds,rg));
    1145:   PetscFunctionReturn(PETSC_SUCCESS);
    1146: }

    1148: static PetscErrorCode DSNEPGetRG_NEP(DS ds,RG *rg)
    1149: {
    1150:   DS_NEP         *ctx = (DS_NEP*)ds->data;

    1152:   PetscFunctionBegin;
    1153:   if (!ctx->rg) {
    1154:     PetscCall(RGCreate(PetscObjectComm((PetscObject)ds),&ctx->rg));
    1155:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->rg,(PetscObject)ds,1));
    1156:     PetscCall(RGSetOptionsPrefix(ctx->rg,((PetscObject)ds)->prefix));
    1157:     PetscCall(RGAppendOptionsPrefix(ctx->rg,"ds_nep_"));
    1158:     PetscCall(PetscObjectSetOptions((PetscObject)ctx->rg,((PetscObject)ds)->options));
    1159:   }
    1160:   *rg = ctx->rg;
    1161:   PetscFunctionReturn(PETSC_SUCCESS);
    1162: }

    1164: /*@
    1165:    DSNEPGetRG - Obtain the region object associated to the `DSNEP` solver.

    1167:    Collective

    1169:    Input Parameter:
    1170: .  ds  - the direct solver context

    1172:    Output Parameter:
    1173: .  rg  - the region context

    1175:    Level: developer

    1177: .seealso: [](sec:ds), `DSNEP`, `DSNEPSetRG()`
    1178: @*/
    1179: PetscErrorCode DSNEPGetRG(DS ds,RG *rg)
    1180: {
    1181:   PetscFunctionBegin;
    1183:   PetscAssertPointer(rg,2);
    1184:   PetscUseMethod(ds,"DSNEPGetRG_C",(DS,RG*),(ds,rg));
    1185:   PetscFunctionReturn(PETSC_SUCCESS);
    1186: }

    1188: static PetscErrorCode DSSetFromOptions_NEP(DS ds,PetscOptionItems PetscOptionsObject)
    1189: {
    1190:   PetscInt       k;
    1191:   PetscBool      flg;
    1192: #if defined(PETSC_USE_COMPLEX)
    1193:   PetscReal      r;
    1194:   PetscBool      flg1;
    1195:   DS_NEP         *ctx = (DS_NEP*)ds->data;
    1196: #endif

    1198:   PetscFunctionBegin;
    1199:   PetscOptionsHeadBegin(PetscOptionsObject,"DS NEP Options");

    1201:     PetscCall(PetscOptionsInt("-ds_nep_minimality","Maximum minimality index","DSNEPSetMinimality",4,&k,&flg));
    1202:     if (flg) PetscCall(DSNEPSetMinimality(ds,k));

    1204:     PetscCall(PetscOptionsInt("-ds_nep_integration_points","Number of integration points","DSNEPSetIntegrationPoints",64,&k,&flg));
    1205:     if (flg) PetscCall(DSNEPSetIntegrationPoints(ds,k));

    1207:     PetscCall(PetscOptionsInt("-ds_nep_sampling_size","Number of sampling columns","DSNEPSetSamplingSize",0,&k,&flg));
    1208:     if (flg) PetscCall(DSNEPSetSamplingSize(ds,k));

    1210: #if defined(PETSC_USE_COMPLEX)
    1211:     r = ctx->rtol;
    1212:     PetscCall(PetscOptionsReal("-ds_nep_refine_tol","Refinement tolerance","DSNEPSetRefine",ctx->rtol,&r,&flg1));
    1213:     k = ctx->Nit;
    1214:     PetscCall(PetscOptionsInt("-ds_nep_refine_its","Number of iterative refinement iterations","DSNEPSetRefine",ctx->Nit,&k,&flg));
    1215:     if (flg1||flg) PetscCall(DSNEPSetRefine(ds,r,k));

    1217:     if (ds->method==1) {
    1218:       if (!ctx->rg) PetscCall(DSNEPGetRG(ds,&ctx->rg));
    1219:       PetscCall(RGSetFromOptions(ctx->rg));
    1220:     }
    1221: #endif

    1223:   PetscOptionsHeadEnd();
    1224:   PetscFunctionReturn(PETSC_SUCCESS);
    1225: }

    1227: static PetscErrorCode DSDestroy_NEP(DS ds)
    1228: {
    1229:   DS_NEP         *ctx = (DS_NEP*)ds->data;
    1230:   PetscInt       i;

    1232:   PetscFunctionBegin;
    1233:   for (i=0;i<ctx->nf;i++) PetscCall(FNDestroy(&ctx->f[i]));
    1234:   PetscCall(RGDestroy(&ctx->rg));
    1235:   PetscCall(PetscLayoutDestroy(&ctx->map));
    1236:   if (ctx->computematrixdestroy) PetscCall((*ctx->computematrixdestroy)(&ctx->computematrixctx));
    1237:   PetscCall(PetscFree(ds->data));
    1238:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",NULL));
    1239:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",NULL));
    1240:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",NULL));
    1241:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",NULL));
    1242:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",NULL));
    1243:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",NULL));
    1244:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",NULL));
    1245:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",NULL));
    1246:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",NULL));
    1247:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",NULL));
    1248:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",NULL));
    1249:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",NULL));
    1250:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",NULL));
    1251:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",NULL));
    1252:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",NULL));
    1253:   PetscFunctionReturn(PETSC_SUCCESS);
    1254: }

    1256: static PetscErrorCode DSMatGetSize_NEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
    1257: {
    1258:   DS_NEP *ctx = (DS_NEP*)ds->data;

    1260:   PetscFunctionBegin;
    1261:   *rows = ds->n;
    1262:   if (t==DS_MAT_Q || t==DS_MAT_Z || t==DS_MAT_U || t==DS_MAT_V) *rows *= ctx->max_mid;
    1263:   *cols = ds->n;
    1264:   if (t==DS_MAT_Q || t==DS_MAT_Z || t==DS_MAT_U || t==DS_MAT_V || t==DS_MAT_X || t==DS_MAT_Y) *cols *= ctx->max_mid;
    1265:   PetscFunctionReturn(PETSC_SUCCESS);
    1266: }

    1268: /*MC
    1269:    DSNEP - Dense Nonlinear Eigenvalue Problem.

    1271:    Notes:
    1272:    The problem is expressed as $T(\lambda)x = 0$, where $T(\lambda)$ is a
    1273:    parameter-dependent matrix written as $T(\lambda) = \sum_i E_i f_i(\lambda)$.
    1274:    The eigenvalues $\lambda$ are the arguments returned by `DSSolve()`.

    1276:    The coefficient matrices $E_i$ are the extra matrices of the `DS`, and
    1277:    the scalar functions $f_i$ are passed via `DSNEPSetFN()`. Optionally, a
    1278:    callback function to fill the $E_i$ matrices can be set with
    1279:    `DSNEPSetComputeMatrixFunction()`.

    1281:    Used DS matrices:
    1282: +  `DS_MAT_E0` to `DS_MAT_E9` - coefficient matrices of the split form of $T(\lambda)$
    1283: .  `DS_MAT_X`  - eigenvectors
    1284: .  `DS_MAT_A`  - (workspace) $T(\lambda)$ evaluated at a given $\lambda$ (SLP only)
    1285: .  `DS_MAT_B`  - (workspace) $T'(\lambda)$ evaluated at a given $\lambda$ (SLP only)
    1286: .  `DS_MAT_Q`  - (workspace) left Hankel matrix (contour only)
    1287: .  `DS_MAT_Z`  - (workspace) right Hankel matrix (contour only)
    1288: .  `DS_MAT_U`  - (workspace) left singular vectors (contour only)
    1289: .  `DS_MAT_V`  - (workspace) right singular vectors (contour only)
    1290: -  `DS_MAT_W`  - (workspace) auxiliary matrix of size $n\times n$

    1292:    Implemented methods:
    1293: +  0 - Successive Linear Problems (SLP), computes just one eigenpair
    1294: -  1 - Contour integral, computes all eigenvalues inside a region

    1296:    Level: beginner

    1298: .seealso: [](sec:ds), `DSCreate()`, `DSSetType()`, `DSType`, `DSNEPSetFN()`, `DSNEPSetComputeMatrixFunction()`
    1299: M*/
    1300: SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS ds)
    1301: {
    1302:   DS_NEP         *ctx;

    1304:   PetscFunctionBegin;
    1305:   PetscCall(PetscNew(&ctx));
    1306:   ds->data = (void*)ctx;
    1307:   ctx->max_mid = 4;
    1308:   ctx->nnod    = 64;
    1309:   ctx->Nit     = 3;
    1310:   ctx->rtol    = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);

    1312:   ds->ops->allocate       = DSAllocate_NEP;
    1313:   ds->ops->setfromoptions = DSSetFromOptions_NEP;
    1314:   ds->ops->view           = DSView_NEP;
    1315:   ds->ops->vectors        = DSVectors_NEP;
    1316:   ds->ops->solve[0]       = DSSolve_NEP_SLP;
    1317: #if defined(PETSC_USE_COMPLEX)
    1318:   ds->ops->solve[1]       = DSSolve_NEP_Contour;
    1319: #endif
    1320:   ds->ops->sort           = DSSort_NEP;
    1321: #if !defined(PETSC_HAVE_MPIUNI)
    1322:   ds->ops->synchronize    = DSSynchronize_NEP;
    1323: #endif
    1324:   ds->ops->destroy        = DSDestroy_NEP;
    1325:   ds->ops->matgetsize     = DSMatGetSize_NEP;

    1327:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",DSNEPSetFN_NEP));
    1328:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",DSNEPGetFN_NEP));
    1329:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",DSNEPGetNumFN_NEP));
    1330:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",DSNEPGetMinimality_NEP));
    1331:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",DSNEPSetMinimality_NEP));
    1332:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",DSNEPGetRefine_NEP));
    1333:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",DSNEPSetRefine_NEP));
    1334:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",DSNEPGetIntegrationPoints_NEP));
    1335:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",DSNEPSetIntegrationPoints_NEP));
    1336:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",DSNEPGetSamplingSize_NEP));
    1337:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",DSNEPSetSamplingSize_NEP));
    1338:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",DSNEPSetRG_NEP));
    1339:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",DSNEPGetRG_NEP));
    1340:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",DSNEPSetComputeMatrixFunction_NEP));
    1341:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",DSNEPGetComputeMatrixFunction_NEP));
    1342:   PetscFunctionReturn(PETSC_SUCCESS);
    1343: }