Actual source code: dspep.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcds.h" I*/
12: #include <slepcblaslapack.h>
14: typedef struct {
15: PetscInt d; /* polynomial degree */
16: PetscReal *pbc; /* polynomial basis coefficients */
17: } DS_PEP;
19: PetscErrorCode DSAllocate_PEP(DS ds,PetscInt ld)
20: {
22: DS_PEP *ctx = (DS_PEP*)ds->data;
23: PetscInt i;
26: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
27: DSAllocateMat_Private(ds,DS_MAT_X);
28: DSAllocateMat_Private(ds,DS_MAT_Y);
29: for (i=0;i<=ctx->d;i++) {
30: DSAllocateMat_Private(ds,DSMatExtra[i]);
31: }
32: PetscFree(ds->perm);
33: PetscMalloc1(ld*ctx->d,&ds->perm);
34: PetscLogObjectMemory((PetscObject)ds,ld*ctx->d*sizeof(PetscInt));
35: return(0);
36: }
38: PetscErrorCode DSView_PEP(DS ds,PetscViewer viewer)
39: {
40: PetscErrorCode ierr;
41: DS_PEP *ctx = (DS_PEP*)ds->data;
42: PetscViewerFormat format;
43: PetscInt i;
46: PetscViewerGetFormat(viewer,&format);
47: PetscViewerASCIIPrintf(viewer,"polynomial degree: %D\n",ctx->d);
48: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
49: for (i=0;i<=ctx->d;i++) {
50: DSViewMat(ds,viewer,DSMatExtra[i]);
51: }
52: if (ds->state>DS_STATE_INTERMEDIATE) {
53: DSViewMat(ds,viewer,DS_MAT_X);
54: }
55: return(0);
56: }
58: PetscErrorCode DSVectors_PEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
59: {
61: if (rnorm) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
62: switch (mat) {
63: case DS_MAT_X:
64: break;
65: case DS_MAT_Y:
66: break;
67: default:
68: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
69: }
70: return(0);
71: }
73: PetscErrorCode DSSort_PEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *kout)
74: {
76: DS_PEP *ctx = (DS_PEP*)ds->data;
77: PetscInt n,i,j,k,p,*perm,told,ld;
78: PetscScalar *A,*X,*Y,rtmp,rtmp2;
81: if (!ds->sc) return(0);
82: n = ds->n*ctx->d;
83: A = ds->mat[DS_MAT_A];
84: perm = ds->perm;
85: for (i=0;i<n;i++) perm[i] = i;
86: told = ds->t;
87: ds->t = n; /* force the sorting routines to consider d*n eigenvalues */
88: if (rr) {
89: DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
90: } else {
91: DSSortEigenvalues_Private(ds,wr,wi,perm,PETSC_FALSE);
92: }
93: ds->t = told; /* restore value of t */
94: for (i=0;i<n;i++) A[i] = wr[perm[i]];
95: for (i=0;i<n;i++) wr[i] = A[i];
96: for (i=0;i<n;i++) A[i] = wi[perm[i]];
97: for (i=0;i<n;i++) wi[i] = A[i];
98: /* cannot use DSPermuteColumns_Private() since matrix is not square */
99: ld = ds->ld;
100: X = ds->mat[DS_MAT_X];
101: Y = ds->mat[DS_MAT_Y];
102: for (i=0;i<n;i++) {
103: p = perm[i];
104: if (p != i) {
105: j = i + 1;
106: while (perm[j] != i) j++;
107: perm[j] = p; perm[i] = i;
108: /* swap columns i and j */
109: for (k=0;k<ds->n;k++) {
110: rtmp = X[k+p*ld]; X[k+p*ld] = X[k+i*ld]; X[k+i*ld] = rtmp;
111: rtmp2 = Y[k+p*ld]; Y[k+p*ld] = Y[k+i*ld]; Y[k+i*ld] = rtmp2;
112: }
113: }
114: }
115: return(0);
116: }
118: PetscErrorCode DSSolve_PEP_QZ(DS ds,PetscScalar *wr,PetscScalar *wi)
119: {
120: #if defined(SLEPC_MISSING_LAPACK_GGEV)
122: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GGEV - Lapack routine is unavailable");
123: #else
125: DS_PEP *ctx = (DS_PEP*)ds->data;
126: PetscInt i,j,k,off;
127: PetscScalar *A,*B,*W,*X,*U,*Y,*E,*work,*beta,norm;
128: PetscReal *ca,*cb,*cg;
129: PetscBLASInt info,n,ldd,nd,lrwork=0,lwork,one=1;
130: #if defined(PETSC_USE_COMPLEX)
131: PetscReal *rwork;
132: #else
133: PetscScalar norm0;
134: #endif
137: if (!ds->mat[DS_MAT_A]) {
138: DSAllocateMat_Private(ds,DS_MAT_A);
139: }
140: if (!ds->mat[DS_MAT_B]) {
141: DSAllocateMat_Private(ds,DS_MAT_B);
142: }
143: if (!ds->mat[DS_MAT_W]) {
144: DSAllocateMat_Private(ds,DS_MAT_W);
145: }
146: if (!ds->mat[DS_MAT_U]) {
147: DSAllocateMat_Private(ds,DS_MAT_U);
148: }
149: PetscBLASIntCast(ds->n*ctx->d,&nd);
150: PetscBLASIntCast(ds->n,&n);
151: PetscBLASIntCast(ds->ld*ctx->d,&ldd);
152: #if defined(PETSC_USE_COMPLEX)
153: PetscBLASIntCast(nd+2*nd,&lwork);
154: PetscBLASIntCast(8*nd,&lrwork);
155: #else
156: PetscBLASIntCast(nd+8*nd,&lwork);
157: #endif
158: DSAllocateWork_Private(ds,lwork,lrwork,0);
159: beta = ds->work;
160: work = ds->work + nd;
161: lwork -= nd;
162: A = ds->mat[DS_MAT_A];
163: B = ds->mat[DS_MAT_B];
164: W = ds->mat[DS_MAT_W];
165: U = ds->mat[DS_MAT_U];
166: X = ds->mat[DS_MAT_X];
167: Y = ds->mat[DS_MAT_Y];
168: E = ds->mat[DSMatExtra[ctx->d]];
170: /* build matrices A and B of the linearization */
171: PetscMemzero(A,ldd*ldd*sizeof(PetscScalar));
172: if (!ctx->pbc) { /* monomial basis */
173: for (i=0;i<nd-ds->n;i++) A[i+(i+ds->n)*ldd] = 1.0;
174: for (i=0;i<ctx->d;i++) {
175: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
176: for (j=0;j<ds->n;j++) {
177: PetscMemcpy(A+off+j*ldd,ds->mat[DSMatExtra[i]]+j*ds->ld,ds->n*sizeof(PetscScalar));
178: }
179: }
180: } else {
181: ca = ctx->pbc;
182: cb = ca+ctx->d+1;
183: cg = cb+ctx->d+1;
184: for (i=0;i<ds->n;i++) {
185: A[i+(i+ds->n)*ldd] = ca[0];
186: A[i+i*ldd] = cb[0];
187: }
188: for (;i<nd-ds->n;i++) {
189: j = i/ds->n;
190: A[i+(i+ds->n)*ldd] = ca[j];
191: A[i+i*ldd] = cb[j];
192: A[i+(i-ds->n)*ldd] = cg[j];
193: }
194: for (i=0;i<ctx->d-2;i++) {
195: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
196: for (j=0;j<ds->n;j++)
197: for (k=0;k<ds->n;k++)
198: *(A+off+j*ldd+k) = *(ds->mat[DSMatExtra[i]]+j*ds->ld+k)*ca[ctx->d-1];
199: }
200: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
201: for (j=0;j<ds->n;j++)
202: for (k=0;k<ds->n;k++)
203: *(A+off+j*ldd+k) = *(ds->mat[DSMatExtra[i]]+j*ds->ld+k)*ca[ctx->d-1]-E[j*ds->ld+k]*cg[ctx->d-1];
204: off = (++i)*ds->n*ldd+(ctx->d-1)*ds->n;
205: for (j=0;j<ds->n;j++)
206: for (k=0;k<ds->n;k++)
207: *(A+off+j*ldd+k) = *(ds->mat[DSMatExtra[i]]+j*ds->ld+k)*ca[ctx->d-1]-E[j*ds->ld+k]*cb[ctx->d-1];
208: }
209: PetscMemzero(B,ldd*ldd*sizeof(PetscScalar));
210: for (i=0;i<nd-ds->n;i++) B[i+i*ldd] = 1.0;
211: off = (ctx->d-1)*ds->n*(ldd+1);
212: for (j=0;j<ds->n;j++) {
213: for (i=0;i<ds->n;i++) B[off+i+j*ldd] = -E[i+j*ds->ld];
214: }
216: /* solve generalized eigenproblem */
217: #if defined(PETSC_USE_COMPLEX)
218: rwork = ds->rwork;
219: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,rwork,&info));
220: #else
221: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
222: #endif
223: SlepcCheckLapackInfo("ggev",info);
225: /* copy eigenvalues */
226: for (i=0;i<nd;i++) {
227: if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
228: else wr[i] /= beta[i];
229: #if !defined(PETSC_USE_COMPLEX)
230: if (beta[i]==0.0) wi[i] = 0.0;
231: else wi[i] /= beta[i];
232: #else
233: if (wi) wi[i] = 0.0;
234: #endif
235: }
237: /* copy and normalize eigenvectors */
238: for (j=0;j<nd;j++) {
239: PetscMemcpy(X+j*ds->ld,W+j*ldd,ds->n*sizeof(PetscScalar));
240: PetscMemcpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n*sizeof(PetscScalar));
241: }
242: for (j=0;j<nd;j++) {
243: #if !defined(PETSC_USE_COMPLEX)
244: if (wi[j] != 0.0) {
245: norm = BLASnrm2_(&n,X+j*ds->ld,&one);
246: norm0 = BLASnrm2_(&n,X+(j+1)*ds->ld,&one);
247: norm = 1.0/SlepcAbsEigenvalue(norm,norm0);
248: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X+j*ds->ld,&one));
249: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X+(j+1)*ds->ld,&one));
250: norm = BLASnrm2_(&n,Y+j*ds->ld,&one);
251: norm0 = BLASnrm2_(&n,Y+(j+1)*ds->ld,&one);
252: norm = 1.0/SlepcAbsEigenvalue(norm,norm0);
253: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,Y+j*ds->ld,&one));
254: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,Y+(j+1)*ds->ld,&one));
255: j++;
256: } else
257: #endif
258: {
259: norm = 1.0/BLASnrm2_(&n,X+j*ds->ld,&one);
260: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X+j*ds->ld,&one));
261: norm = 1.0/BLASnrm2_(&n,Y+j*ds->ld,&one);
262: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,Y+j*ds->ld,&one));
263: }
264: }
265: return(0);
266: #endif
267: }
269: PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
270: {
272: DS_PEP *ctx = (DS_PEP*)ds->data;
273: PetscInt ld=ds->ld,k=0;
274: PetscMPIInt ldnd,rank,off=0,size,dn;
277: if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
278: if (eigr) k += ctx->d*ds->n;
279: if (eigi) k += ctx->d*ds->n;
280: DSAllocateWork_Private(ds,k,0,0);
281: PetscMPIIntCast(k*sizeof(PetscScalar),&size);
282: PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd);
283: PetscMPIIntCast(ctx->d*ds->n,&dn);
284: MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
285: if (!rank) {
286: if (ds->state>=DS_STATE_CONDENSED) {
287: MPI_Pack(ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
288: MPI_Pack(ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
289: }
290: if (eigr) {
291: MPI_Pack(eigr,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
292: }
293: if (eigi) {
294: MPI_Pack(eigi,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
295: }
296: }
297: MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
298: if (rank) {
299: if (ds->state>=DS_STATE_CONDENSED) {
300: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
301: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
302: }
303: if (eigr) {
304: MPI_Unpack(ds->work,size,&off,eigr,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
305: }
306: if (eigi) {
307: MPI_Unpack(ds->work,size,&off,eigi,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
308: }
309: }
310: return(0);
311: }
313: static PetscErrorCode DSPEPSetDegree_PEP(DS ds,PetscInt d)
314: {
315: DS_PEP *ctx = (DS_PEP*)ds->data;
318: if (d<0) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The degree must be a non-negative integer");
319: if (d>=DS_NUM_EXTRA) SETERRQ1(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Only implemented for polynomials of degree at most %D",DS_NUM_EXTRA-1);
320: ctx->d = d;
321: return(0);
322: }
324: /*@
325: DSPEPSetDegree - Sets the polynomial degree for a DSPEP.
327: Logically Collective on DS
329: Input Parameters:
330: + ds - the direct solver context
331: - d - the degree
333: Level: intermediate
335: .seealso: DSPEPGetDegree()
336: @*/
337: PetscErrorCode DSPEPSetDegree(DS ds,PetscInt d)
338: {
344: PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
345: return(0);
346: }
348: static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
349: {
350: DS_PEP *ctx = (DS_PEP*)ds->data;
353: *d = ctx->d;
354: return(0);
355: }
357: /*@
358: DSPEPGetDegree - Returns the polynomial degree for a DSPEP.
360: Not collective
362: Input Parameter:
363: . ds - the direct solver context
365: Output Parameters:
366: . d - the degree
368: Level: intermediate
370: .seealso: DSPEPSetDegree()
371: @*/
372: PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
373: {
379: PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
380: return(0);
381: }
383: static PetscErrorCode DSPEPSetCoefficients_PEP(DS ds,PetscReal *pbc)
384: {
386: DS_PEP *ctx = (DS_PEP*)ds->data;
387: PetscInt i;
390: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
391: if (ctx->pbc) { PetscFree(ctx->pbc); }
392: PetscMalloc1(3*(ctx->d+1),&ctx->pbc);
393: for (i=0;i<3*(ctx->d+1);i++) ctx->pbc[i] = pbc[i];
394: ds->state = DS_STATE_RAW;
395: return(0);
396: }
398: /*@C
399: DSPEPSetCoefficients - Sets the polynomial basis coefficients for a DSPEP.
401: Logically Collective on DS
403: Input Parameters:
404: + ds - the direct solver context
405: - pbc - the polynomial basis coefficients
407: Notes:
408: This function is required only in the case of a polynomial specified in a
409: non-monomial basis, to provide the coefficients that will be used
410: during the linearization, multiplying the identity blocks on the three main
411: diagonal blocks. Depending on the polynomial basis (Chebyshev, Legendre, ...)
412: the coefficients must be different.
414: There must be a total of 3*(d+1) coefficients, where d is the degree of the
415: polynomial. The coefficients are arranged in three groups: alpha, beta, and
416: gamma, according to the definition of the three-term recurrence. In the case
417: of the monomial basis, alpha=1 and beta=gamma=0, in which case it is not
418: necessary to invoke this function.
420: Level: advanced
422: .seealso: DSPEPGetCoefficients(), DSPEPSetDegree()
423: @*/
424: PetscErrorCode DSPEPSetCoefficients(DS ds,PetscReal *pbc)
425: {
430: PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
431: return(0);
432: }
434: static PetscErrorCode DSPEPGetCoefficients_PEP(DS ds,PetscReal **pbc)
435: {
437: DS_PEP *ctx = (DS_PEP*)ds->data;
438: PetscInt i;
441: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
442: PetscCalloc1(3*(ctx->d+1),pbc);
443: if (ctx->pbc) for (i=0;i<3*(ctx->d+1);i++) (*pbc)[i] = ctx->pbc[i];
444: else for (i=0;i<ctx->d+1;i++) (*pbc)[i] = 1.0;
445: return(0);
446: }
448: /*@C
449: DSPEPGetCoefficients - Returns the polynomial basis coefficients for a DSPEP.
451: Not collective
453: Input Parameter:
454: . ds - the direct solver context
456: Output Parameters:
457: . pbc - the polynomial basis coefficients
459: Note:
460: The returned array has length 3*(d+1) and should be freed by the user.
462: Fortran Note:
463: The calling sequence from Fortran is
464: .vb
465: DSPEPGetCoefficients(eps,pbc,ierr)
466: double precision pbc(d+1) output
467: .ve
469: Level: advanced
471: .seealso: DSPEPSetCoefficients()
472: @*/
473: PetscErrorCode DSPEPGetCoefficients(DS ds,PetscReal **pbc)
474: {
480: PetscUseMethod(ds,"DSPEPGetCoefficients_C",(DS,PetscReal**),(ds,pbc));
481: return(0);
482: }
484: PetscErrorCode DSDestroy_PEP(DS ds)
485: {
487: DS_PEP *ctx = (DS_PEP*)ds->data;
490: if (ctx->pbc) { PetscFree(ctx->pbc); }
491: PetscFree(ds->data);
492: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL);
493: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL);
494: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL);
495: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL);
496: return(0);
497: }
499: PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
500: {
501: DS_PEP *ctx = (DS_PEP*)ds->data;
504: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
505: *rows = ds->n;
506: if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
507: *cols = ds->n;
508: if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U || t==DS_MAT_X || t==DS_MAT_Y) *cols *= ctx->d;
509: return(0);
510: }
512: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
513: {
514: DS_PEP *ctx;
518: PetscNewLog(ds,&ctx);
519: ds->data = (void*)ctx;
521: ds->ops->allocate = DSAllocate_PEP;
522: ds->ops->view = DSView_PEP;
523: ds->ops->vectors = DSVectors_PEP;
524: ds->ops->solve[0] = DSSolve_PEP_QZ;
525: ds->ops->sort = DSSort_PEP;
526: ds->ops->synchronize = DSSynchronize_PEP;
527: ds->ops->destroy = DSDestroy_PEP;
528: ds->ops->matgetsize = DSMatGetSize_PEP;
529: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP);
530: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP);
531: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP);
532: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP);
533: return(0);
534: }