Actual source code: dsnep.c

slepc-3.16.2 2022-02-01
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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 <slepc/private/rgimpl.h>
 13: #include <slepcblaslapack.h>

 15: typedef struct {
 16:   PetscInt       nf;                 /* number of functions in f[] */
 17:   FN             f[DS_NUM_EXTRA];    /* functions defining the nonlinear operator */
 18:   PetscInt       max_mid;            /* maximum minimality index */
 19:   PetscInt       nnod;               /* number of nodes for quadrature rules */
 20:   PetscInt       spls;               /* number of sampling columns for quadrature rules */
 21:   PetscInt       Nit;                /* number of refinement iterations */
 22:   PetscReal      rtol;               /* tolerance of Newton refinement */
 23:   RG             rg;                 /* region for contour integral */
 24:   PetscLayout    map;                /* used to distribute work among MPI processes */
 25:   void           *computematrixctx;
 26:   PetscErrorCode (*computematrix)(DS,PetscScalar,PetscBool,DSMatType,void*);
 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: {
 37:   DS_NEP         *ctx = (DS_NEP*)ds->data;
 38:   PetscScalar    *T,*E,alpha;
 39:   PetscInt       i,ld,n;
 40:   PetscBLASInt   k,inc=1;

 43:   PetscLogEventBegin(DS_Other,ds,0,0,0);
 44:   if (ctx->computematrix) {
 45:     (*ctx->computematrix)(ds,lambda,deriv,mat,ctx->computematrixctx);
 46:   } else {
 47:     DSGetDimensions(ds,&n,NULL,NULL,NULL);
 48:     DSGetLeadingDimension(ds,&ld);
 49:     PetscBLASIntCast(ld*n,&k);
 50:     DSGetArray(ds,mat,&T);
 51:     PetscArrayzero(T,k);
 52:     for (i=0;i<ctx->nf;i++) {
 53:       if (deriv) {
 54:         FNEvaluateDerivative(ctx->f[i],lambda,&alpha);
 55:       } else {
 56:         FNEvaluateFunction(ctx->f[i],lambda,&alpha);
 57:       }
 58:       E = ds->mat[DSMatExtra[i]];
 59:       PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&k,&alpha,E,&inc,T,&inc));
 60:     }
 61:     DSRestoreArray(ds,mat,&T);
 62:   }
 63:   PetscLogEventEnd(DS_Other,ds,0,0,0);
 64:   return(0);
 65: }

 67: PetscErrorCode DSAllocate_NEP(DS ds,PetscInt ld)
 68: {
 70:   DS_NEP         *ctx = (DS_NEP*)ds->data;
 71:   PetscInt       i;

 74:   DSAllocateMat_Private(ds,DS_MAT_X);
 75:   for (i=0;i<ctx->nf;i++) {
 76:     DSAllocateMat_Private(ds,DSMatExtra[i]);
 77:   }
 78:   PetscFree(ds->perm);
 79:   PetscMalloc1(ld*ctx->max_mid,&ds->perm);
 80:   PetscLogObjectMemory((PetscObject)ds,ld*ctx->max_mid*sizeof(PetscInt));
 81:   return(0);
 82: }

 84: PetscErrorCode DSView_NEP(DS ds,PetscViewer viewer)
 85: {
 86:   PetscErrorCode    ierr;
 87:   DS_NEP            *ctx = (DS_NEP*)ds->data;
 88:   PetscViewerFormat format;
 89:   PetscInt          i;
 90:   const char        *methodname[] = {
 91:                      "Successive Linear Problems",
 92:                      "Contour Integral"
 93:   };
 94:   const int         nmeth=sizeof(methodname)/sizeof(methodname[0]);

 97:   PetscViewerGetFormat(viewer,&format);
 98:   if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
 99:     if (ds->method<nmeth) {
100:       PetscViewerASCIIPrintf(viewer,"solving the problem with: %s\n",methodname[ds->method]);
101:     }
102: #if defined(PETSC_USE_COMPLEX)
103:     if (ds->method==1) {  /* contour integral method */
104:       PetscViewerASCIIPrintf(viewer,"number of integration points: %D\n",ctx->nnod);
105:       PetscViewerASCIIPrintf(viewer,"maximum minimality index: %D\n",ctx->max_mid);
106:       if (ctx->spls) { PetscViewerASCIIPrintf(viewer,"number of sampling columns for quadrature: %D\n",ctx->spls); }
107:       if (ctx->Nit) { PetscViewerASCIIPrintf(viewer,"doing iterative refinement (%D its, tolerance %g)\n",ctx->Nit,(double)ctx->rtol); }
108:       RGView(ctx->rg,viewer);
109:     }
110: #endif
111:     if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
112:       PetscViewerASCIIPrintf(viewer,"number of functions: %D\n",ctx->nf);
113:     }
114:     return(0);
115:   }
116:   for (i=0;i<ctx->nf;i++) {
117:     FNView(ctx->f[i],viewer);
118:     DSViewMat(ds,viewer,DSMatExtra[i]);
119:   }
120:   if (ds->state>DS_STATE_INTERMEDIATE) { DSViewMat(ds,viewer,DS_MAT_X); }
121:   return(0);
122: }

124: PetscErrorCode DSVectors_NEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
125: {
127:   if (rnorm) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
128:   switch (mat) {
129:     case DS_MAT_X:
130:       break;
131:     case DS_MAT_Y:
132:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
133:     default:
134:       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
135:   }
136:   return(0);
137: }

139: PetscErrorCode DSSort_NEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *dummy)
140: {
142:   DS_NEP         *ctx = (DS_NEP*)ds->data;
143:   PetscInt       n,l,i,*perm,lds;
144:   PetscScalar    *A;

147:   if (!ds->sc) return(0);
148:   n = ds->n*ctx->max_mid;
149:   lds = ds->ld*ctx->max_mid;
150:   l = ds->l;
151:   A = ds->mat[DS_MAT_A];
152:   perm = ds->perm;
153:   for (i=0;i<n;i++) perm[i] = i;
154:   if (rr) {
155:     DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
156:   } else {
157:     DSSortEigenvalues_Private(ds,wr,NULL,perm,PETSC_FALSE);
158:   }
159:   for (i=l;i<ds->t;i++) A[i+i*lds] = wr[perm[i]];
160:   for (i=l;i<ds->t;i++) wr[i] = A[i+i*lds];
161:   /* n != ds->n */
162:   DSPermuteColumns_Private(ds,0,ds->t,ds->n,DS_MAT_X,perm);
163:   return(0);
164: }

