Actual source code: dsnep.c

slepc-3.17.2 2022-08-09
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 <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: {
 36:   DS_NEP         *ctx = (DS_NEP*)ds->data;
 37:   PetscScalar    *T,*E,alpha;
 38:   PetscInt       i,ld,n;
 39:   PetscBLASInt   k,inc=1;

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

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

 66:   DSAllocateMat_Private(ds,DS_MAT_X);
 67:   for (i=0;i<ctx->nf;i++) DSAllocateMat_Private(ds,DSMatExtra[i]);
 68:   PetscFree(ds->perm);
 69:   PetscMalloc1(ld*ctx->max_mid,&ds->perm);
 70:   PetscLogObjectMemory((PetscObject)ds,ld*ctx->max_mid*sizeof(PetscInt));
 71:   PetscFunctionReturn(0);
 72: }

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

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

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

122: PetscErrorCode DSSort_NEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *dummy)
123: {
124:   DS_NEP         *ctx = (DS_NEP*)ds->data;
125:   PetscInt       n,l,i,*perm,lds;
126:   PetscScalar    *A;

128:   if (!ds->sc) PetscFunctionReturn(0);
129:   n = ds->n*ctx->max_mid;
130:   lds = ds->ld*ctx->max_mid;
131:   l = ds->l;
132:   A = ds->mat[DS_MAT_A];
133:   perm = ds->perm;
134:   for (i=0;i<n;i++) perm[i] = i;
135:   if (rr) DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
136:   else DSSortEigenvalues_Private(ds,wr,NULL,perm,PETSC_FALSE);
137:   for (i=l;i<ds->t;i++) A[i+i*lds] = wr[perm[i]];
138:   for (i=l;i<ds->t;i++) wr[i] = A[i+i*lds];
139:   /* n != ds->n */
140:   DSPermuteColumns_Private(ds,0,ds->t,ds->n,DS_MAT_X,perm);
141:   PetscFunctionReturn(0);
142: }

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

157:   if (!ds->mat[DS_MAT_A]) DSAllocateMat_Private(ds,DS_MAT_A);
158:   if (!ds->mat[DS_MAT_B]) DSAllocateMat_Private(ds,DS_MAT_B);
159:   if (!ds->mat[DS_MAT_W]) DSAllocateMat_Private(ds,DS_MAT_W);
160:   PetscBLASIntCast(ds->n,&n);
161:   PetscBLASIntCast(ds->ld,&ld);
162: #if defined(PETSC_USE_COMPLEX)
163:   PetscBLASIntCast(2*ds->n+2*ds->n,&lwork);
164:   PetscBLASIntCast(8*ds->n,&lrwork);
165: #else
166:   PetscBLASIntCast(3*ds->n+8*ds->n,&lwork);
167: #endif
168:   DSAllocateWork_Private(ds,lwork,lrwork,0);
169:   alpha = ds->work;
170:   beta = ds->work + ds->n;
171: #if defined(PETSC_USE_COMPLEX)
172:   work = ds->work + 2*ds->n;
173:   lwork -= 2*ds->n;
174: #else
175:   alphai = ds->work + 2*ds->n;
176:   work = ds->work + 3*ds->n;
177:   lwork -= 3*ds->n;
178: #endif
179:   A = ds->mat[DS_MAT_A];
180:   B = ds->mat[DS_MAT_B];
181:   W = ds->mat[DS_MAT_W];
182:   X = ds->mat[DS_MAT_X];

184:   sigma = 0.0;
185:   if (ds->sc->comparison==SlepcCompareTargetMagnitude || ds->sc->comparison==SlepcCompareTargetReal) sigma = *(PetscScalar*)ds->sc->comparisonctx;
186:   lambda = sigma;
187:   tol = n*PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);

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

191:     /* evaluate T and T' */
192:     DSNEPComputeMatrix(ds,lambda,PETSC_FALSE,DS_MAT_A);
193:     if (it) {
194:       PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&n,&n,&sone,A,&ld,X,&one,&szero,X+ld,&one));
195:       norm = BLASnrm2_(&n,X+ld,&one);
196:       if (norm/PetscAbsScalar(lambda)<=tol) break;
197:     }
198:     DSNEPComputeMatrix(ds,lambda,PETSC_TRUE,DS_MAT_B);

200:     /* compute eigenvalue correction mu and eigenvector u */
201: #if defined(PETSC_USE_COMPLEX)
202:     rwork = ds->rwork;
203:     PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,beta,NULL,&ld,W,&ld,work,&lwork,rwork,&info));
204: #else
205:     PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&n,A,&ld,B,&ld,alpha,alphai,beta,NULL,&ld,W,&ld,work,&lwork,&info));
206: #endif
207:     SlepcCheckLapackInfo("ggev",info);

