Actual source code: dsnep.c

slepc-3.19.2 2023-09-05
Report Typos and Errors
  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:   void           *computematrixctx;
 25:   PetscErrorCode (*computematrix)(DS,PetscScalar,PetscBool,DSMatType,void*);
 26: } DS_NEP;

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

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

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

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

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

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

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

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

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

151: PetscErrorCode DSSolve_NEP_SLP(DS ds,PetscScalar *wr,PetscScalar *wi)
152: {
153:   PetscScalar    *A,*B,*W,*X,*work,*alpha,*beta;
154:   PetscScalar    sigma,lambda,mu,re,re2,sone=1.0,szero=0.0;
155:   PetscBLASInt   info,n,ld,lrwork=0,lwork,one=1,zero=0;
156:   PetscInt       it,pos,j,maxit=100,result;
157:   PetscReal      norm,tol,done=1.0;
158: #if defined(PETSC_USE_COMPLEX)
159:   PetscReal      *rwork;
160: #else
161:   PetscReal      *alphai,im,im2;
162: #endif

164:   PetscFunctionBegin;
165:   PetscCall(PetscBLASIntCast(ds->n,&n));
166:   PetscCall(PetscBLASIntCast(ds->ld,&ld));
167: #if defined(PETSC_USE_COMPLEX)
168:   PetscCall(PetscBLASIntCast(2*ds->n+2*ds->n,&lwork));
169:   PetscCall(PetscBLASIntCast(8*ds->n,&lrwork));
170: #else
171:   PetscCall(PetscBLASIntCast(3*ds->n+8*ds->n,&lwork));
172: #endif
173:   PetscCall(DSAllocateWork_Private(ds,lwork,lrwork,0));
174:   alpha = ds->work;
175:   beta = ds->work + ds->n;
176: #if defined(PETSC_USE_COMPLEX)
177:   work = ds->work + 2*ds->n;
178:   lwork -= 2*ds->n;
179: #else
180:   alphai = ds->work + 2*ds->n;
181:   work = ds->work + 3*ds->n;
182:   lwork -= 3*ds->n;
183: #endif
184:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
185:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
186:   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
187:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
188:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
189:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
190:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));

192:   sigma = 0.0;
193:   if (ds->sc->comparison==SlepcCompareTargetMagnitude || ds->sc->comparison==SlepcCompareTargetReal) sigma = *(PetscScalar*)ds->sc->comparisonctx;
194:   lambda = sigma;
195:   tol = n*PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);

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

199:     /* evaluate T and T' */
200:     PetscCall(DSNEPComputeMatrix(ds,lambda,PETSC_FALSE,DS_MAT_A));
201:     if (it) {
202:       PetscCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,A,&ld,X,&one,&szero,X+ld,&one));
203:       norm = BLASnrm2_(&n,X+ld,&one);
204:       if (norm/PetscAbsScalar(lambda)<=tol) break;
205:     }
206:     PetscCall(DSNEPComputeMatrix(ds,lambda,PETSC_TRUE,DS_MAT_B));

208:     /* compute eigenvalue correction mu and eigenvector u */
209: #if defined(PETSC_USE_COMPLEX)
210:     rwork = ds->rwork;
211:     PetscCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,beta,NULL,&ld,W,&ld,work,&lwork,rwork,&info));
212: #else
213:     PetscCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,alphai,beta,NULL,&ld,W,&ld,work,&lwork,&info));
214: #endif
215:     SlepcCheckLapackInfo("ggev",info);

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

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

254:     /* correct eigenvalue approximation */
255:     lambda = lambda - mu;
256:   }
257:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
258:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
259:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
260:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));

262:   PetscCheck(it<maxit,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"DSNEP did not converge");
263:   ds->t = 1;
264:   wr[0] = lambda;
265:   if (wi) wi[0] = 0.0;
266:   PetscFunctionReturn(PETSC_SUCCESS);
267: }

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

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

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

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

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

403:   sigma = ds->rwork;
404:   rwork = ds->rwork+mid*n;
405:   perm  = ds->iwork;
406:   z     = ds->work+nwu;    nwu += nnod;         /* quadrature points */
407:   zn    = ds->work+nwu;    nwu += nnod;         /* normalized quadrature points */
408:   w     = ds->work+nwu;    nwu += nnod;         /* quadrature weights */
409:   Rc    = ds->work+nwu;    nwu += n*p;
410:   R     = ds->work+nwu;    nwu += n*p;
411:   alpha = ds->work+nwu;    nwu += mid*n;
412:   beta  = ds->work+nwu;    nwu += mid*n;
413:   S     = ds->work+nwu;    nwu += 2*mid*n*p;
414:   work  = ds->work+nwu;    /*nwu += mid*n*5;*/

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