166: PetscErrorCode DSSolve_NEP_SLP(DS ds,PetscScalar *wr,PetscScalar *wi)
167: {
169:   PetscScalar    *A,*B,*W,*X,*work,*alpha,*beta;
170:   PetscScalar    sigma,lambda,mu,re,re2,sone=1.0,szero=0.0;
171:   PetscBLASInt   info,n,ld,lrwork=0,lwork,one=1,zero=0;
172:   PetscInt       it,pos,j,maxit=100,result;
173:   PetscReal      norm,tol,done=1.0;
174: #if defined(PETSC_USE_COMPLEX)
175:   PetscReal      *rwork;
176: #else
177:   PetscReal      *alphai,im,im2;
178: #endif

181:   if (!ds->mat[DS_MAT_A]) {
182:     DSAllocateMat_Private(ds,DS_MAT_A);
183:   }
184:   if (!ds->mat[DS_MAT_B]) {
185:     DSAllocateMat_Private(ds,DS_MAT_B);
186:   }
187:   if (!ds->mat[DS_MAT_W]) {
188:     DSAllocateMat_Private(ds,DS_MAT_W);
189:   }
190:   PetscBLASIntCast(ds->n,&n);
191:   PetscBLASIntCast(ds->ld,&ld);
192: #if defined(PETSC_USE_COMPLEX)
193:   PetscBLASIntCast(2*ds->n+2*ds->n,&lwork);
194:   PetscBLASIntCast(8*ds->n,&lrwork);
195: #else
196:   PetscBLASIntCast(3*ds->n+8*ds->n,&lwork);
197: #endif
198:   DSAllocateWork_Private(ds,lwork,lrwork,0);
199:   alpha = ds->work;
200:   beta = ds->work + ds->n;
201: #if defined(PETSC_USE_COMPLEX)
202:   work = ds->work + 2*ds->n;
203:   lwork -= 2*ds->n;
204: #else
205:   alphai = ds->work + 2*ds->n;
206:   work = ds->work + 3*ds->n;
207:   lwork -= 3*ds->n;
208: #endif
209:   A = ds->mat[DS_MAT_A];
210:   B = ds->mat[DS_MAT_B];
211:   W = ds->mat[DS_MAT_W];
212:   X = ds->mat[DS_MAT_X];

214:   sigma = 0.0;
215:   if (ds->sc->comparison==SlepcCompareTargetMagnitude || ds->sc->comparison==SlepcCompareTargetReal) sigma = *(PetscScalar*)ds->sc->comparisonctx;
216:   lambda = sigma;
217:   tol = n*PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);

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

221:     /* evaluate T and T' */
222:     DSNEPComputeMatrix(ds,lambda,PETSC_FALSE,DS_MAT_A);
223:     if (it) {
224:       PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,A,&ld,X,&one,&szero,X+ld,&one));
225:       norm = BLASnrm2_(&n,X+ld,&one);
226:       if (norm/PetscAbsScalar(lambda)<=tol) break;
227:     }
228:     DSNEPComputeMatrix(ds,lambda,PETSC_TRUE,DS_MAT_B);

230:     /* compute eigenvalue correction mu and eigenvector u */
231: #if defined(PETSC_USE_COMPLEX)
232:     rwork = ds->rwork;
233:     PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,beta,NULL,&ld,W,&ld,work,&lwork,rwork,&info));
234: #else
235:     PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,alphai,beta,NULL,&ld,W,&ld,work,&lwork,&info));
236: #endif
237:     SlepcCheckLapackInfo("ggev",info);

239:     /* find smallest eigenvalue */
240:     j = 0;
241:     if (beta[j]==0.0) re = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
242:     else re = alpha[j]/beta[j];
243: #if !defined(PETSC_USE_COMPLEX)
244:     if (beta[j]==0.0) im = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
245:     else im = alphai[j]/beta[j];
246: #endif
247:     pos = 0;
248:     for (j=1;j<n;j++) {
249:       if (beta[j]==0.0) re2 = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
250:       else re2 = alpha[j]/beta[j];
251: #if !defined(PETSC_USE_COMPLEX)
252:       if (beta[j]==0.0) im2 = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
253:       else im2 = alphai[j]/beta[j];
254:       SlepcCompareSmallestMagnitude(re,im,re2,im2,&result,NULL);
255: #else
256:       SlepcCompareSmallestMagnitude(re,0.0,re2,0.0,&result,NULL);
257: #endif
258:       if (result > 0) {
259:         re = re2;
260: #if !defined(PETSC_USE_COMPLEX)
261:         im = im2;
262: #endif
263:         pos = j;
264:       }
265:     }

267: #if !defined(PETSC_USE_COMPLEX)
268:     if (im!=0.0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"DSNEP found a complex eigenvalue; try rerunning with complex scalars");
269: #endif
270:     mu = alpha[pos]/beta[pos];
271:     PetscArraycpy(X,W+pos*ld,n);
272:     norm = BLASnrm2_(&n,X,&one);
273:     PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&one,X,&n,&info));
274:     SlepcCheckLapackInfo("lascl",info);

276:     /* correct eigenvalue approximation */
277:     lambda = lambda - mu;
278:   }

280:   if (it==maxit) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"DSNEP did not converge");
281:   ds->t = 1;
282:   wr[0] = lambda;
283:   if (wi) wi[0] = 0.0;
284:   return(0);
285: }