209:     /* find smallest eigenvalue */
210:     j = 0;
211:     if (beta[j]==0.0) re = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
212:     else re = alpha[j]/beta[j];
213: #if !defined(PETSC_USE_COMPLEX)
214:     if (beta[j]==0.0) im = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
215:     else im = alphai[j]/beta[j];
216: #endif
217:     pos = 0;
218:     for (j=1;j<n;j++) {
219:       if (beta[j]==0.0) re2 = (PetscRealPart(alpha[j])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
220:       else re2 = alpha[j]/beta[j];
221: #if !defined(PETSC_USE_COMPLEX)
222:       if (beta[j]==0.0) im2 = (alphai[j]>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
223:       else im2 = alphai[j]/beta[j];
224:       SlepcCompareSmallestMagnitude(re,im,re2,im2,&result,NULL);
225: #else
226:       SlepcCompareSmallestMagnitude(re,0.0,re2,0.0,&result,NULL);
227: #endif
228:       if (result > 0) {
229:         re = re2;
230: #if !defined(PETSC_USE_COMPLEX)
231:         im = im2;
232: #endif
233:         pos = j;
234:       }
235:     }

237: #if !defined(PETSC_USE_COMPLEX)
239: #endif
240:     mu = alpha[pos]/beta[pos];
241:     PetscArraycpy(X,W+pos*ld,n);
242:     norm = BLASnrm2_(&n,X,&one);
243:     PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&one,X,&n,&info));
244:     SlepcCheckLapackInfo("lascl",info);

246:     /* correct eigenvalue approximation */
247:     lambda = lambda - mu;
248:   }

251:   ds->t = 1;
252:   wr[0] = lambda;
253:   if (wi) wi[0] = 0.0;
254:   PetscFunctionReturn(0);
255: }

257: #if defined(PETSC_USE_COMPLEX)
258: /*
259:   Newton refinement for eigenpairs computed with contour integral.
260:   k  - number of eigenpairs to refine
261:   wr - eigenvalues (eigenvectors are stored in DS_MAT_X)
262: */
263: static PetscErrorCode DSNEPNewtonRefine(DS ds,PetscInt k,PetscScalar *wr)
264: {
265:   DS_NEP         *ctx = (DS_NEP*)ds->data;
266:   PetscScalar    *X,*W,*U,*R,sone=1.0,szero=0.0;
267:   PetscReal      norm;
268:   PetscInt       i,j,ii,nwu=0,*p,jstart=0,jend=k;
269:   const PetscInt *range;
270:   PetscBLASInt   n,*perm,info,ld,one=1,n1;
271:   PetscMPIInt    len,size,root;
272:   PetscLayout    map;

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

341: PetscErrorCode DSSolve_NEP_Contour(DS ds,PetscScalar *wr,PetscScalar *wi)
342: {
343:   DS_NEP         *ctx = (DS_NEP*)ds->data;
344:   PetscScalar    *alpha,*beta,*A,*B,*X,*W,*work,*Rc,*R,*w,*z,*zn,*S,*U,*V;
345:   PetscScalar    sone=1.0,szero=0.0,center;
346:   PetscReal      *rwork,norm,radius,vscale,rgscale,*sigma;
347:   PetscBLASInt   info,n,*perm,p,pp,ld,lwork,k_,rk_,colA,rowA,one=1;
348:   PetscInt       mid,lds,nnod=ctx->nnod,k,i,ii,jj,j,s,off,rk,nwu=0,nw,lrwork,*inside,kstart=0,kend=nnod;
349:   PetscMPIInt    len;
350:   PetscBool      isellipse;
351:   PetscRandom    rand;

354:   /* Contour parameters */
355:   PetscObjectTypeCompare((PetscObject)ctx->rg,RGELLIPSE,&isellipse);
357:   RGEllipseGetParameters(ctx->rg,&center,&radius,&vscale);
358:   RGGetScale(ctx->rg,&rgscale);
359:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {
360:     if (!ctx->map) PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)ds),PETSC_DECIDE,ctx->nnod,1,&ctx->map);
361:     PetscLayoutGetRange(ctx->map,&kstart,&kend);
362:   }

364:   if (!ds->mat[DS_MAT_A]) DSAllocateMat_Private(ds,DS_MAT_A); /* size mid*n */
365:   if (!ds->mat[DS_MAT_B]) DSAllocateMat_Private(ds,DS_MAT_B); /* size mid*n */
366:   if (!ds->mat[DS_MAT_W]) DSAllocateMat_Private(ds,DS_MAT_W); /* size mid*n */
367:   if (!ds->mat[DS_MAT_U]) DSAllocateMat_Private(ds,DS_MAT_U); /* size mid*n */
368:   if (!ds->mat[DS_MAT_V]) DSAllocateMat_Private(ds,DS_MAT_V); /* size n */
369:   A = ds->mat[DS_MAT_A];
370:   B = ds->mat[DS_MAT_B];
371:   W = ds->mat[DS_MAT_W];
372:   U = ds->mat[DS_MAT_U];
373:   V = ds->mat[DS_MAT_V];
374:   X = ds->mat[DS_MAT_X];
375:   mid  = ctx->max_mid;
376:   PetscBLASIntCast(ds->n,&n);
377:   p    = n;   /* maximum number of columns for the probing matrix */
378:   PetscBLASIntCast(ds->ld,&ld);
379:   PetscBLASIntCast(mid*n,&rowA);
380:   PetscBLASIntCast(5*rowA,&lwork);
381:   nw   = n*(2*p+7*mid)+3*nnod+2*mid*n*p;
382:   lrwork = mid*n*6+8*n;
383:   DSAllocateWork_Private(ds,nw,lrwork,n+1);