419:   /* Set random matrix */
420:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)ds),&rand));
421:   PetscCall(PetscRandomSetSeed(rand,0x12345678));
422:   PetscCall(PetscRandomSeed(rand));
423:   for (j=0;j<p;j++)
424:     for (i=0;i<n;i++) PetscCall(PetscRandomGetValue(rand,Rc+i+j*n));
425:   PetscCall(PetscArrayzero(S,2*mid*n*p));
426:   /* Loop of integration points */
427:   for (k=kstart;k<kend;k++) {
428:     PetscCall(PetscInfo(NULL,"Solving integration point %" PetscInt_FMT "\n",k));
429:     PetscCall(PetscArraycpy(R,Rc,p*n));
430:     PetscCall(DSNEPComputeMatrix(ds,z[k],PETSC_FALSE,DS_MAT_W));

432:     /* LU factorization */
433:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
434:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
435:     PetscCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n,&n,W,&ld,perm,&info));
436:     SlepcCheckLapackInfo("getrf",info);
437:     PetscCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&p,W,&ld,perm,R,&n,&info));
438:     SlepcCheckLapackInfo("getrs",info);
439:     PetscCall(PetscFPTrapPop());
440:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));

442:     /* Moments computation */
443:     for (s=0;s<2*ctx->max_mid;s++) {
444:       off = s*n*p;
445:       for (j=0;j<p;j++)
446:         for (i=0;i<n;i++) S[off+i+j*n] += w[k]*R[j*n+i];
447:       w[k] *= zn[k];
448:     }
449:   }

451:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {  /* compute final S via reduction */
452:     PetscCall(PetscMPIIntCast(2*mid*n*p,&len));
453:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE,S,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ds)));
454:   }
455:   p = ctx->spls?PetscMin(ctx->spls,n):n;
456:   pp = p;
457:   do {
458:     p = pp;
459:     PetscCall(PetscBLASIntCast(mid*p,&colA));

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

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

522:   /* Newton refinement */
523:   if (ctx->Nit) PetscCall(DSNEPNewtonRefine(ds,k,wr));
524:   ds->t = k;
525:   PetscCall(PetscRandomDestroy(&rand));
526:   PetscFunctionReturn(PETSC_SUCCESS);
527: }
528: #endif