287: #if defined(PETSC_USE_COMPLEX)
288: /*
289:   Newton refinement for eigenpairs computed with contour integral.
290:   k  - number of eigenpairs to refine
291:   wr - eigenvalues (eigenvectors are stored in DS_MAT_X)
292: */
293: static PetscErrorCode DSNEPNewtonRefine(DS ds,PetscInt k,PetscScalar *wr)
294: {
296:   DS_NEP         *ctx = (DS_NEP*)ds->data;
297:   PetscScalar    *X,*W,*U,*R,sone=1.0,szero=0.0;
298:   PetscReal      norm;
299:   PetscInt       i,j,ii,nwu=0,*p,jstart=0,jend=k;
300:   const PetscInt *range;
301:   PetscBLASInt   n,*perm,info,ld,one=1,n1;
302:   PetscMPIInt    len,size,root;
303:   PetscLayout    map;

306:   X = ds->mat[DS_MAT_X];
307:   W = ds->mat[DS_MAT_W];
308:   PetscBLASIntCast(ds->n,&n);
309:   PetscBLASIntCast(ds->ld,&ld);
310:   n1 = n+1;
311:   p  = ds->perm;
312:   PetscArrayzero(p,k);
313:   DSAllocateWork_Private(ds,(n+2)*(n+1),0,n+1);
314:   U    = ds->work+nwu;    nwu += (n+1)*(n+1);
315:   R    = ds->work+nwu;    nwu += n+1;
316:   perm = ds->iwork;
317:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {
318:     PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)ds),PETSC_DECIDE,k,1,&map);
319:     PetscLayoutGetRange(map,&jstart,&jend);
320:   }
321:   for (ii=0;ii<ctx->Nit;ii++) {
322:     for (j=jstart;j<jend;j++) {
323:       if (p[j]<2) {
324:         DSNEPComputeMatrix(ds,wr[j],PETSC_FALSE,DS_MAT_W);
325:         PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,W,&ld,X+ld*j,&one,&szero,R,&one));
326:         norm = BLASnrm2_(&n,R,&one);
327:         if (norm/PetscAbsScalar(wr[j]) > ctx->rtol) {
328:           PetscInfo2(NULL,"Refining eigenpair %D, residual=%g\n",j,(double)norm/PetscAbsScalar(wr[j]));
329:           p[j] = 1;
330:           R[n] = 0.0;
331:           for (i=0;i<n;i++) {
332:             PetscArraycpy(U+i*n1,W+i*ld,n);
333:             U[n+i*n1] = PetscConj(X[j*ld+i]);
334:           }
335:           U[n+n*n1] = 0.0;
336:           DSNEPComputeMatrix(ds,wr[j],PETSC_TRUE,DS_MAT_W);
337:           PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,W,&ld,X+ld*j,&one,&szero,U+n*(n+1),&one));
338:           /* solve system  */
339:           PetscFPTrapPush(PETSC_FP_TRAP_OFF);
340:           PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n1,&n1,U,&n1,perm,&info));
341:           SlepcCheckLapackInfo("getrf",info);
342:           PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n1,&one,U,&n1,perm,R,&n1,&info));
343:           SlepcCheckLapackInfo("getrs",info);
344:           PetscFPTrapPop();
345:           wr[j] -= R[n];
346:           for (i=0;i<n;i++) X[j*ld+i] -= R[i];
347:           /* normalization */
348:           norm = BLASnrm2_(&n,X+ld*j,&one);
349:           for (i=0;i<n;i++) X[ld*j+i] /= norm;
350:         } else p[j] = 2;
351:       }
352:     }
353:   }
354:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {  /* communicate results */
355:     PetscMPIIntCast(k,&len);
356: #if defined(PETSC_HAVE_MPI_IN_PLACE)
357:     MPIU_Allreduce(MPI_IN_PLACE,p,len,MPIU_INT,MPIU_SUM,PetscObjectComm((PetscObject)ds));
358: #else
359:     {
360:       PetscInt *buffer;
361:       PetscMalloc1(len,&buffer);
362:       PetscArraycpy(buffer,p,len);
363:       MPIU_Allreduce(buffer,p,len,MPIU_INT,MPIU_SUM,PetscObjectComm((PetscObject)ds));
364:       PetscFree(buffer);
365:     }
366: #endif
367:     MPI_Comm_size(PetscObjectComm((PetscObject)ds),&size);
368:     PetscLayoutGetRanges(map,&range);
369:     for (j=0;j<k;j++) {
370:       if (p[j]) {  /* j-th eigenpair has been refined */
371:         for (root=0;root<size;root++) if (range[root+1]>j) break;
372:         PetscMPIIntCast(1,&len);
373:         MPI_Bcast(wr+j,len,MPIU_SCALAR,root,PetscObjectComm((PetscObject)ds));
374:         PetscMPIIntCast(n,&len);
375:         MPI_Bcast(X+ld*j,len,MPIU_SCALAR,root,PetscObjectComm((PetscObject)ds));
376:       }
377:     }
378:     PetscLayoutDestroy(&map);
379:   }
380:   return(0);
381: }

383: PetscErrorCode DSSolve_NEP_Contour(DS ds,PetscScalar *wr,PetscScalar *wi)
384: {
386:   DS_NEP         *ctx = (DS_NEP*)ds->data;
387:   PetscScalar    *alpha,*beta,*A,*B,*X,*W,*work,*Rc,*R,*w,*z,*zn,*S,*U,*V;
388:   PetscScalar    sone=1.0,szero=0.0,center;
389:   PetscReal      *rwork,norm,radius,vscale,rgscale,*sigma;
390:   PetscBLASInt   info,n,*perm,p,pp,ld,lwork,k_,rk_,colA,rowA,one=1;
391:   PetscInt       mid,lds,nnod=ctx->nnod,k,i,ii,jj,j,s,off,rk,nwu=0,nw,lrwork,*inside,kstart=0,kend=nnod;
392:   PetscMPIInt    len;
393:   PetscBool      isellipse;
394:   PetscRandom    rand;

397:   if (!ctx->rg) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ORDER,"The contour solver requires a region passed with DSNEPSetRG()");
398:   /* Contour parameters */
399:   PetscObjectTypeCompare((PetscObject)ctx->rg,RGELLIPSE,&isellipse);
400:   if (isellipse) {
401:     RGEllipseGetParameters(ctx->rg,&center,&radius,&vscale);
402:   } else SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Region must be Ellipse");
403:   RGGetScale(ctx->rg,&rgscale);
404:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {
405:     if (!ctx->map) { PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)ds),PETSC_DECIDE,ctx->nnod,1,&ctx->map); }
406:     PetscLayoutGetRange(ctx->map,&kstart,&kend);
407:   }

409:   if (!ds->mat[DS_MAT_A]) {
410:     DSAllocateMat_Private(ds,DS_MAT_A);
411:   } /* size mid*n */
412:   if (!ds->mat[DS_MAT_B]) {
413:     DSAllocateMat_Private(ds,DS_MAT_B);
414:   } /* size mid*n */
415:   if (!ds->mat[DS_MAT_W]) {
416:     DSAllocateMat_Private(ds,DS_MAT_W);
417:   } /* size mid*n */
418:   if (!ds->mat[DS_MAT_U]) {
419:     DSAllocateMat_Private(ds,DS_MAT_U);
420:   } /* size mid*n */
421:   if (!ds->mat[DS_MAT_V]) {
422:     DSAllocateMat_Private(ds,DS_MAT_V);
423:   } /* size n */
424:   A = ds->mat[DS_MAT_A];
425:   B = ds->mat[DS_MAT_B];
426:   W = ds->mat[DS_MAT_W];
427:   U = ds->mat[DS_MAT_U];
428:   V = ds->mat[DS_MAT_V];
429:   X = ds->mat[DS_MAT_X];
430:   mid  = ctx->max_mid;
431:   PetscBLASIntCast(ds->n,&n);
432:   p    = n;   /* maximum number of columns for the probing matrix */
433:   PetscBLASIntCast(ds->ld,&ld);
434:   PetscBLASIntCast(mid*n,&rowA);
435:   PetscBLASIntCast(5*rowA,&lwork);
436:   nw   = n*(2*p+7*mid)+3*nnod+2*mid*n*p;
437:   lrwork = mid*n*6+8*n;
438:   DSAllocateWork_Private(ds,nw,lrwork,n+1);