385:   sigma = ds->rwork;
386:   rwork = ds->rwork+mid*n;
387:   perm  = ds->iwork;
388:   z     = ds->work+nwu;    nwu += nnod;         /* quadrature points */
389:   zn    = ds->work+nwu;    nwu += nnod;         /* normalized quadrature points */
390:   w     = ds->work+nwu;    nwu += nnod;         /* quadrature weights */
391:   Rc    = ds->work+nwu;    nwu += n*p;
392:   R     = ds->work+nwu;    nwu += n*p;
393:   alpha = ds->work+nwu;    nwu += mid*n;
394:   beta  = ds->work+nwu;    nwu += mid*n;
395:   S     = ds->work+nwu;    nwu += 2*mid*n*p;
396:   work  = ds->work+nwu;    /*nwu += mid*n*5;*/

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

401:   /* Set random matrix */
402:   PetscRandomCreate(PetscObjectComm((PetscObject)ds),&rand);
403:   PetscRandomSetSeed(rand,0x12345678);
404:   PetscRandomSeed(rand);
405:   for (j=0;j<p;j++)
406:     for (i=0;i<n;i++) PetscRandomGetValue(rand,Rc+i+j*n);
407:   PetscArrayzero(S,2*mid*n*p);
408:   /* Loop of integration points */
409:   for (k=kstart;k<kend;k++) {
410:     PetscInfo(NULL,"Solving integration point %" PetscInt_FMT "\n",k);
411:     PetscArraycpy(R,Rc,p*n);
412:     DSNEPComputeMatrix(ds,z[k],PETSC_FALSE,DS_MAT_V);

414:     /* LU factorization */
415:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
416:     PetscStackCallBLAS("LAPACKgetrf",LAPACKgetrf_(&n,&n,V,&ld,perm,&info));
417:     SlepcCheckLapackInfo("getrf",info);
418:     PetscStackCallBLAS("LAPACKgetrs",LAPACKgetrs_("N",&n,&p,V,&ld,perm,R,&n,&info));
419:     SlepcCheckLapackInfo("getrs",info);
420:     PetscFPTrapPop();

422:     /* Moments computation */
423:     for (s=0;s<2*ctx->max_mid;s++) {
424:       off = s*n*p;
425:       for (j=0;j<p;j++)
426:         for (i=0;i<n;i++) S[off+i+j*n] += w[k]*R[j*n+i];
427:       w[k] *= zn[k];
428:     }
429:   }

431:   if (ds->pmode==DS_PARALLEL_DISTRIBUTED) {  /* compute final S via reduction */
432:     PetscMPIIntCast(2*mid*n*p,&len);
433:     MPIU_Allreduce(MPI_IN_PLACE,S,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ds));
434:   }
435:   p = ctx->spls?PetscMin(ctx->spls,n):n;
436:   pp = p;
437:   do {
438:     p = pp;
439:     PetscBLASIntCast(mid*p,&colA);

441:     PetscInfo(ds,"Computing SVD of size %" PetscBLASInt_FMT "x%" PetscBLASInt_FMT "\n",rowA,colA);
442:     for (jj=0;jj<mid;jj++) {
443:       for (ii=0;ii<mid;ii++) {
444:         off = jj*p*rowA+ii*n;
445:         for (j=0;j<p;j++)
446:           for (i=0;i<n;i++) A[off+j*rowA+i] = S[((jj+ii)*n+j)*n+i];
447:       }
448:     }
449:     PetscFPTrapPush(PETSC_FP_TRAP_OFF);
450:     PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&rowA,&colA,A,&rowA,sigma,U,&rowA,W,&colA,work,&lwork,rwork,&info));
451:     SlepcCheckLapackInfo("gesvd",info);
452:     PetscFPTrapPop();

454:     rk = colA;
455:     for (i=1;i<colA;i++) if (sigma[i]/sigma[0]<PETSC_MACHINE_EPSILON*1e4) {rk = i; break;}
456:     if (rk<colA || p==n) break;
457:     pp *= 2;
458:   } while (pp<=n);
459:   PetscInfo(ds,"Solving generalized eigenproblem of size %" PetscInt_FMT "\n",rk);
460:   for (jj=0;jj<mid;jj++) {
461:     for (ii=0;ii<mid;ii++) {
462:       off = jj*p*rowA+ii*n;
463:       for (j=0;j<p;j++)
464:         for (i=0;i<n;i++) A[off+j*rowA+i] = S[((jj+ii+1)*n+j)*n+i];
465:     }
466:   }
467:   PetscBLASIntCast(rk,&rk_);
468:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","C",&rowA,&rk_,&colA,&sone,A,&rowA,W,&colA,&szero,B,&rowA));
469:   PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&rk_,&rk_,&rowA,&sone,U,&rowA,B,&rowA,&szero,A,&rk_));
470:   PetscArrayzero(B,n*mid*n*mid);
471:   for (j=0;j<rk;j++) B[j+j*rk_] = sigma[j];
472:   PetscStackCallBLAS("LAPACKggev",LAPACKggev_("N","V",&rk_,A,&rk_,B,&rk_,alpha,beta,NULL,&ld,W,&rk_,work,&lwork,rwork,&info));
473:   for (i=0;i<rk;i++) wr[i] = (center+radius*alpha[i]/beta[i])*rgscale;
474:   PetscMalloc1(rk,&inside);
475:   RGCheckInside(ctx->rg,rk,wr,wi,inside);
476:   k=0;
477:   for (i=0;i<rk;i++)
478:     if (inside[i]==1) inside[k++] = i;
479:   /* Discard values outside region */
480:   lds = ld*mid;
481:   PetscArrayzero(A,lds*lds);
482:   PetscArrayzero(B,lds*lds);
483:   for (i=0;i<k;i++) A[i+i*lds] = (center*beta[inside[i]]+radius*alpha[inside[i]])*rgscale;
484:   for (i=0;i<k;i++) B[i+i*lds] = beta[inside[i]];
485:   for (i=0;i<k;i++) wr[i] = A[i+i*lds]/B[i+i*lds];
486:   for (j=0;j<k;j++) for (i=0;i<rk;i++) W[j*rk+i] = sigma[i]*W[inside[j]*rk+i];
487:   PetscBLASIntCast(k,&k_);
488:   PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&n,&k_,&rk_,&sone,U,&rowA,W,&rk_,&szero,X,&ld));

