Actual source code: dsnep.c
slepc-3.17.2 2022-08-09
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,¢er,&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: }