440:   sigma = ds->rwork;
441:   rwork = ds->rwork+mid*n;
442:   perm  = ds->iwork;
443:   z     = ds->work+nwu;    nwu += nnod;         /* quadrature points */
444:   zn    = ds->work+nwu;    nwu += nnod;         /* normalized quadrature points */
445:   w     = ds->work+nwu;    nwu += nnod;         /* quadrature weights */
446:   Rc    = ds->work+nwu;    nwu += n*p;
447:   R     = ds->work+nwu;    nwu += n*p;
448:   alpha = ds->work+nwu;    nwu += mid*n;
449:   beta  = ds->work+nwu;    nwu += mid*n;
450:   S     = ds->work+nwu;    nwu += 2*mid*n*p;
451:   work  = ds->work+nwu;    nwu += mid*n*5;

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

456:   /* Set random matrix */
457:   PetscRandomCreate(PetscObjectComm((PetscObject)ds),&rand);
458:   PetscRandomSetSeed(rand,0x12345678);
459:   PetscRandomSeed(rand);
460:   for (j=0;j<p;j++)
461:     for (i=0;i<n;i++) { PetscRandomGetValue(rand,Rc+i+j*n); }
462:   PetscArrayzero(S,2*mid*n*p);
463:   /* Loop of integration points */
464:   for (k=kstart;k<kend;k++) {
465:     PetscInfo1(NULL,"Solving integration point %D\n",k);
466:     PetscArraycpy(R,Rc,p*n);
467:     DSNEPComputeMatrix(ds,z[k],PETSC_FALSE,DS_MAT_V);

469:     /* LU factorization */
470:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
471:     PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n,&n,V,&ld,perm,&info));
472:     SlepcCheckLapackInfo("getrf",info);
473:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&p,V,&ld,perm,R,&n,&info));
474:     SlepcCheckLapackInfo("getrs",info);
475:     PetscFPTrapPop();

477:     /* Moments computation */
478:     for (s=0;s<2*ctx->max_mid;s++) {
479:       off = s*n*p;
480:       for (j=0;j<p;j++)
481:         for (i=0;i<n;i++) S[off+i+j*n] += w[k]*R[j*n+i];
482:       w[k] *= zn[k];
483:     }
484:   }

486:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {  /* compute final S via reduction */
487:     PetscMPIIntCast(2*mid*n*p,&len);
488: #if defined(PETSC_HAVE_MPI_IN_PLACE)
489:     MPIU_Allreduce(MPI_IN_PLACE,S,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ds));
490: #else
491:     {
492:       PetscScalar *buffer;
493:       PetscMalloc1(len,&buffer);
494:       PetscArraycpy(buffer,S,len);
495:       MPIU_Allreduce(buffer,S,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ds));
496:       PetscFree(buffer);
497:     }
498: #endif
499:   }
500:   p = ctx->spls?PetscMin(ctx->spls,n):n;
501:   pp = p;
502:   do {
503:     p = pp;
504:     PetscBLASIntCast(mid*p,&colA);

506:     PetscInfo2(ds,"Computing SVD of size %Dx%D\n",rowA,colA);
507:     for (jj=0;jj<mid;jj++) {
508:       for (ii=0;ii<mid;ii++) {
509:         off = jj*p*rowA+ii*n;
510:         for (j=0;j<p;j++)
511:           for (i=0;i<n;i++) A[off+j*rowA+i] = S[((jj+ii)*n+j)*n+i];
512:       }
513:     }
514:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
515:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&rowA,&colA,A,&rowA,sigma,U,&rowA,W,&colA,work,&lwork,rwork,&info));
516:     SlepcCheckLapackInfo("gesvd",info);
517:     PetscFPTrapPop();

519:     rk = colA;
520:     for (i=1;i<colA;i++) if (sigma[i]/sigma[0]<PETSC_MACHINE_EPSILON*1e4) {rk = i; break;}
521:     if (rk<colA || p==n) break;
522:     pp *= 2;
523:   } while (pp<=n);
524:   PetscInfo1(ds,"Solving generalized eigenproblem of size %D\n",rk);
525:   for (jj=0;jj<mid;jj++) {
526:     for (ii=0;ii<mid;ii++) {
527:       off = jj*p*rowA+ii*n;
528:       for (j=0;j<p;j++)
529:         for (i=0;i<n;i++) A[off+j*rowA+i] = S[((jj+ii+1)*n+j)*n+i];
530:     }
531:   }
532:   PetscBLASIntCast(rk,&rk_);
533:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&rowA,&rk_,&colA,&sone,A,&rowA,W,&colA,&szero,B,&rowA));
534:   PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&rk_,&rk_,&rowA,&sone,U,&rowA,B,&rowA,&szero,A,&rk_));
535:   PetscArrayzero(B,n*mid*n*mid);
536:   for (j=0;j<rk;j++) B[j+j*rk_] = sigma[j];
537:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&rk_,A,&rk_,B,&rk_,alpha,beta,NULL,&ld,W,&rk_,work,&lwork,rwork,&info));
538:   for (i=0;i<rk;i++) wr[i] = (center+radius*alpha[i]/beta[i])*rgscale;
539:   PetscMalloc1(rk,&inside);
540:   RGCheckInside(ctx->rg,rk,wr,wi,inside);
541:   k=0;
542:   for (i=0;i<rk;i++)
543:     if (inside[i]==1) inside[k++] = i;
544:   /* Discard values outside region */
545:   lds = ld*mid;
546:   PetscArrayzero(A,lds*lds);
547:   PetscArrayzero(B,lds*lds);
548:   for (i=0;i<k;i++) A[i+i*lds] = (center*beta[inside[i]]+radius*alpha[inside[i]])*rgscale;
549:   for (i=0;i<k;i++) B[i+i*lds] = beta[inside[i]];
550:   for (i=0;i<k;i++) wr[i] = A[i+i*lds]/B[i+i*lds];
551:   for (j=0;j<k;j++) for (i=0;i<rk;i++) W[j*rk+i] = sigma[i]*W[inside[j]*rk+i];
552:   PetscBLASIntCast(k,&k_);
553:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k_,&rk_,&sone,U,&rowA,W,&rk_,&szero,X,&ld));