490:   /* Normalize */
491:   for (j=0;j<k;j++) {
492:     norm = BLASnrm2_(&n,X+ld*j,&one);
493:     for (i=0;i<n;i++) X[ld*j+i] /= norm;
494:   }
495:   PetscFree(inside);
496:   /* Newton refinement */
497:   DSNEPNewtonRefine(ds,k,wr);
498:   ds->t = k;
499:   PetscRandomDestroy(&rand);
500:   PetscFunctionReturn(0);
501: }
502: #endif

504: PetscErrorCode DSSynchronize_NEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
505: {
506:   DS_NEP         *ctx = (DS_NEP*)ds->data;
507:   PetscInt       ld=ds->ld,k=0;
508:   PetscMPIInt    n,n2,rank,size,off=0;

510:   if (!ds->method) { /* SLP */
511:     if (ds->state>=DS_STATE_CONDENSED) k += ds->n;
512:     if (eigr) k += 1;
513:     if (eigi) k += 1;
514:     PetscMPIIntCast(1,&n);
515:     PetscMPIIntCast(ds->n,&n2);
516:   } else { /* Contour */
517:     if (ds->state>=DS_STATE_CONDENSED) k += ctx->max_mid*ds->n*ld;
518:     if (eigr) k += ctx->max_mid*ds->n;
519:     if (eigi) k += ctx->max_mid*ds->n;
520:     PetscMPIIntCast(ctx->max_mid*ds->n,&n);
521:     PetscMPIIntCast(ctx->max_mid*ds->n*ld,&n2);
522:   }
523:   DSAllocateWork_Private(ds,k,0,0);
524:   PetscMPIIntCast(k*sizeof(PetscScalar),&size);
525:   MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
526:   if (!rank) {
527:     if (ds->state>=DS_STATE_CONDENSED) MPI_Pack(ds->mat[DS_MAT_X],n2,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
528:     if (eigr) MPI_Pack(eigr,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
529: #if !defined(PETSC_USE_COMPLEX)
530:     if (eigi) MPI_Pack(eigi,n,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
531: #endif
532:   }
533:   MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
534:   if (rank) {
535:     if (ds->state>=DS_STATE_CONDENSED) MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_X],n2,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
536:     if (eigr) MPI_Unpack(ds->work,size,&off,eigr,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
537: #if !defined(PETSC_USE_COMPLEX)
538:     if (eigi) MPI_Unpack(ds->work,size,&off,eigi,n,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
539: #endif
540:   }
541:   PetscFunctionReturn(0);
542: }

544: static PetscErrorCode DSNEPSetFN_NEP(DS ds,PetscInt n,FN fn[])
545: {
546:   DS_NEP         *ctx = (DS_NEP*)ds->data;
547:   PetscInt       i;

551:   if (ds->ld) PetscInfo(ds,"DSNEPSetFN() called after DSAllocate()\n");
552:   for (i=0;i<n;i++) PetscObjectReference((PetscObject)fn[i]);
553:   for (i=0;i<ctx->nf;i++) FNDestroy(&ctx->f[i]);
554:   for (i=0;i<n;i++) ctx->f[i] = fn[i];
555:   ctx->nf = n;
556:   PetscFunctionReturn(0);
557: }

559: /*@
560:    DSNEPSetFN - Sets a number of functions that define the nonlinear
561:    eigenproblem.

563:    Collective on ds

565:    Input Parameters:
566: +  ds - the direct solver context
567: .  n  - number of functions
568: -  fn - array of functions

570:    Notes:
571:    The nonlinear eigenproblem is defined in terms of the split nonlinear
572:    operator T(lambda) = sum_i A_i*f_i(lambda).

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

578:    Level: advanced

580: .seealso: DSNEPGetFN(), DSAllocate()
581:  @*/
582: PetscErrorCode DSNEPSetFN(DS ds,PetscInt n,FN fn[])
583: {
584:   PetscInt       i;

589:   for (i=0;i<n;i++) {
592:   }
593:   PetscTryMethod(ds,"DSNEPSetFN_C",(DS,PetscInt,FN[]),(ds,n,fn));
594:   PetscFunctionReturn(0);
595: }

597: static PetscErrorCode DSNEPGetFN_NEP(DS ds,PetscInt k,FN *fn)
598: {
599:   DS_NEP *ctx = (DS_NEP*)ds->data;

602:   *fn = ctx->f[k];
603:   PetscFunctionReturn(0);
604: }

606: /*@
607:    DSNEPGetFN - Gets the functions associated with the nonlinear DS.

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

611:    Input Parameters:
612: +  ds - the direct solver context
613: -  k  - the index of the requested function (starting in 0)

615:    Output Parameter:
616: .  fn - the function

618:    Level: advanced

620: .seealso: DSNEPSetFN()
621: @*/
622: PetscErrorCode DSNEPGetFN(DS ds,PetscInt k,FN *fn)
623: {
626:   PetscUseMethod(ds,"DSNEPGetFN_C",(DS,PetscInt,FN*),(ds,k,fn));
627:   PetscFunctionReturn(0);
628: }

630: static PetscErrorCode DSNEPGetNumFN_NEP(DS ds,PetscInt *n)
631: {
632:   DS_NEP *ctx = (DS_NEP*)ds->data;

634:   *n = ctx->nf;
635:   PetscFunctionReturn(0);
636: }

638: /*@
639:    DSNEPGetNumFN - Returns the number of functions stored internally by
640:    the DS.

642:    Not collective

644:    Input Parameter:
645: .  ds - the direct solver context

647:    Output Parameters:
648: .  n - the number of functions passed in DSNEPSetFN()

650:    Level: advanced

652: .seealso: DSNEPSetFN()
653: @*/
654: PetscErrorCode DSNEPGetNumFN(DS ds,PetscInt *n)
655: {
658:   PetscUseMethod(ds,"DSNEPGetNumFN_C",(DS,PetscInt*),(ds,n));
659:   PetscFunctionReturn(0);
660: }

662: static PetscErrorCode DSNEPSetMinimality_NEP(DS ds,PetscInt n)
663: {
664:   DS_NEP *ctx = (DS_NEP*)ds->data;

666:   if (n == PETSC_DECIDE || n == PETSC_DEFAULT) ctx->max_mid = 4;
667:   else {
669:     ctx->max_mid = n;
670:   }
671:   PetscFunctionReturn(0);
672: }

674: /*@
675:    DSNEPSetMinimality - Sets the maximum minimality index used internally by
676:    the DSNEP.

678:    Logically Collective on ds

680:    Input Parameters:
681: +  ds - the direct solver context
682: -  n  - the maximum minimality index

684:    Options Database Key:
685: .  -ds_nep_minimality <n> - sets the maximum minimality index

687:    Notes:
688:    The maximum minimality index is used only in the contour integral method,
689:    and is related to the highest momemts used in the method. The default
690:    value is 1, an larger value might give better accuracy in some cases, but
691:    at a higher cost.

693:    Level: advanced

695: .seealso: DSNEPGetMinimality()
696: @*/
697: PetscErrorCode DSNEPSetMinimality(DS ds,PetscInt n)
698: {
701:   PetscTryMethod(ds,"DSNEPSetMinimality_C",(DS,PetscInt),(ds,n));
702:   PetscFunctionReturn(0);
703: }

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

709:   *n = ctx->max_mid;
710:   PetscFunctionReturn(0);
711: }

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

717:    Not collective

719:    Input Parameter:
720: .  ds - the direct solver context

722:    Output Parameters:
723: .  n - the maximum minimality index passed in DSNEPSetMinimality()

725:    Level: advanced

727: .seealso: DSNEPSetMinimality()
728: @*/
729: PetscErrorCode DSNEPGetMinimality(DS ds,PetscInt *n)
730: {
733:   PetscUseMethod(ds,"DSNEPGetMinimality_C",(DS,PetscInt*),(ds,n));
734:   PetscFunctionReturn(0);
735: }

737: static PetscErrorCode DSNEPSetRefine_NEP(DS ds,PetscReal tol,PetscInt its)
738: {
739:   DS_NEP *ctx = (DS_NEP*)ds->data;

741:   if (tol == PETSC_DEFAULT) ctx->rtol = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);
742:   else {
744:     ctx->rtol = tol;
745:   }
746:   if (its == PETSC_DECIDE || its == PETSC_DEFAULT) ctx->Nit = 3;
747:   else {
749:     ctx->Nit = its;
750:   }
751:   PetscFunctionReturn(0);
752: }

754: /*@
755:    DSNEPSetRefine - Sets the tolerance and the number of iterations of Newton iterative
756:    refinement for eigenpairs.

758:    Logically Collective on ds

760:    Input Parameters:
761: +  ds  - the direct solver context
762: .  tol - the tolerance
763: -  its - the number of iterations

765:    Options Database Key:
766: +  -ds_nep_refine_tol <tol> - sets the tolerance
767: -  -ds_nep_refine_its <its> - sets the number of Newton iterations

769:    Notes:
770:    Iterative refinement of eigenpairs is currently used only in the contour
771:    integral method.

773:    Level: advanced

775: .seealso: DSNEPGetRefine()
776: @*/
777: PetscErrorCode DSNEPSetRefine(DS ds,PetscReal tol,PetscInt its)
778: {
782:   PetscTryMethod(ds,"DSNEPSetRefine_C",(DS,PetscReal,PetscInt),(ds,tol,its));
783:   PetscFunctionReturn(0);
784: }

786: static PetscErrorCode DSNEPGetRefine_NEP(DS ds,PetscReal *tol,PetscInt *its)
787: {
788:   DS_NEP *ctx = (DS_NEP*)ds->data;

790:   if (tol) *tol = ctx->rtol;
791:   if (its) *its = ctx->Nit;
792:   PetscFunctionReturn(0);
793: }

795: /*@
796:    DSNEPGetRefine - Returns the tolerance and the number of iterations of Newton iterative
797:    refinement for eigenpairs.

799:    Not collective

801:    Input Parameter:
802: .  ds - the direct solver context

804:    Output Parameters:
805: +  tol - the tolerance
806: -  its - the number of iterations

808:    Level: advanced

810: .seealso: DSNEPSetRefine()
811: @*/
812: PetscErrorCode DSNEPGetRefine(DS ds,PetscReal *tol,PetscInt *its)
813: {
815:   PetscUseMethod(ds,"DSNEPGetRefine_C",(DS,PetscReal*,PetscInt*),(ds,tol,its));
816:   PetscFunctionReturn(0);
817: }

819: static PetscErrorCode DSNEPSetIntegrationPoints_NEP(DS ds,PetscInt ip)
820: {
821:   DS_NEP         *ctx = (DS_NEP*)ds->data;

823:   if (ip == PETSC_DECIDE || ip == PETSC_DEFAULT) ctx->nnod = 64;
824:   else {
826:     ctx->nnod = ip;
827:   }
828:   PetscLayoutDestroy(&ctx->map);  /* need to redistribute at next solve */
829:   PetscFunctionReturn(0);
830: }

832: /*@
833:    DSNEPSetIntegrationPoints - Sets the number of integration points to be
834:    used in the contour integral method.

836:    Logically Collective on ds

838:    Input Parameters:
839: +  ds - the direct solver context
840: -  ip - the number of integration points

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

845:    Notes:
846:    This parameter is relevant only in the contour integral method.

848:    Level: advanced

850: .seealso: DSNEPGetIntegrationPoints()
851: @*/
852: PetscErrorCode DSNEPSetIntegrationPoints(DS ds,PetscInt ip)
853: {
856:   PetscTryMethod(ds,"DSNEPSetIntegrationPoints_C",(DS,PetscInt),(ds,ip));
857:   PetscFunctionReturn(0);
858: }

860: static PetscErrorCode DSNEPGetIntegrationPoints_NEP(DS ds,PetscInt *ip)
861: {
862:   DS_NEP *ctx = (DS_NEP*)ds->data;

864:   *ip = ctx->nnod;
865:   PetscFunctionReturn(0);
866: }

868: /*@
869:    DSNEPGetIntegrationPoints - Returns the number of integration points used
870:    in the contour integral method.

872:    Not collective

874:    Input Parameter:
875: .  ds - the direct solver context

877:    Output Parameters:
878: .  ip - the number of integration points

880:    Level: advanced

882: .seealso: DSNEPSetIntegrationPoints()
883: @*/
884: PetscErrorCode DSNEPGetIntegrationPoints(DS ds,PetscInt *ip)
885: {
888:   PetscUseMethod(ds,"DSNEPGetIntegrationPoints_C",(DS,PetscInt*),(ds,ip));
889:   PetscFunctionReturn(0);
890: }

892: static PetscErrorCode DSNEPSetSamplingSize_NEP(DS ds,PetscInt p)
893: {
894:   DS_NEP *ctx = (DS_NEP*)ds->data;

896:   if (p == PETSC_DECIDE || p == PETSC_DEFAULT) ctx->spls = 0;
897:   else {
899:     ctx->spls = p;
900:   }
901:   PetscFunctionReturn(0);
902: }

904: /*@
905:    DSNEPSetSamplingSize - Sets the number of sampling columns to be
906:    used in the contour integral method.

908:    Logically Collective on ds

910:    Input Parameters:
911: +  ds - the direct solver context
912: -  p - the number of columns for the sampling matrix

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

917:    Notes:
918:    This parameter is relevant only in the contour integral method.

920:    Level: advanced

922: .seealso: DSNEPGetSamplingSize()
923: @*/
924: PetscErrorCode DSNEPSetSamplingSize(DS ds,PetscInt p)
925: {
928:   PetscTryMethod(ds,"DSNEPSetSamplingSize_C",(DS,PetscInt),(ds,p));
929:   PetscFunctionReturn(0);
930: }

932: static PetscErrorCode DSNEPGetSamplingSize_NEP(DS ds,PetscInt *p)
933: {
934:   DS_NEP *ctx = (DS_NEP*)ds->data;

936:   *p = ctx->spls;
937:   PetscFunctionReturn(0);
938: }

940: /*@
941:    DSNEPGetSamplingSize - Returns the number of sampling columns used
942:    in the contour integral method.

944:    Not collective

946:    Input Parameter:
947: .  ds - the direct solver context

949:    Output Parameters:
950: .  p -  the number of columns for the sampling matrix

952:    Level: advanced

954: .seealso: DSNEPSetSamplingSize()
955: @*/
956: PetscErrorCode DSNEPGetSamplingSize(DS ds,PetscInt *p)
957: {
960:   PetscUseMethod(ds,"DSNEPGetSamplingSize_C",(DS,PetscInt*),(ds,p));
961:   PetscFunctionReturn(0);
962: }

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

968:   dsctx->computematrix    = fun;
969:   dsctx->computematrixctx = ctx;
970:   PetscFunctionReturn(0);
971: }

973: /*@C
974:    DSNEPSetComputeMatrixFunction - Sets a user-provided subroutine to compute
975:    the matrices T(lambda) or T'(lambda).

977:    Logically Collective on ds

979:    Input Parameters:
980: +  ds  - the direct solver context
981: .  fun - a pointer to the user function
982: -  ctx - a context pointer (the last parameter to the user function)

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

987: +   ds     - the direct solver object
988: .   lambda - point where T(lambda) or T'(lambda) must be evaluated
989: .   deriv  - if true compute T'(lambda), otherwise compute T(lambda)
990: .   mat    - the DS matrix where the result must be stored
991: -   ctx    - optional context, as set by DSNEPSetComputeMatrixFunction()

993:    Note:
994:    The result is computed as T(lambda) = sum_i E_i*f_i(lambda), and similarly
995:    for the derivative.

997:    Level: developer

999: .seealso: DSNEPGetComputeMatrixFunction()
1000: @*/
1001: PetscErrorCode DSNEPSetComputeMatrixFunction(DS ds,PetscErrorCode (*fun)(DS,PetscScalar,PetscBool,DSMatType,void*),void *ctx)
1002: {
1004:   PetscTryMethod(ds,"DSNEPSetComputeMatrixFunction_C",(DS,PetscErrorCode (*)(DS,PetscScalar,PetscBool,DSMatType,void*),void*),(ds,fun,ctx));
1005:   PetscFunctionReturn(0);
1006: }

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

1012:   if (fun) *fun = dsctx->computematrix;
1013:   if (ctx) *ctx = dsctx->computematrixctx;
1014:   PetscFunctionReturn(0);
1015: }