530: #if !defined(PETSC_HAVE_MPIUNI)
531: PetscErrorCode DSSynchronize_NEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
532: {
533:   DS_NEP         *ctx = (DS_NEP*)ds->data;
534:   PetscInt       ld=ds->ld,k=0;
535:   PetscMPIInt    n,n2,rank,size,off=0;
536:   PetscScalar    *X;

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

576: static PetscErrorCode DSNEPSetFN_NEP(DS ds,PetscInt n,FN fn[])
577: {
578:   DS_NEP         *ctx = (DS_NEP*)ds->data;
579:   PetscInt       i;

581:   PetscFunctionBegin;
582:   PetscCheck(n>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more functions, you have %" PetscInt_FMT,n);
583:   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);
584:   if (ds->ld) PetscCall(PetscInfo(ds,"DSNEPSetFN() called after DSAllocate()\n"));
585:   for (i=0;i<n;i++) PetscCall(PetscObjectReference((PetscObject)fn[i]));
586:   for (i=0;i<ctx->nf;i++) PetscCall(FNDestroy(&ctx->f[i]));
587:   for (i=0;i<n;i++) ctx->f[i] = fn[i];
588:   ctx->nf = n;
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

592: /*@
593:    DSNEPSetFN - Sets a number of functions that define the nonlinear
594:    eigenproblem.

596:    Collective

598:    Input Parameters:
599: +  ds - the direct solver context
600: .  n  - number of functions
601: -  fn - array of functions

603:    Notes:
604:    The nonlinear eigenproblem is defined in terms of the split nonlinear
605:    operator T(lambda) = sum_i A_i*f_i(lambda).

607:    This function must be called before DSAllocate(). Then DSAllocate()
608:    will allocate an extra matrix A_i per each function, that can be
609:    filled in the usual way.

611:    Level: advanced

613: .seealso: DSNEPGetFN(), DSAllocate()
614: @*/
615: PetscErrorCode DSNEPSetFN(DS ds,PetscInt n,FN fn[])
616: {
617:   PetscInt       i;

619:   PetscFunctionBegin;
623:   for (i=0;i<n;i++) {
625:     PetscCheckSameComm(ds,1,fn[i],3);
626:   }
627:   PetscTryMethod(ds,"DSNEPSetFN_C",(DS,PetscInt,FN[]),(ds,n,fn));
628:   PetscFunctionReturn(PETSC_SUCCESS);
629: }

631: static PetscErrorCode DSNEPGetFN_NEP(DS ds,PetscInt k,FN *fn)
632: {
633:   DS_NEP *ctx = (DS_NEP*)ds->data;

635:   PetscFunctionBegin;
636:   PetscCheck(k>=0 && k<ctx->nf,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,ctx->nf-1);
637:   *fn = ctx->f[k];
638:   PetscFunctionReturn(PETSC_SUCCESS);
639: }

641: /*@
642:    DSNEPGetFN - Gets the functions associated with the nonlinear DS.

644:    Not Collective

646:    Input Parameters:
647: +  ds - the direct solver context
648: -  k  - the index of the requested function (starting in 0)

650:    Output Parameter:
651: .  fn - the function

653:    Level: advanced

655: .seealso: DSNEPSetFN()
656: @*/
657: PetscErrorCode DSNEPGetFN(DS ds,PetscInt k,FN *fn)
658: {
659:   PetscFunctionBegin;
662:   PetscUseMethod(ds,"DSNEPGetFN_C",(DS,PetscInt,FN*),(ds,k,fn));
663:   PetscFunctionReturn(PETSC_SUCCESS);
664: }

666: static PetscErrorCode DSNEPGetNumFN_NEP(DS ds,PetscInt *n)
667: {
668:   DS_NEP *ctx = (DS_NEP*)ds->data;

670:   PetscFunctionBegin;
671:   *n = ctx->nf;
672:   PetscFunctionReturn(PETSC_SUCCESS);
673: }

675: /*@
676:    DSNEPGetNumFN - Returns the number of functions stored internally by
677:    the DS.

679:    Not Collective

681:    Input Parameter:
682: .  ds - the direct solver context

684:    Output Parameters:
685: .  n - the number of functions passed in DSNEPSetFN()

687:    Level: advanced

689: .seealso: DSNEPSetFN()
690: @*/
691: PetscErrorCode DSNEPGetNumFN(DS ds,PetscInt *n)
692: {
693:   PetscFunctionBegin;
696:   PetscUseMethod(ds,"DSNEPGetNumFN_C",(DS,PetscInt*),(ds,n));
697:   PetscFunctionReturn(PETSC_SUCCESS);
698: }

700: static PetscErrorCode DSNEPSetMinimality_NEP(DS ds,PetscInt n)
701: {
702:   DS_NEP *ctx = (DS_NEP*)ds->data;

704:   PetscFunctionBegin;
705:   if (n == PETSC_DECIDE || n == PETSC_DEFAULT) ctx->max_mid = 4;
706:   else {
707:     PetscCheck(n>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The minimality value must be > 0");
708:     ctx->max_mid = n;
709:   }
710:   PetscFunctionReturn(PETSC_SUCCESS);
711: }

713: /*@
714:    DSNEPSetMinimality - Sets the maximum minimality index used internally by
715:    the DSNEP.

717:    Logically Collective

719:    Input Parameters:
720: +  ds - the direct solver context
721: -  n  - the maximum minimality index

723:    Options Database Key:
724: .  -ds_nep_minimality <n> - sets the maximum minimality index

726:    Notes:
727:    The maximum minimality index is used only in the contour integral method,
728:    and is related to the highest momemts used in the method. The default
729:    value is 1, an larger value might give better accuracy in some cases, but
730:    at a higher cost.

732:    Level: advanced

734: .seealso: DSNEPGetMinimality()
735: @*/
736: PetscErrorCode DSNEPSetMinimality(DS ds,PetscInt n)
737: {
738:   PetscFunctionBegin;
741:   PetscTryMethod(ds,"DSNEPSetMinimality_C",(DS,PetscInt),(ds,n));
742:   PetscFunctionReturn(PETSC_SUCCESS);
743: }

745: static PetscErrorCode DSNEPGetMinimality_NEP(DS ds,PetscInt *n)
746: {
747:   DS_NEP *ctx = (DS_NEP*)ds->data;

749:   PetscFunctionBegin;
750:   *n = ctx->max_mid;
751:   PetscFunctionReturn(PETSC_SUCCESS);
752: }

754: /*@
755:    DSNEPGetMinimality - Returns the maximum minimality index used internally by
756:    the DSNEP.

758:    Not Collective

760:    Input Parameter:
761: .  ds - the direct solver context

763:    Output Parameters:
764: .  n - the maximum minimality index passed in DSNEPSetMinimality()

766:    Level: advanced

768: .seealso: DSNEPSetMinimality()
769: @*/
770: PetscErrorCode DSNEPGetMinimality(DS ds,PetscInt *n)
771: {
772:   PetscFunctionBegin;
775:   PetscUseMethod(ds,"DSNEPGetMinimality_C",(DS,PetscInt*),(ds,n));
776:   PetscFunctionReturn(PETSC_SUCCESS);
777: }

779: static PetscErrorCode DSNEPSetRefine_NEP(DS ds,PetscReal tol,PetscInt its)
780: {
781:   DS_NEP *ctx = (DS_NEP*)ds->data;

783:   PetscFunctionBegin;
784:   if (tol == (PetscReal)PETSC_DEFAULT) ctx->rtol = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);
785:   else {
786:     PetscCheck(tol>0.0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The tolerance must be > 0");
787:     ctx->rtol = tol;
788:   }
789:   if (its == PETSC_DECIDE || its == PETSC_DEFAULT) ctx->Nit = 3;
790:   else {
791:     PetscCheck(its>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The number of iterations must be >= 0");
792:     ctx->Nit = its;
793:   }
794:   PetscFunctionReturn(PETSC_SUCCESS);
795: }

797: /*@
798:    DSNEPSetRefine - Sets the tolerance and the number of iterations of Newton iterative
799:    refinement for eigenpairs.

801:    Logically Collective

803:    Input Parameters:
804: +  ds  - the direct solver context
805: .  tol - the tolerance
806: -  its - the number of iterations

808:    Options Database Key:
809: +  -ds_nep_refine_tol <tol> - sets the tolerance
810: -  -ds_nep_refine_its <its> - sets the number of Newton iterations

812:    Notes:
813:    Iterative refinement of eigenpairs is currently used only in the contour
814:    integral method.

816:    Level: advanced

818: .seealso: DSNEPGetRefine()
819: @*/
820: PetscErrorCode DSNEPSetRefine(DS ds,PetscReal tol,PetscInt its)
821: {
822:   PetscFunctionBegin;
826:   PetscTryMethod(ds,"DSNEPSetRefine_C",(DS,PetscReal,PetscInt),(ds,tol,its));
827:   PetscFunctionReturn(PETSC_SUCCESS);
828: }

830: static PetscErrorCode DSNEPGetRefine_NEP(DS ds,PetscReal *tol,PetscInt *its)
831: {
832:   DS_NEP *ctx = (DS_NEP*)ds->data;

834:   PetscFunctionBegin;
835:   if (tol) *tol = ctx->rtol;
836:   if (its) *its = ctx->Nit;
837:   PetscFunctionReturn(PETSC_SUCCESS);
838: }

840: /*@
841:    DSNEPGetRefine - Returns the tolerance and the number of iterations of Newton iterative
842:    refinement for eigenpairs.

844:    Not Collective

846:    Input Parameter:
847: .  ds - the direct solver context

849:    Output Parameters:
850: +  tol - the tolerance
851: -  its - the number of iterations

853:    Level: advanced

855: .seealso: DSNEPSetRefine()
856: @*/
857: PetscErrorCode DSNEPGetRefine(DS ds,PetscReal *tol,PetscInt *its)
858: {
859:   PetscFunctionBegin;
861:   PetscUseMethod(ds,"DSNEPGetRefine_C",(DS,PetscReal*,PetscInt*),(ds,tol,its));
862:   PetscFunctionReturn(PETSC_SUCCESS);
863: }

865: static PetscErrorCode DSNEPSetIntegrationPoints_NEP(DS ds,PetscInt ip)
866: {
867:   DS_NEP         *ctx = (DS_NEP*)ds->data;

869:   PetscFunctionBegin;
870:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) ctx->nnod = 64;
871:   else {
872:     PetscCheck(ip>0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The number of integration points must be > 0");
873:     ctx->nnod = ip;
874:   }
875:   PetscCall(PetscLayoutDestroy(&ctx->map));  /* need to redistribute at next solve */
876:   PetscFunctionReturn(PETSC_SUCCESS);
877: }

879: /*@
880:    DSNEPSetIntegrationPoints - Sets the number of integration points to be
881:    used in the contour integral method.

883:    Logically Collective

885:    Input Parameters:
886: +  ds - the direct solver context
887: -  ip - the number of integration points

889:    Options Database Key:
890: .  -ds_nep_integration_points <ip> - sets the number of integration points

892:    Notes:
893:    This parameter is relevant only in the contour integral method.

895:    Level: advanced

897: .seealso: DSNEPGetIntegrationPoints()
898: @*/
899: PetscErrorCode DSNEPSetIntegrationPoints(DS ds,PetscInt ip)
900: {
901:   PetscFunctionBegin;
904:   PetscTryMethod(ds,"DSNEPSetIntegrationPoints_C",(DS,PetscInt),(ds,ip));
905:   PetscFunctionReturn(PETSC_SUCCESS);
906: }

908: static PetscErrorCode DSNEPGetIntegrationPoints_NEP(DS ds,PetscInt *ip)
909: {
910:   DS_NEP *ctx = (DS_NEP*)ds->data;

912:   PetscFunctionBegin;
913:   *ip = ctx->nnod;
914:   PetscFunctionReturn(PETSC_SUCCESS);
915: }

917: /*@
918:    DSNEPGetIntegrationPoints - Returns the number of integration points used
919:    in the contour integral method.

921:    Not Collective

923:    Input Parameter:
924: .  ds - the direct solver context

926:    Output Parameters:
927: .  ip - the number of integration points

929:    Level: advanced

931: .seealso: DSNEPSetIntegrationPoints()
932: @*/
933: PetscErrorCode DSNEPGetIntegrationPoints(DS ds,PetscInt *ip)
934: {
935:   PetscFunctionBegin;
938:   PetscUseMethod(ds,"DSNEPGetIntegrationPoints_C",(DS,PetscInt*),(ds,ip));
939:   PetscFunctionReturn(PETSC_SUCCESS);
940: }

942: static PetscErrorCode DSNEPSetSamplingSize_NEP(DS ds,PetscInt p)
943: {
944:   DS_NEP *ctx = (DS_NEP*)ds->data;

946:   PetscFunctionBegin;
947:   if (p == PETSC_DECIDE || p == PETSC_DEFAULT) ctx->spls = 0;
948:   else {
949:     PetscCheck(p>=20,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The sample size cannot be smaller than 20");
950:     ctx->spls = p;
951:   }
952:   PetscFunctionReturn(PETSC_SUCCESS);
953: }

955: /*@
956:    DSNEPSetSamplingSize - Sets the number of sampling columns to be
957:    used in the contour integral method.

959:    Logically Collective

961:    Input Parameters:
962: +  ds - the direct solver context
963: -  p - the number of columns for the sampling matrix

965:    Options Database Key:
966: .  -ds_nep_sampling_size <p> - sets the number of sampling columns

968:    Notes:
969:    This parameter is relevant only in the contour integral method.

971:    Level: advanced

973: .seealso: DSNEPGetSamplingSize()
974: @*/
975: PetscErrorCode DSNEPSetSamplingSize(DS ds,PetscInt p)
976: {
977:   PetscFunctionBegin;
980:   PetscTryMethod(ds,"DSNEPSetSamplingSize_C",(DS,PetscInt),(ds,p));
981:   PetscFunctionReturn(PETSC_SUCCESS);
982: }

984: static PetscErrorCode DSNEPGetSamplingSize_NEP(DS ds,PetscInt *p)
985: {
986:   DS_NEP *ctx = (DS_NEP*)ds->data;

988:   PetscFunctionBegin;
989:   *p = ctx->spls;
990:   PetscFunctionReturn(PETSC_SUCCESS);
991: }

993: /*@
994:    DSNEPGetSamplingSize - Returns the number of sampling columns used
995:    in the contour integral method.

997:    Not Collective

999:    Input Parameter:
1000: .  ds - the direct solver context

1002:    Output Parameters:
1003: .  p -  the number of columns for the sampling matrix

1005:    Level: advanced

1007: .seealso: DSNEPSetSamplingSize()
1008: @*/
1009: PetscErrorCode DSNEPGetSamplingSize(DS ds,PetscInt *p)
1010: {
1011:   PetscFunctionBegin;
1014:   PetscUseMethod(ds,"DSNEPGetSamplingSize_C",(DS,PetscInt*),(ds,p));
1015:   PetscFunctionReturn(PETSC_SUCCESS);
1016: }

1018: static PetscErrorCode DSNEPSetComputeMatrixFunction_NEP(DS ds,PetscErrorCode (*fun)(DS,PetscScalar,PetscBool,DSMatType,void*),void *ctx)
1019: {
1020:   DS_NEP *dsctx = (DS_NEP*)ds->data;

1022:   PetscFunctionBegin;
1023:   dsctx->computematrix    = fun;
1024:   dsctx->computematrixctx = ctx;
1025:   PetscFunctionReturn(PETSC_SUCCESS);
1026: }

1028: /*@C
1029:    DSNEPSetComputeMatrixFunction - Sets a user-provided subroutine to compute
1030:    the matrices T(lambda) or T'(lambda).

1032:    Logically Collective

1034:    Input Parameters:
1035: +  ds  - the direct solver context
1036: .  fun - a pointer to the user function
1037: -  ctx - a context pointer (the last parameter to the user function)

1039:    Calling sequence of fun:
1040: $  PetscErrorCode fun(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat,void *ctx)
1041: +   ds     - the direct solver object
1042: .   lambda - point where T(lambda) or T'(lambda) must be evaluated
1043: .   deriv  - if true compute T'(lambda), otherwise compute T(lambda)
1044: .   mat    - the DS matrix where the result must be stored
1045: -   ctx    - optional context, as set by DSNEPSetComputeMatrixFunction()

1047:    Note:
1048:    The result is computed as T(lambda) = sum_i E_i*f_i(lambda), and similarly
1049:    for the derivative.

1051:    Level: developer

1053: .seealso: DSNEPGetComputeMatrixFunction()
1054: @*/
1055: PetscErrorCode DSNEPSetComputeMatrixFunction(DS ds,PetscErrorCode (*fun)(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat,void *ctx),void *ctx)
1056: {
1057:   PetscFunctionBegin;
1059:   PetscTryMethod(ds,"DSNEPSetComputeMatrixFunction_C",(DS,PetscErrorCode (*)(DS,PetscScalar,PetscBool,DSMatType,void*),void*),(ds,fun,ctx));
1060:   PetscFunctionReturn(PETSC_SUCCESS);
1061: }

1063: static PetscErrorCode DSNEPGetComputeMatrixFunction_NEP(DS ds,PetscErrorCode (**fun)(DS,PetscScalar,PetscBool,DSMatType,void*),void **ctx)
1064: {
1065:   DS_NEP *dsctx = (DS_NEP*)ds->data;

1067:   PetscFunctionBegin;
1068:   if (fun) *fun = dsctx->computematrix;
1069:   if (ctx) *ctx = dsctx->computematrixctx;
1070:   PetscFunctionReturn(PETSC_SUCCESS);
1071: }

1073: /*@C
1074:    DSNEPGetComputeMatrixFunction - Returns the user-provided callback function
1075:    set in DSNEPSetComputeMatrixFunction().

1077:    Not Collective

1079:    Input Parameter:
1080: .  ds  - the direct solver context

1082:    Output Parameters:
1083: +  fun - the pointer to the user function
1084: -  ctx - the context pointer

1086:    Level: developer

1088: .seealso: DSNEPSetComputeMatrixFunction()
1089: @*/
1090: PetscErrorCode DSNEPGetComputeMatrixFunction(DS ds,PetscErrorCode (**fun)(DS,PetscScalar,PetscBool,DSMatType,void*),void **ctx)
1091: {
1092:   PetscFunctionBegin;
1094:   PetscUseMethod(ds,"DSNEPGetComputeMatrixFunction_C",(DS,PetscErrorCode (**)(DS,PetscScalar,PetscBool,DSMatType,void*),void**),(ds,fun,ctx));
1095:   PetscFunctionReturn(PETSC_SUCCESS);
1096: }

1098: static PetscErrorCode DSNEPSetRG_NEP(DS ds,RG rg)
1099: {
1100:   DS_NEP         *dsctx = (DS_NEP*)ds->data;

1102:   PetscFunctionBegin;
1103:   PetscCall(PetscObjectReference((PetscObject)rg));
1104:   PetscCall(RGDestroy(&dsctx->rg));
1105:   dsctx->rg = rg;
1106:   PetscFunctionReturn(PETSC_SUCCESS);
1107: }

1109: /*@
1110:    DSNEPSetRG - Associates a region object to the DSNEP solver.

1112:    Collective

1114:    Input Parameters:
1115: +  ds  - the direct solver context
1116: -  rg  - the region context

1118:    Notes:
1119:    The region is used only in the contour integral method, and
1120:    should enclose the wanted eigenvalues.

1122:    Level: developer

1124: .seealso: DSNEPGetRG()
1125: @*/
1126: PetscErrorCode DSNEPSetRG(DS ds,RG rg)
1127: {
1128:   PetscFunctionBegin;
1130:   if (rg) {
1132:     PetscCheckSameComm(ds,1,rg,2);
1133:   }
1134:   PetscTryMethod(ds,"DSNEPSetRG_C",(DS,RG),(ds,rg));
1135:   PetscFunctionReturn(PETSC_SUCCESS);
1136: }

1138: static PetscErrorCode DSNEPGetRG_NEP(DS ds,RG *rg)
1139: {
1140:   DS_NEP         *ctx = (DS_NEP*)ds->data;

1142:   PetscFunctionBegin;
1143:   if (!ctx->rg) {
1144:     PetscCall(RGCreate(PetscObjectComm((PetscObject)ds),&ctx->rg));
1145:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->rg,(PetscObject)ds,1));
1146:     PetscCall(RGSetOptionsPrefix(ctx->rg,((PetscObject)ds)->prefix));
1147:     PetscCall(RGAppendOptionsPrefix(ctx->rg,"ds_nep_"));
1148:     PetscCall(PetscObjectSetOptions((PetscObject)ctx->rg,((PetscObject)ds)->options));
1149:   }
1150:   *rg = ctx->rg;
1151:   PetscFunctionReturn(PETSC_SUCCESS);
1152: }

1154: /*@
1155:    DSNEPGetRG - Obtain the region object associated to the DSNEP solver.

1157:    Collective

1159:    Input Parameter:
1160: .  ds  - the direct solver context

1162:    Output Parameter:
1163: .  rg  - the region context

1165:    Level: developer

1167: .seealso: DSNEPSetRG()
1168: @*/
1169: PetscErrorCode DSNEPGetRG(DS ds,RG *rg)
1170: {
1171:   PetscFunctionBegin;
1174:   PetscUseMethod(ds,"DSNEPGetRG_C",(DS,RG*),(ds,rg));
1175:   PetscFunctionReturn(PETSC_SUCCESS);
1176: }

1178: PetscErrorCode DSSetFromOptions_NEP(DS ds,PetscOptionItems *PetscOptionsObject)
1179: {
1180:   PetscInt       k;
1181:   PetscBool      flg;
1182: #if defined(PETSC_USE_COMPLEX)
1183:   PetscReal      r;
1184:   PetscBool      flg1;
1185:   DS_NEP         *ctx = (DS_NEP*)ds->data;
1186: #endif

1188:   PetscFunctionBegin;
1189:   PetscOptionsHeadBegin(PetscOptionsObject,"DS NEP Options");

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

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

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

1200: #if defined(PETSC_USE_COMPLEX)
1201:     r = ctx->rtol;
1202:     PetscCall(PetscOptionsReal("-ds_nep_refine_tol","Refinement tolerance","DSNEPSetRefine",ctx->rtol,&r,&flg1));
1203:     k = ctx->Nit;
1204:     PetscCall(PetscOptionsInt("-ds_nep_refine_its","Number of iterative refinement iterations","DSNEPSetRefine",ctx->Nit,&k,&flg));
1205:     if (flg1||flg) PetscCall(DSNEPSetRefine(ds,r,k));

1207:     if (ds->method==1) {
1208:       if (!ctx->rg) PetscCall(DSNEPGetRG(ds,&ctx->rg));
1209:       PetscCall(RGSetFromOptions(ctx->rg));
1210:     }
1211: #endif

1213:   PetscOptionsHeadEnd();
1214:   PetscFunctionReturn(PETSC_SUCCESS);
1215: }

1217: PetscErrorCode DSDestroy_NEP(DS ds)
1218: {
1219:   DS_NEP         *ctx = (DS_NEP*)ds->data;
1220:   PetscInt       i;

1222:   PetscFunctionBegin;
1223:   for (i=0;i<ctx->nf;i++) PetscCall(FNDestroy(&ctx->f[i]));
1224:   PetscCall(RGDestroy(&ctx->rg));
1225:   PetscCall(PetscLayoutDestroy(&ctx->map));
1226:   PetscCall(PetscFree(ds->data));
1227:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",NULL));
1228:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",NULL));
1229:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",NULL));
1230:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",NULL));
1231:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",NULL));
1232:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",NULL));
1233:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",NULL));
1234:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",NULL));
1235:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",NULL));
1236:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",NULL));
1237:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",NULL));
1238:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",NULL));
1239:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",NULL));
1240:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",NULL));
1241:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",NULL));
1242:   PetscFunctionReturn(PETSC_SUCCESS);
1243: }

1245: PetscErrorCode DSMatGetSize_NEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
1246: {
1247:   DS_NEP *ctx = (DS_NEP*)ds->data;

1249:   PetscFunctionBegin;
1250:   *rows = ds->n;
1251:   if (t==DS_MAT_Q || t==DS_MAT_Z || t==DS_MAT_U || t==DS_MAT_V) *rows *= ctx->max_mid;
1252:   *cols = ds->n;
1253:   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;
1254:   PetscFunctionReturn(PETSC_SUCCESS);
1255: }

1257: /*MC
1258:    DSNEP - Dense Nonlinear Eigenvalue Problem.

1260:    Level: beginner

1262:    Notes:
1263:    The problem is expressed as T(lambda)*x = 0, where T(lambda) is a
1264:    parameter-dependent matrix written as T(lambda) = sum_i E_i*f_i(lambda).
1265:    The eigenvalues lambda are the arguments returned by DSSolve()..

1267:    The coefficient matrices E_i are the extra matrices of the DS, and
1268:    the scalar functions f_i are passed via DSNEPSetFN(). Optionally, a
1269:    callback function to fill the E_i matrices can be set with
1270:    DSNEPSetComputeMatrixFunction().

1272:    Used DS matrices:
1273: +  DS_MAT_Ex - coefficient matrices of the split form of T(lambda)
1274: .  DS_MAT_A  - (workspace) T(lambda) evaluated at a given lambda (SLP only)
1275: .  DS_MAT_B  - (workspace) T'(lambda) evaluated at a given lambda (SLP only)
1276: .  DS_MAT_Q  - (workspace) left Hankel matrix (contour only)
1277: .  DS_MAT_Z  - (workspace) right Hankel matrix (contour only)
1278: .  DS_MAT_U  - (workspace) left singular vectors (contour only)
1279: .  DS_MAT_V  - (workspace) right singular vectors (contour only)
1280: -  DS_MAT_W  - (workspace) auxiliary matrix of size nxn

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

1286: .seealso: DSCreate(), DSSetType(), DSType, DSNEPSetFN(), DSNEPSetComputeMatrixFunction()
1287: M*/
1288: SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS ds)
1289: {
1290:   DS_NEP         *ctx;

1292:   PetscFunctionBegin;
1293:   PetscCall(PetscNew(&ctx));
1294:   ds->data = (void*)ctx;
1295:   ctx->max_mid = 4;
1296:   ctx->nnod    = 64;
1297:   ctx->Nit     = 3;
1298:   ctx->rtol    = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);

1300:   ds->ops->allocate       = DSAllocate_NEP;
1301:   ds->ops->setfromoptions = DSSetFromOptions_NEP;
1302:   ds->ops->view           = DSView_NEP;
1303:   ds->ops->vectors        = DSVectors_NEP;
1304:   ds->ops->solve[0]       = DSSolve_NEP_SLP;
1305: #if defined(PETSC_USE_COMPLEX)
1306:   ds->ops->solve[1]       = DSSolve_NEP_Contour;
1307: #endif
1308:   ds->ops->sort           = DSSort_NEP;
1309: #if !defined(PETSC_HAVE_MPIUNI)
1310:   ds->ops->synchronize    = DSSynchronize_NEP;
1311: #endif
1312:   ds->ops->destroy        = DSDestroy_NEP;
1313:   ds->ops->matgetsize     = DSMatGetSize_NEP;

1315:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",DSNEPSetFN_NEP));
1316:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",DSNEPGetFN_NEP));
1317:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",DSNEPGetNumFN_NEP));
1318:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",DSNEPGetMinimality_NEP));
1319:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",DSNEPSetMinimality_NEP));
1320:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",DSNEPGetRefine_NEP));
1321:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",DSNEPSetRefine_NEP));
1322:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",DSNEPGetIntegrationPoints_NEP));
1323:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",DSNEPSetIntegrationPoints_NEP));
1324:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",DSNEPGetSamplingSize_NEP));
1325:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",DSNEPSetSamplingSize_NEP));
1326:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",DSNEPSetRG_NEP));
1327:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",DSNEPGetRG_NEP));
1328:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",DSNEPSetComputeMatrixFunction_NEP));
1329:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",DSNEPGetComputeMatrixFunction_NEP));
1330:   PetscFunctionReturn(PETSC_SUCCESS);
1331: }