555:   /* Normalize */
556:   for (j=0;j<k;j++) {
557:     norm = BLASnrm2_(&n,X+ld*j,&one);
558:     for (i=0;i<n;i++) X[ld*j+i] /= norm;
559:   }
560:   PetscFree(inside);
561:   /* Newton refinement */
562:   DSNEPNewtonRefine(ds,k,wr);
563:   ds->t = k;
564:   PetscRandomDestroy(&rand);
565:   return(0);
566: }
567: #endif

569: PetscErrorCode DSSynchronize_NEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
570: {
572:   PetscInt       k=0;
573:   PetscMPIInt    n,rank,size,off=0;

576:   if (ds->state>=DS_STATE_CONDENSED) k += ds->n;
577:   if (eigr) k += 1;
578:   if (eigi) k += 1;
579:   DSAllocateWork_Private(ds,k,0,0);
580:   PetscMPIIntCast(k*sizeof(PetscScalar),&size);
581:   PetscMPIIntCast(ds->n,&n);
582:   MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
583:   if (!rank) {
584:     if (ds->state>=DS_STATE_CONDENSED) {
585:       MPI_Pack(ds->mat[DS_MAT_X],n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
586:     }
587:     if (eigr) {
588:       MPI_Pack(eigr,1,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
589:     }
590: #if !defined(PETSC_USE_COMPLEX)
591:     if (eigi) {
592:       MPI_Pack(eigi,1,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
593:     }
594: #endif
595:   }
596:   MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
597:   if (rank) {
598:     if (ds->state>=DS_STATE_CONDENSED) {
599:       MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_X],n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
600:     }
601:     if (eigr) {
602:       MPI_Unpack(ds->work,size,&off,eigr,1,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
603:     }
604: #if !defined(PETSC_USE_COMPLEX)
605:     if (eigi) {
606:       MPI_Unpack(ds->work,size,&off,eigi,1,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
607:     }
608: #endif
609:   }
610:   return(0);
611: }

613: static PetscErrorCode DSNEPSetFN_NEP(DS ds,PetscInt n,FN fn[])
614: {
616:   DS_NEP         *ctx = (DS_NEP*)ds->data;
617:   PetscInt       i;

620:   if (n<=0) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more functions, you have %D",n);
621:   if (n>DS_NUM_EXTRA) SETERRQ2(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Too many functions, you specified %D but the limit is %D",n,DS_NUM_EXTRA);
622:   if (ds->ld) { PetscInfo(ds,"DSNEPSetFN() called after DSAllocate()\n"); }
623:   for (i=0;i<n;i++) {
624:     PetscObjectReference((PetscObject)fn[i]);
625:   }
626:   for (i=0;i<ctx->nf;i++) {
627:     FNDestroy(&ctx->f[i]);
628:   }
629:   for (i=0;i<n;i++) ctx->f[i] = fn[i];
630:   ctx->nf = n;
631:   return(0);
632: }

634: /*@
635:    DSNEPSetFN - Sets a number of functions that define the nonlinear
636:    eigenproblem.

638:    Collective on ds

640:    Input Parameters:
641: +  ds - the direct solver context
642: .  n  - number of functions
643: -  fn - array of functions

645:    Notes:
646:    The nonlinear eigenproblem is defined in terms of the split nonlinear
647:    operator T(lambda) = sum_i A_i*f_i(lambda).

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

653:    Level: advanced

655: .seealso: DSNEPGetFN(), DSAllocate()
656:  @*/
657: PetscErrorCode DSNEPSetFN(DS ds,PetscInt n,FN fn[])
658: {
659:   PetscInt       i;

666:   for (i=0;i<n;i++) {
669:   }
670:   PetscTryMethod(ds,"DSNEPSetFN_C",(DS,PetscInt,FN[]),(ds,n,fn));
671:   return(0);
672: }

674: static PetscErrorCode DSNEPGetFN_NEP(DS ds,PetscInt k,FN *fn)
675: {
676:   DS_NEP *ctx = (DS_NEP*)ds->data;

679:   if (k<0 || k>=ctx->nf) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",ctx->nf-1);
680:   *fn = ctx->f[k];
681:   return(0);
682: }

684: /*@
685:    DSNEPGetFN - Gets the functions associated with the nonlinear DS.

687:    Not collective, though parallel FNs are returned if the DS is parallel

689:    Input Parameters:
690: +  ds - the direct solver context
691: -  k  - the index of the requested function (starting in 0)

693:    Output Parameter:
694: .  fn - the function

696:    Level: advanced

698: .seealso: DSNEPSetFN()
699: @*/
700: PetscErrorCode DSNEPGetFN(DS ds,PetscInt k,FN *fn)
701: {

707:   PetscUseMethod(ds,"DSNEPGetFN_C",(DS,PetscInt,FN*),(ds,k,fn));
708:   return(0);
709: }

711: static PetscErrorCode DSNEPGetNumFN_NEP(DS ds,PetscInt *n)
712: {
713:   DS_NEP *ctx = (DS_NEP*)ds->data;

716:   *n = ctx->nf;
717:   return(0);
718: }

720: /*@
721:    DSNEPGetNumFN - Returns the number of functions stored internally by
722:    the DS.

724:    Not collective

726:    Input Parameter:
727: .  ds - the direct solver context

729:    Output Parameters:
730: .  n - the number of functions passed in DSNEPSetFN()

732:    Level: advanced

734: .seealso: DSNEPSetFN()
735: @*/
736: PetscErrorCode DSNEPGetNumFN(DS ds,PetscInt *n)
737: {

743:   PetscUseMethod(ds,"DSNEPGetNumFN_C",(DS,PetscInt*),(ds,n));
744:   return(0);
745: }

747: static PetscErrorCode DSNEPSetMinimality_NEP(DS ds,PetscInt n)
748: {
749:   DS_NEP *ctx = (DS_NEP*)ds->data;

752:   if (n == PETSC_DECIDE || n == PETSC_DEFAULT) ctx->max_mid = 4;
753:   else {
754:     if (n<1) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The minimality value must be > 0");
755:     ctx->max_mid = n;
756:   }
757:   return(0);
758: }

760: /*@
761:    DSNEPSetMinimality - Sets the maximum minimality index used internally by
762:    the DSNEP.

764:    Logically Collective on ds

766:    Input Parameters:
767: +  ds - the direct solver context
768: -  n  - the maximum minimality index

770:    Options Database Key:
771: .  -ds_nep_minimality <n> - sets the maximum minimality index

773:    Notes:
774:    The maximum minimality index is used only in the contour integral method,
775:    and is related to the highest momemts used in the method. The default
776:    value is 1, an larger value might give better accuracy in some cases, but
777:    at a higher cost.

779:    Level: advanced

781: .seealso: DSNEPGetMinimality()
782: @*/
783: PetscErrorCode DSNEPSetMinimality(DS ds,PetscInt n)
784: {

790:   PetscTryMethod(ds,"DSNEPSetMinimality_C",(DS,PetscInt),(ds,n));
791:   return(0);
792: }

794: static PetscErrorCode DSNEPGetMinimality_NEP(DS ds,PetscInt *n)
795: {
796:   DS_NEP *ctx = (DS_NEP*)ds->data;

799:   *n = ctx->max_mid;
800:   return(0);
801: }

803: /*@
804:    DSNEPGetMinimality - Returns the maximum minimality index used internally by
805:    the DSNEP.

807:    Not collective

809:    Input Parameter:
810: .  ds - the direct solver context

812:    Output Parameters:
813: .  n - the maximum minimality index passed in DSNEPSetMinimality()

815:    Level: advanced

817: .seealso: DSNEPSetMinimality()
818: @*/
819: PetscErrorCode DSNEPGetMinimality(DS ds,PetscInt *n)
820: {

826:   PetscUseMethod(ds,"DSNEPGetMinimality_C",(DS,PetscInt*),(ds,n));
827:   return(0);
828: }

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

835:   if (tol == PETSC_DEFAULT) ctx->rtol = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);
836:   else {
837:     if (tol<=0.0) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The tolerance must be > 0");
838:     ctx->rtol = tol;
839:   }
840:   if (its == PETSC_DECIDE || its == PETSC_DEFAULT) ctx->Nit = 3;
841:   else {
842:     if (its<0) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The number of iterations must be >= 0");
843:     ctx->Nit = its;
844:   }
845:   return(0);
846: }

848: /*@
849:    DSNEPSetRefine - Sets the tolerance and the number of iterations of Newton iterative
850:    refinement for eigenpairs.

852:    Logically Collective on ds

854:    Input Parameters:
855: +  ds  - the direct solver context
856: .  tol - the tolerance
857: -  its - the number of iterations

859:    Options Database Key:
860: +  -ds_nep_refine_tol <tol> - sets the tolerance
861: -  -ds_nep_refine_its <its> - sets the number of Newton iterations

863:    Notes:
864:    Iterative refinement of eigenpairs is currently used only in the contour
865:    integral method.

867:    Level: advanced

869: .seealso: DSNEPGetRefine()
870: @*/
871: PetscErrorCode DSNEPSetRefine(DS ds,PetscReal tol,PetscInt its)
872: {

879:   PetscTryMethod(ds,"DSNEPSetRefine_C",(DS,PetscReal,PetscInt),(ds,tol,its));
880:   return(0);
881: }

883: static PetscErrorCode DSNEPGetRefine_NEP(DS ds,PetscReal *tol,PetscInt *its)
884: {
885:   DS_NEP *ctx = (DS_NEP*)ds->data;

888:   if (tol) *tol = ctx->rtol;
889:   if (its) *its = ctx->Nit;
890:   return(0);
891: }

893: /*@
894:    DSNEPGetRefine - Returns the tolerance and the number of iterations of Newton iterative
895:    refinement for eigenpairs.

897:    Not collective

899:    Input Parameter:
900: .  ds - the direct solver context

902:    Output Parameters:
903: +  tol - the tolerance
904: -  its - the number of iterations

906:    Level: advanced

908: .seealso: DSNEPSetRefine()
909: @*/
910: PetscErrorCode DSNEPGetRefine(DS ds,PetscReal *tol,PetscInt *its)
911: {

916:   PetscUseMethod(ds,"DSNEPGetRefine_C",(DS,PetscReal*,PetscInt*),(ds,tol,its));
917:   return(0);
918: }

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

926:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) ctx->nnod = 64;
927:   else {
928:     if (ip<1) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The number of integration points must be > 0");
929:     ctx->nnod = ip;
930:   }
931:   PetscLayoutDestroy(&ctx->map);  /* need to redistribute at next solve */
932:   return(0);
933: }

935: /*@
936:    DSNEPSetIntegrationPoints - Sets the number of integration points to be
937:    used in the contour integral method.

939:    Logically Collective on ds

941:    Input Parameters:
942: +  ds - the direct solver context
943: -  ip - the number of integration points

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

948:    Notes:
949:    This parameter is relevant only in the contour integral method.

951:    Level: advanced

953: .seealso: DSNEPGetIntegrationPoints()
954: @*/
955: PetscErrorCode DSNEPSetIntegrationPoints(DS ds,PetscInt ip)
956: {

962:   PetscTryMethod(ds,"DSNEPSetIntegrationPoints_C",(DS,PetscInt),(ds,ip));
963:   return(0);
964: }

966: static PetscErrorCode DSNEPGetIntegrationPoints_NEP(DS ds,PetscInt *ip)
967: {
968:   DS_NEP *ctx = (DS_NEP*)ds->data;

971:   *ip = ctx->nnod;
972:   return(0);
973: }

975: /*@
976:    DSNEPGetIntegrationPoints - Returns the number of integration points used
977:    in the contour integral method.

979:    Not collective

981:    Input Parameter:
982: .  ds - the direct solver context

984:    Output Parameters:
985: .  ip - the number of integration points

987:    Level: advanced

989: .seealso: DSNEPSetIntegrationPoints()
990: @*/
991: PetscErrorCode DSNEPGetIntegrationPoints(DS ds,PetscInt *ip)
992: {

998:   PetscUseMethod(ds,"DSNEPGetIntegrationPoints_C",(DS,PetscInt*),(ds,ip));
999:   return(0);
1000: }

1002: static PetscErrorCode DSNEPSetSamplingSize_NEP(DS ds,PetscInt p)
1003: {
1004:   DS_NEP *ctx = (DS_NEP*)ds->data;

1007:   if (p == PETSC_DECIDE || p == PETSC_DEFAULT) ctx->spls = 0;
1008:   else {
1009:     if (p<20) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The sample size can not be smaller than 20");
1010:     ctx->spls = p;
1011:   }
1012:   return(0);
1013: }

1015: /*@
1016:    DSNEPSetSamplingSize - Sets the number of sampling columns to be
1017:    used in the contour integral method.

1019:    Logically Collective on ds

1021:    Input Parameters:
1022: +  ds - the direct solver context
1023: -  p - the number of columns for the sampling matrix

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

1028:    Notes:
1029:    This parameter is relevant only in the contour integral method.

1031:    Level: advanced

1033: .seealso: DSNEPGetSamplingSize()
1034: @*/
1035: PetscErrorCode DSNEPSetSamplingSize(DS ds,PetscInt p)
1036: {

1042:   PetscTryMethod(ds,"DSNEPSetSamplingSize_C",(DS,PetscInt),(ds,p));
1043:   return(0);
1044: }

1046: static PetscErrorCode DSNEPGetSamplingSize_NEP(DS ds,PetscInt *p)
1047: {
1048:   DS_NEP *ctx = (DS_NEP*)ds->data;

1051:   *p = ctx->spls;
1052:   return(0);
1053: }

1055: /*@
1056:    DSNEPGetSamplingSize - Returns the number of sampling columns used
1057:    in the contour integral method.

1059:    Not collective

1061:    Input Parameter:
1062: .  ds - the direct solver context

1064:    Output Parameters:
1065: .  p -  the number of columns for the sampling matrix

1067:    Level: advanced

1069: .seealso: DSNEPSetSamplingSize()
1070: @*/
1071: PetscErrorCode DSNEPGetSamplingSize(DS ds,PetscInt *p)
1072: {

1078:   PetscUseMethod(ds,"DSNEPGetSamplingSize_C",(DS,PetscInt*),(ds,p));
1079:   return(0);
1080: }

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

1087:   dsctx->computematrix    = fun;
1088:   dsctx->computematrixctx = ctx;
1089:   return(0);
1090: }