1017: /*@C
1018:    DSNEPGetComputeMatrixFunction - Returns the user-provided callback function
1019:    set in DSNEPSetComputeMatrixFunction().

1021:    Not Collective

1023:    Input Parameter:
1024: .  ds  - the direct solver context

1026:    Output Parameters:
1027: +  fun - the pointer to the user function
1028: -  ctx - the context pointer

1030:    Level: developer

1032: .seealso: DSNEPSetComputeMatrixFunction()
1033: @*/
1034: PetscErrorCode DSNEPGetComputeMatrixFunction(DS ds,PetscErrorCode (**fun)(DS,PetscScalar,PetscBool,DSMatType,void*),void **ctx)
1035: {
1037:   PetscUseMethod(ds,"DSNEPGetComputeMatrixFunction_C",(DS,PetscErrorCode (**)(DS,PetscScalar,PetscBool,DSMatType,void*),void**),(ds,fun,ctx));
1038:   PetscFunctionReturn(0);
1039: }

1041: static PetscErrorCode DSNEPSetRG_NEP(DS ds,RG rg)
1042: {
1043:   DS_NEP         *dsctx = (DS_NEP*)ds->data;

1045:   PetscObjectReference((PetscObject)rg);
1046:   RGDestroy(&dsctx->rg);
1047:   dsctx->rg = rg;
1048:   PetscLogObjectParent((PetscObject)ds,(PetscObject)dsctx->rg);
1049:   PetscFunctionReturn(0);
1050: }

