Actual source code: dspep.c
slepc-3.14.2 2021-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/dsimpl.h>
12: #include <slepcblaslapack.h>
14: typedef struct {
15: PetscInt 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) { DSViewMat(ds,viewer,DS_MAT_X); }
53: return(0);
54: }
56: PetscErrorCode DSVectors_PEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
57: {
59: if (rnorm) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
60: switch (mat) {
61: case DS_MAT_X:
62: break;
63: case DS_MAT_Y:
64: break;
65: default:
66: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
67: }
68: return(0);
69: }
71: PetscErrorCode DSSort_PEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *kout)
72: {
74: DS_PEP *ctx = (DS_PEP*)ds->data;
75: PetscInt n,i,j,k,p,*perm,told,ld;
76: PetscScalar *A,*X,*Y,rtmp,rtmp2;
79: if (!ds->sc) return(0);
80: n = ds->n*ctx->d;
81: A = ds->mat[DS_MAT_A];
82: perm = ds->perm;
83: for (i=0;i<n;i++) perm[i] = i;
84: told = ds->t;
85: ds->t = n; /* force the sorting routines to consider d*n eigenvalues */
86: if (rr) {
87: DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
88: } else {
89: DSSortEigenvalues_Private(ds,wr,wi,perm,PETSC_FALSE);
90: }
91: ds->t = told; /* restore value of t */
92: for (i=0;i<n;i++) A[i] = wr[perm[i]];
93: for (i=0;i<n;i++) wr[i] = A[i];
94: for (i=0;i<n;i++) A[i] = wi[perm[i]];
95: for (i=0;i<n;i++) wi[i] = A[i];
96: /* cannot use DSPermuteColumns_Private() since matrix is not square */
97: ld = ds->ld;
98: X = ds->mat[DS_MAT_X];
99: Y = ds->mat[DS_MAT_Y];
100: for (i=0;i<n;i++) {
101: p = perm[i];
102: if (p != i) {
103: j = i + 1;
104: while (perm[j] != i) j++;
105: perm[j] = p; perm[i] = i;
106: /* swap columns i and j */
107: for (k=0;k<ds->n;k++) {
108: rtmp = X[k+p*ld]; X[k+p*ld] = X[k+i*ld]; X[k+i*ld] = rtmp;
109: rtmp2 = Y[k+p*ld]; Y[k+p*ld] = Y[k+i*ld]; Y[k+i*ld] = rtmp2;
110: }
111: }
112: }
113: return(0);
114: }
116: PetscErrorCode DSSolve_PEP_QZ(DS ds,PetscScalar *wr,PetscScalar *wi)
117: {
119: DS_PEP *ctx = (DS_PEP*)ds->data;
120: PetscInt i,j,k,off;
121: PetscScalar *A,*B,*W,*X,*U,*Y,*E,*work,*beta,norm;
122: PetscReal *ca,*cb,*cg;
123: PetscBLASInt info,n,ldd,nd,lrwork=0,lwork,one=1;
124: #if defined(PETSC_USE_COMPLEX)
125: PetscReal *rwork;
126: #else
127: PetscScalar norm0;
128: #endif
131: if (!ds->mat[DS_MAT_A]) {
132: DSAllocateMat_Private(ds,DS_MAT_A);
133: }
134: if (!ds->mat[DS_MAT_B]) {
135: DSAllocateMat_Private(ds,DS_MAT_B);
136: }
137: if (!ds->mat[DS_MAT_W]) {
138: DSAllocateMat_Private(ds,DS_MAT_W);
139: }
140: if (!ds->mat[DS_MAT_U]) {
141: DSAllocateMat_Private(ds,DS_MAT_U);
142: }
143: PetscBLASIntCast(ds->n*ctx->d,&nd);
144: PetscBLASIntCast(ds->n,&n);
145: PetscBLASIntCast(ds->ld*ctx->d,&ldd);
146: #if defined(PETSC_USE_COMPLEX)
147: PetscBLASIntCast(nd+2*nd,&lwork);
148: PetscBLASIntCast(8*nd,&lrwork);
149: #else
150: PetscBLASIntCast(nd+8*nd,&lwork);
151: #endif
152: DSAllocateWork_Private(ds,lwork,lrwork,0);
153: beta = ds->work;
154: work = ds->work + nd;
155: lwork -= nd;
156: A = ds->mat[DS_MAT_A];
157: B = ds->mat[DS_MAT_B];
158: W = ds->mat[DS_MAT_W];
159: U = ds->mat[DS_MAT_U];
160: X = ds->mat[DS_MAT_X];
161: Y = ds->mat[DS_MAT_Y];
162: E = ds->mat[DSMatExtra[ctx->d]];
164: /* build matrices A and B of the linearization */
165: PetscArrayzero(A,ldd*ldd);
166: if (!ctx->pbc) { /* monomial basis */
167: for (i=0;i<nd-ds->n;i++) A[i+(i+ds->n)*ldd] = 1.0;
168: for (i=0;i<ctx->d;i++) {
169: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
170: for (j=0;j<ds->n;j++) {
171: PetscArraycpy(A+off+j*ldd,ds->mat[DSMatExtra[i]]+j*ds->ld,ds->n);
172: }
173: }
174: } else {
175: ca = ctx->pbc;
176: cb = ca+ctx->d+1;
177: cg = cb+ctx->d+1;
178: for (i=0;i<ds->n;i++) {
179: A[i+(i+ds->n)*ldd] = ca[0];
180: A[i+i*ldd] = cb[0];
181: }
182: for (;i<nd-ds->n;i++) {
183: j = i/ds->n;
184: A[i+(i+ds->n)*ldd] = ca[j];
185: A[i+i*ldd] = cb[j];
186: A[i+(i-ds->n)*ldd] = cg[j];
187: }
188: for (i=0;i<ctx->d-2;i++) {
189: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
190: for (j=0;j<ds->n;j++)
191: for (k=0;k<ds->n;k++)
192: *(A+off+j*ldd+k) = *(ds->mat[DSMatExtra[i]]+j*ds->ld+k)*ca[ctx->d-1];
193: }
194: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
195: for (j=0;j<ds->n;j++)
196: for (k=0;k<ds->n;k++)
197: *(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];
198: off = (++i)*ds->n*ldd+(ctx->d-1)*ds->n;
199: for (j=0;j<ds->n;j++)
200: for (k=0;k<ds->n;k++)
201: *(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];
202: }
203: PetscArrayzero(B,ldd*ldd);
204: for (i=0;i<nd-ds->n;i++) B[i+i*ldd] = 1.0;
205: off = (ctx->d-1)*ds->n*(ldd+1);
206: for (j=0;j<ds->n;j++) {
207: for (i=0;i<ds->n;i++) B[off+i+j*ldd] = -E[i+j*ds->ld];
208: }
210: /* solve generalized eigenproblem */
211: #if defined(PETSC_USE_COMPLEX)
212: rwork = ds->rwork;
213: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,rwork,&info));
214: #else
215: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
216: #endif
217: SlepcCheckLapackInfo("ggev",info);
219: /* copy eigenvalues */
220: for (i=0;i<nd;i++) {
221: if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
222: else wr[i] /= beta[i];
223: #if !defined(PETSC_USE_COMPLEX)
224: if (beta[i]==0.0) wi[i] = 0.0;
225: else wi[i] /= beta[i];
226: #else
227: if (wi) wi[i] = 0.0;
228: #endif
229: }
231: /* copy and normalize eigenvectors */
232: for (j=0;j<nd;j++) {
233: PetscArraycpy(X+j*ds->ld,W+j*ldd,ds->n);
234: PetscArraycpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n);
235: }
236: for (j=0;j<nd;j++) {
237: #if !defined(PETSC_USE_COMPLEX)
238: if (wi[j] != 0.0) {
239: norm = BLASnrm2_(&n,X+j*ds->ld,&one);
240: norm0 = BLASnrm2_(&n,X+(j+1)*ds->ld,&one);
241: norm = 1.0/SlepcAbsEigenvalue(norm,norm0);
242: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X+j*ds->ld,&one));
243: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X+(j+1)*ds->ld,&one));
244: norm = BLASnrm2_(&n,Y+j*ds->ld,&one);
245: norm0 = BLASnrm2_(&n,Y+(j+1)*ds->ld,&one);
246: norm = 1.0/SlepcAbsEigenvalue(norm,norm0);
247: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,Y+j*ds->ld,&one));
248: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,Y+(j+1)*ds->ld,&one));
249: j++;
250: } else
251: #endif
252: {
253: norm = 1.0/BLASnrm2_(&n,X+j*ds->ld,&one);
254: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,X+j*ds->ld,&one));
255: norm = 1.0/BLASnrm2_(&n,Y+j*ds->ld,&one);
256: PetscStackCallBLAS("BLASscal",BLASscal_(&n,&norm,Y+j*ds->ld,&one));
257: }
258: }
259: return(0);
260: }
262: PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
263: {
265: DS_PEP *ctx = (DS_PEP*)ds->data;
266: PetscInt ld=ds->ld,k=0;
267: PetscMPIInt ldnd,rank,off=0,size,dn;
270: if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
271: if (eigr) k += ctx->d*ds->n;
272: if (eigi) k += ctx->d*ds->n;
273: DSAllocateWork_Private(ds,k,0,0);
274: PetscMPIIntCast(k*sizeof(PetscScalar),&size);
275: PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd);
276: PetscMPIIntCast(ctx->d*ds->n,&dn);
277: MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
278: if (!rank) {
279: if (ds->state>=DS_STATE_CONDENSED) {
280: MPI_Pack(ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
281: MPI_Pack(ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
282: }
283: if (eigr) {
284: MPI_Pack(eigr,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
285: }
286: if (eigi) {
287: MPI_Pack(eigi,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
288: }
289: }
290: MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
291: if (rank) {
292: if (ds->state>=DS_STATE_CONDENSED) {
293: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
294: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
295: }
296: if (eigr) {
297: MPI_Unpack(ds->work,size,&off,eigr,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
298: }
299: if (eigi) {
300: MPI_Unpack(ds->work,size,&off,eigi,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
301: }
302: }
303: return(0);
304: }
306: static PetscErrorCode DSPEPSetDegree_PEP(DS ds,PetscInt d)
307: {
308: DS_PEP *ctx = (DS_PEP*)ds->data;
311: if (d<0) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The degree must be a non-negative integer");
312: 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);
313: ctx->d = d;
314: return(0);
315: }
317: /*@
318: DSPEPSetDegree - Sets the polynomial degree for a DSPEP.
320: Logically Collective on ds
322: Input Parameters:
323: + ds - the direct solver context
324: - d - the degree
326: Level: intermediate
328: .seealso: DSPEPGetDegree()
329: @*/
330: PetscErrorCode DSPEPSetDegree(DS ds,PetscInt d)
331: {
337: PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
338: return(0);
339: }
341: static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
342: {
343: DS_PEP *ctx = (DS_PEP*)ds->data;
346: *d = ctx->d;
347: return(0);
348: }
350: /*@
351: DSPEPGetDegree - Returns the polynomial degree for a DSPEP.
353: Not collective
355: Input Parameter:
356: . ds - the direct solver context
358: Output Parameters:
359: . d - the degree
361: Level: intermediate
363: .seealso: DSPEPSetDegree()
364: @*/
365: PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
366: {
372: PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
373: return(0);
374: }
376: static PetscErrorCode DSPEPSetCoefficients_PEP(DS ds,PetscReal *pbc)
377: {
379: DS_PEP *ctx = (DS_PEP*)ds->data;
380: PetscInt i;
383: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
384: if (ctx->pbc) { PetscFree(ctx->pbc); }
385: PetscMalloc1(3*(ctx->d+1),&ctx->pbc);
386: for (i=0;i<3*(ctx->d+1);i++) ctx->pbc[i] = pbc[i];
387: ds->state = DS_STATE_RAW;
388: return(0);
389: }
391: /*@C
392: DSPEPSetCoefficients - Sets the polynomial basis coefficients for a DSPEP.
394: Logically Collective on ds
396: Input Parameters:
397: + ds - the direct solver context
398: - pbc - the polynomial basis coefficients
400: Notes:
401: This function is required only in the case of a polynomial specified in a
402: non-monomial basis, to provide the coefficients that will be used
403: during the linearization, multiplying the identity blocks on the three main
404: diagonal blocks. Depending on the polynomial basis (Chebyshev, Legendre, ...)
405: the coefficients must be different.
407: There must be a total of 3*(d+1) coefficients, where d is the degree of the
408: polynomial. The coefficients are arranged in three groups: alpha, beta, and
409: gamma, according to the definition of the three-term recurrence. In the case
410: of the monomial basis, alpha=1 and beta=gamma=0, in which case it is not
411: necessary to invoke this function.
413: Level: advanced
415: .seealso: DSPEPGetCoefficients(), DSPEPSetDegree()
416: @*/
417: PetscErrorCode DSPEPSetCoefficients(DS ds,PetscReal *pbc)
418: {
423: PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
424: return(0);
425: }
427: static PetscErrorCode DSPEPGetCoefficients_PEP(DS ds,PetscReal **pbc)
428: {
430: DS_PEP *ctx = (DS_PEP*)ds->data;
431: PetscInt i;
434: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
435: PetscCalloc1(3*(ctx->d+1),pbc);
436: if (ctx->pbc) for (i=0;i<3*(ctx->d+1);i++) (*pbc)[i] = ctx->pbc[i];
437: else for (i=0;i<ctx->d+1;i++) (*pbc)[i] = 1.0;
438: return(0);
439: }
441: /*@C
442: DSPEPGetCoefficients - Returns the polynomial basis coefficients for a DSPEP.
444: Not collective
446: Input Parameter:
447: . ds - the direct solver context
449: Output Parameters:
450: . pbc - the polynomial basis coefficients
452: Note:
453: The returned array has length 3*(d+1) and should be freed by the user.
455: Fortran Note:
456: The calling sequence from Fortran is
457: .vb
458: DSPEPGetCoefficients(eps,pbc,ierr)
459: double precision pbc(d+1) output
460: .ve
462: Level: advanced
464: .seealso: DSPEPSetCoefficients()
465: @*/
466: PetscErrorCode DSPEPGetCoefficients(DS ds,PetscReal **pbc)
467: {
473: PetscUseMethod(ds,"DSPEPGetCoefficients_C",(DS,PetscReal**),(ds,pbc));
474: return(0);
475: }
477: PetscErrorCode DSDestroy_PEP(DS ds)
478: {
480: DS_PEP *ctx = (DS_PEP*)ds->data;
483: if (ctx->pbc) { PetscFree(ctx->pbc); }
484: PetscFree(ds->data);
485: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL);
486: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL);
487: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL);
488: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL);
489: return(0);
490: }
492: PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
493: {
494: DS_PEP *ctx = (DS_PEP*)ds->data;
497: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
498: *rows = ds->n;
499: if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
500: *cols = ds->n;
501: 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;
502: return(0);
503: }
505: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
506: {
507: DS_PEP *ctx;
511: PetscNewLog(ds,&ctx);
512: ds->data = (void*)ctx;
514: ds->ops->allocate = DSAllocate_PEP;
515: ds->ops->view = DSView_PEP;
516: ds->ops->vectors = DSVectors_PEP;
517: ds->ops->solve[0] = DSSolve_PEP_QZ;
518: ds->ops->sort = DSSort_PEP;
519: ds->ops->synchronize = DSSynchronize_PEP;
520: ds->ops->destroy = DSDestroy_PEP;
521: ds->ops->matgetsize = DSMatGetSize_PEP;
522: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP);
523: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP);
524: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP);
525: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP);
526: return(0);
527: }