1092: /*@C
1093:    DSNEPSetComputeMatrixFunction - Sets a user-provided subroutine to compute
1094:    the matrices T(lambda) or T'(lambda).

1096:    Logically Collective on ds

1098:    Input Parameters:
1099: +  ds  - the direct solver context
1100: .  fun - a pointer to the user function
1101: -  ctx - a context pointer (the last parameter to the user function)

1103:    Calling Sequence of fun:
1104: $   fun(DS ds,PetscScalar lambda,PetscBool deriv,DSMatType mat,void *ctx)

1106: +   ds     - the direct solver object
1107: .   lambda - point where T(lambda) or T'(lambda) must be evaluated
1108: .   deriv  - if true compute T'(lambda), otherwise compute T(lambda)
1109: .   mat    - the DS matrix where the result must be stored
1110: -   ctx    - optional context, as set by DSNEPSetComputeMatrixFunction()

1112:    Note:
1113:    The result is computed as T(lambda) = sum_i E_i*f_i(lambda), and similarly
1114:    for the derivative.

1116:    Level: developer

1118: .seealso: DSNEPGetComputeMatrixFunction()
1119: @*/
1120: PetscErrorCode DSNEPSetComputeMatrixFunction(DS ds,PetscErrorCode (*fun)(DS,PetscScalar,PetscBool,DSMatType,void*),void *ctx)
1121: {

1126:   PetscTryMethod(ds,"DSNEPSetComputeMatrixFunction_C",(DS,PetscErrorCode (*)(DS,PetscScalar,PetscBool,DSMatType,void*),void*),(ds,fun,ctx));
1127:   return(0);
1128: }

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