1052: /*@
1053:    DSNEPSetRG - Associates a region object to the DSNEP solver.

1055:    Logically Collective on ds

1057:    Input Parameters:
1058: +  ds  - the direct solver context
1059: -  rg  - the region context

1061:    Notes:
1062:    The region is used only in the contour integral method, and
1063:    should enclose the wanted eigenvalues.

1065:    Level: developer

1067: .seealso: DSNEPGetRG()
1068: @*/
1069: PetscErrorCode DSNEPSetRG(DS ds,RG rg)
1070: {
1072:   if (rg) {
1075:   }
1076:   PetscTryMethod(ds,"DSNEPSetRG_C",(DS,RG),(ds,rg));
1077:   PetscFunctionReturn(0);
1078: }

1080: static PetscErrorCode DSNEPGetRG_NEP(DS ds,RG *rg)
1081: {
1082:   DS_NEP         *ctx = (DS_NEP*)ds->data;

1084:   if (!ctx->rg) {
1085:     RGCreate(PetscObjectComm((PetscObject)ds),&ctx->rg);
1086:     PetscObjectIncrementTabLevel((PetscObject)ctx->rg,(PetscObject)ds,1);
1087:     RGSetOptionsPrefix(ctx->rg,((PetscObject)ds)->prefix);
1088:     RGAppendOptionsPrefix(ctx->rg,"ds_nep_");
1089:     PetscLogObjectParent((PetscObject)ds,(PetscObject)ctx->rg);
1090:     PetscObjectSetOptions((PetscObject)ctx->rg,((PetscObject)ds)->options);
1091:   }
1092:   *rg = ctx->rg;
1093:   PetscFunctionReturn(0);
1094: }