1135:   if (fun) *fun = dsctx->computematrix;
1136:   if (ctx) *ctx = dsctx->computematrixctx;
1137:   return(0);
1138: }

1140: /*@C
1141:    DSNEPGetComputeMatrixFunction - Returns the user-provided callback function
1142:    set in DSNEPSetComputeMatrixFunction().

1144:    Not Collective

1146:    Input Parameter:
1147: .  ds  - the direct solver context

1149:    Output Parameters:
1150: +  fun - the pointer to the user function
1151: -  ctx - the context pointer

1153:    Level: developer

1155: .seealso: DSNEPSetComputeMatrixFunction()
1156: @*/
1157: PetscErrorCode DSNEPGetComputeMatrixFunction(DS ds,PetscErrorCode (**fun)(DS,PetscScalar,PetscBool,DSMatType,void*),void **ctx)
1158: {

1163:   PetscUseMethod(ds,"DSNEPGetComputeMatrixFunction_C",(DS,PetscErrorCode (**)(DS,PetscScalar,PetscBool,DSMatType,void*),void**),(ds,fun,ctx));
1164:   return(0);
1165: }

1167: static PetscErrorCode DSNEPSetRG_NEP(DS ds,RG rg)
1168: {
1170:   DS_NEP         *dsctx = (DS_NEP*)ds->data;

1173:   PetscObjectReference((PetscObject)rg);
1174:   RGDestroy(&dsctx->rg);
1175:   dsctx->rg = rg;
1176:   PetscLogObjectParent((PetscObject)ds,(PetscObject)dsctx->rg);
1177:   return(0);
1178: }

1180: /*@
1181:    DSNEPSetRG - Associates a region object to the DSNEP solver.

1183:    Logically Collective on ds

1185:    Input Parameters:
1186: +  ds  - the direct solver context
1187: -  rg  - the region context

1189:    Notes:
1190:    The region is used only in the contour integral method, and
1191:    should enclose the wanted eigenvalues.

1193:    Level: developer

1195: .seealso: DSNEPGetRG()
1196: @*/
1197: PetscErrorCode DSNEPSetRG(DS ds,RG rg)
1198: {

1203:   if (rg) {
1206:   }
1207:   PetscTryMethod(ds,"DSNEPSetRG_C",(DS,RG),(ds,rg));
1208:   return(0);
1209: }

1211: static PetscErrorCode DSNEPGetRG_NEP(DS ds,RG *rg)
1212: {
1214:   DS_NEP         *ctx = (DS_NEP*)ds->data;

1217:   if (!ctx->rg) {
1218:     RGCreate(PetscObjectComm((PetscObject)ds),&ctx->rg);
1219:     PetscObjectIncrementTabLevel((PetscObject)ctx->rg,(PetscObject)ds,1);
1220:     RGSetOptionsPrefix(ctx->rg,((PetscObject)ds)->prefix);
1221:     RGAppendOptionsPrefix(ctx->rg,"ds_nep_");
1222:     PetscLogObjectParent((PetscObject)ds,(PetscObject)ctx->rg);
1223:     PetscObjectSetOptions((PetscObject)ctx->rg,((PetscObject)ds)->options);
1224:   }
1225:   *rg = ctx->rg;
1226:   return(0);
1227: }

1229: /*@
1230:    DSNEPGetRG - Obtain the region object associated to the DSNEP solver.

1232:    Not Collective

1234:    Input Parameter:
1235: .  ds  - the direct solver context

1237:    Output Parameter:
1238: .  rg  - the region context

1240:    Level: developer

1242: .seealso: DSNEPSetRG()
1243: @*/
1244: PetscErrorCode DSNEPGetRG(DS ds,RG *rg)
1245: {

1251:   PetscUseMethod(ds,"DSNEPGetRG_C",(DS,RG*),(ds,rg));
1252:   return(0);
1253: }

1255: PetscErrorCode DSSetFromOptions_NEP(PetscOptionItems *PetscOptionsObject,DS ds)
1256: {
1258:   PetscInt       k;
1259:   PetscBool      flg;
1260: #if defined(PETSC_USE_COMPLEX)
1261:   PetscReal      r;
1262:   PetscBool      flg1;
1263:   DS_NEP         *ctx = (DS_NEP*)ds->data;
1264: #endif

1267:   PetscOptionsHead(PetscOptionsObject,"DS NEP Options");

1269:     PetscOptionsInt("-ds_nep_minimality","Maximum minimality index","DSNEPSetMinimality",4,&k,&flg);
1270:     if (flg) { DSNEPSetMinimality(ds,k); }

1272:     PetscOptionsInt("-ds_nep_integration_points","Number of integration points","DSNEPSetIntegrationPoints",64,&k,&flg);
1273:     if (flg) { DSNEPSetIntegrationPoints(ds,k); }

1275:     PetscOptionsInt("-ds_nep_sampling_size","Number of sampling columns","DSNEPSetSamplingSize",0,&k,&flg);
1276:     if (flg) { DSNEPSetSamplingSize(ds,k); }

1278: #if defined(PETSC_USE_COMPLEX)
1279:     r = ctx->rtol;
1280:     PetscOptionsReal("-ds_nep_refine_tol","Refinement tolerance","DSNEPSetRefine",ctx->rtol,&r,&flg1);
1281:     k = ctx->Nit;
1282:     PetscOptionsInt("-ds_nep_refine_its","Number of iterative refinement iterations","DSNEPSetRefine",ctx->Nit,&k,&flg);
1283:     if (flg1||flg) { DSNEPSetRefine(ds,r,k); }

1285:     if (ds->method==1) {
1286:       if (!ctx->rg) { DSNEPGetRG(ds,&ctx->rg); }
1287:       RGSetFromOptions(ctx->rg);
1288:     }
1289: #endif

1291:   PetscOptionsTail();
1292:   return(0);
1293: }

1295: PetscErrorCode DSDestroy_NEP(DS ds)
1296: {
1298:   DS_NEP         *ctx = (DS_NEP*)ds->data;
1299:   PetscInt       i;

1302:   for (i=0;i<ctx->nf;i++) {
1303:     FNDestroy(&ctx->f[i]);
1304:   }
1305:   RGDestroy(&ctx->rg);
1306:   PetscLayoutDestroy(&ctx->map);
1307:   PetscFree(ds->data);
1308:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",NULL);
1309:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",NULL);
1310:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",NULL);
1311:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",NULL);
1312:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",NULL);
1313:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",NULL);
1314:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",NULL);
1315:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",NULL);
1316:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",NULL);
1317:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",NULL);
1318:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",NULL);
1319:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",NULL);
1320:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",NULL);
1321:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",NULL);
1322:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",NULL);
1323:   return(0);
1324: }

1326: PetscErrorCode DSMatGetSize_NEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
1327: {
1328:   DS_NEP *ctx = (DS_NEP*)ds->data;

1331:   *rows = ds->n;
1332:   *cols = ds->n;
1333:   if (t==DS_MAT_X || t==DS_MAT_Y) *cols *= ctx->max_mid;
1334:   return(0);
1335: }

1337: /*MC
1338:    DSNEP - Dense Nonlinear Eigenvalue Problem.

1340:    Level: beginner

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

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

1352:    Used DS matrices:
1353: +  DS_MAT_Ex - coefficient matrices of the split form of T(lambda)
1354: .  DS_MAT_A  - (workspace) T(lambda) evaluated at a given lambda
1355: .  DS_MAT_B  - (workspace) T'(lambda) evaluated at a given lambda
1356: -  DS_MAT_W  - (workspace) eigenvectors of linearization in SLP

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

1362: .seealso: DSCreate(), DSSetType(), DSType, DSNEPSetFN(), DSNEPSetComputeMatrixFunction()
1363: M*/
1364: SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS ds)
1365: {
1366:   DS_NEP         *ctx;

1370:   PetscNewLog(ds,&ctx);
1371:   ds->data = (void*)ctx;
1372:   ctx->max_mid = 4;
1373:   ctx->nnod    = 64;
1374:   ctx->Nit     = 3;
1375:   ctx->rtol    = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);

1377:   ds->ops->allocate       = DSAllocate_NEP;
1378:   ds->ops->setfromoptions = DSSetFromOptions_NEP;
1379:   ds->ops->view           = DSView_NEP;
1380:   ds->ops->vectors        = DSVectors_NEP;
1381:   ds->ops->solve[0]       = DSSolve_NEP_SLP;
1382: #if defined(PETSC_USE_COMPLEX)
1383:   ds->ops->solve[1]       = DSSolve_NEP_Contour;
1384: #endif
1385:   ds->ops->sort           = DSSort_NEP;
1386:   ds->ops->synchronize    = DSSynchronize_NEP;
1387:   ds->ops->destroy        = DSDestroy_NEP;
1388:   ds->ops->matgetsize     = DSMatGetSize_NEP;

1390:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",DSNEPSetFN_NEP);
1391:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",DSNEPGetFN_NEP);
1392:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",DSNEPGetNumFN_NEP);
1393:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",DSNEPGetMinimality_NEP);
1394:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",DSNEPSetMinimality_NEP);
1395:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",DSNEPGetRefine_NEP);
1396:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",DSNEPSetRefine_NEP);
1397:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",DSNEPGetIntegrationPoints_NEP);
1398:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",DSNEPSetIntegrationPoints_NEP);
1399:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",DSNEPGetSamplingSize_NEP);
1400:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",DSNEPSetSamplingSize_NEP);
1401:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",DSNEPSetRG_NEP);
1402:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",DSNEPGetRG_NEP);
1403:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",DSNEPSetComputeMatrixFunction_NEP);
1404:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",DSNEPGetComputeMatrixFunction_NEP);
1405:   return(0);
1406: }