1096: /*@
1097:    DSNEPGetRG - Obtain the region object associated to the DSNEP solver.

1099:    Not Collective

1101:    Input Parameter:
1102: .  ds  - the direct solver context

1104:    Output Parameter:
1105: .  rg  - the region context

1107:    Level: developer

1109: .seealso: DSNEPSetRG()
1110: @*/
1111: PetscErrorCode DSNEPGetRG(DS ds,RG *rg)
1112: {
1115:   PetscUseMethod(ds,"DSNEPGetRG_C",(DS,RG*),(ds,rg));
1116:   PetscFunctionReturn(0);
1117: }

1119: PetscErrorCode DSSetFromOptions_NEP(PetscOptionItems *PetscOptionsObject,DS ds)
1120: {
1121:   PetscInt       k;
1122:   PetscBool      flg;
1123: #if defined(PETSC_USE_COMPLEX)
1124:   PetscReal      r;
1125:   PetscBool      flg1;
1126:   DS_NEP         *ctx = (DS_NEP*)ds->data;
1127: #endif

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

1131:     PetscOptionsInt("-ds_nep_minimality","Maximum minimality index","DSNEPSetMinimality",4,&k,&flg);
1132:     if (flg) DSNEPSetMinimality(ds,k);

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

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

1140: #if defined(PETSC_USE_COMPLEX)
1141:     r = ctx->rtol;
1142:     PetscOptionsReal("-ds_nep_refine_tol","Refinement tolerance","DSNEPSetRefine",ctx->rtol,&r,&flg1);
1143:     k = ctx->Nit;
1144:     PetscOptionsInt("-ds_nep_refine_its","Number of iterative refinement iterations","DSNEPSetRefine",ctx->Nit,&k,&flg);
1145:     if (flg1||flg) DSNEPSetRefine(ds,r,k);

1147:     if (ds->method==1) {
1148:       if (!ctx->rg) DSNEPGetRG(ds,&ctx->rg);
1149:       RGSetFromOptions(ctx->rg);
1150:     }
1151: #endif

1153:   PetscOptionsTail();
1154:   PetscFunctionReturn(0);
1155: }

1157: PetscErrorCode DSDestroy_NEP(DS ds)
1158: {
1159:   DS_NEP         *ctx = (DS_NEP*)ds->data;
1160:   PetscInt       i;

1162:   for (i=0;i<ctx->nf;i++) FNDestroy(&ctx->f[i]);
1163:   RGDestroy(&ctx->rg);
1164:   PetscLayoutDestroy(&ctx->map);
1165:   PetscFree(ds->data);
1166:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",NULL);
1167:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",NULL);
1168:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",NULL);
1169:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",NULL);
1170:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",NULL);
1171:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",NULL);
1172:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",NULL);
1173:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",NULL);
1174:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",NULL);
1175:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",NULL);
1176:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",NULL);
1177:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",NULL);
1178:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",NULL);
1179:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",NULL);
1180:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",NULL);
1181:   PetscFunctionReturn(0);
1182: }

1184: PetscErrorCode DSMatGetSize_NEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
1185: {
1186:   DS_NEP *ctx = (DS_NEP*)ds->data;

1188:   *rows = ds->n;
1189:   *cols = ds->n;
1190:   if (t==DS_MAT_X || t==DS_MAT_Y) *cols *= ctx->max_mid;
1191:   PetscFunctionReturn(0);
1192: }

1194: /*MC
1195:    DSNEP - Dense Nonlinear Eigenvalue Problem.

1197:    Level: beginner

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

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

1209:    Used DS matrices:
1210: +  DS_MAT_Ex - coefficient matrices of the split form of T(lambda)
1211: .  DS_MAT_A  - (workspace) T(lambda) evaluated at a given lambda
1212: .  DS_MAT_B  - (workspace) T'(lambda) evaluated at a given lambda
1213: -  DS_MAT_W  - (workspace) eigenvectors of linearization in SLP

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

1219: .seealso: DSCreate(), DSSetType(), DSType, DSNEPSetFN(), DSNEPSetComputeMatrixFunction()
1220: M*/
1221: SLEPC_EXTERN PetscErrorCode DSCreate_NEP(DS ds)
1222: {
1223:   DS_NEP         *ctx;

1225:   PetscNewLog(ds,&ctx);
1226:   ds->data = (void*)ctx;
1227:   ctx->max_mid = 4;
1228:   ctx->nnod    = 64;
1229:   ctx->Nit     = 3;
1230:   ctx->rtol    = PETSC_MACHINE_EPSILON/PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON);

1232:   ds->ops->allocate       = DSAllocate_NEP;
1233:   ds->ops->setfromoptions = DSSetFromOptions_NEP;
1234:   ds->ops->view           = DSView_NEP;
1235:   ds->ops->vectors        = DSVectors_NEP;
1236:   ds->ops->solve[0]       = DSSolve_NEP_SLP;
1237: #if defined(PETSC_USE_COMPLEX)
1238:   ds->ops->solve[1]       = DSSolve_NEP_Contour;
1239: #endif
1240:   ds->ops->sort           = DSSort_NEP;
1241:   ds->ops->synchronize    = DSSynchronize_NEP;
1242:   ds->ops->destroy        = DSDestroy_NEP;
1243:   ds->ops->matgetsize     = DSMatGetSize_NEP;

1245:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetFN_C",DSNEPSetFN_NEP);
1246:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetFN_C",DSNEPGetFN_NEP);
1247:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetNumFN_C",DSNEPGetNumFN_NEP);
1248:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetMinimality_C",DSNEPGetMinimality_NEP);
1249:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetMinimality_C",DSNEPSetMinimality_NEP);
1250:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRefine_C",DSNEPGetRefine_NEP);
1251:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRefine_C",DSNEPSetRefine_NEP);
1252:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetIntegrationPoints_C",DSNEPGetIntegrationPoints_NEP);
1253:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetIntegrationPoints_C",DSNEPSetIntegrationPoints_NEP);
1254:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetSamplingSize_C",DSNEPGetSamplingSize_NEP);
1255:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetSamplingSize_C",DSNEPSetSamplingSize_NEP);
1256:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetRG_C",DSNEPSetRG_NEP);
1257:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetRG_C",DSNEPGetRG_NEP);
1258:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPSetComputeMatrixFunction_C",DSNEPSetComputeMatrixFunction_NEP);
1259:   PetscObjectComposeFunction((PetscObject)ds,"DSNEPGetComputeMatrixFunction_C",DSNEPGetComputeMatrixFunction_NEP);
1260:   PetscFunctionReturn(0);
1261: }