Actual source code: dspep.c
slepc-3.16.0 2021-09-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/dsimpl.h>
12: #include <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: if (format == PETSC_VIEWER_ASCII_INFO) return(0);
48: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
49: PetscViewerASCIIPrintf(viewer,"polynomial degree: %D\n",ctx->d);
50: return(0);
51: }
52: for (i=0;i<=ctx->d;i++) {
53: DSViewMat(ds,viewer,DSMatExtra[i]);
54: }
55: if (ds->state>DS_STATE_INTERMEDIATE) { DSViewMat(ds,viewer,DS_MAT_X); }
56: return(0);
57: }
59: PetscErrorCode DSVectors_PEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
60: {
62: if (rnorm) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
63: switch (mat) {
64: case DS_MAT_X:
65: break;
66: case DS_MAT_Y:
67: break;
68: default:
69: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
70: }
71: return(0);
72: }
74: PetscErrorCode DSSort_PEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *kout)
75: {
77: DS_PEP *ctx = (DS_PEP*)ds->data;
78: PetscInt n,i,*perm,told;
79: PetscScalar *A;
82: if (!ds->sc) return(0);
83: n = ds->n*ctx->d;
84: A = ds->mat[DS_MAT_A];
85: perm = ds->perm;
86: for (i=0;i<n;i++) perm[i] = i;
87: told = ds->t;
88: ds->t = n; /* force the sorting routines to consider d*n eigenvalues */
89: if (rr) {
90: DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
91: } else {
92: DSSortEigenvalues_Private(ds,wr,wi,perm,PETSC_FALSE);
93: }
94: ds->t = told; /* restore value of t */
95: for (i=0;i<n;i++) A[i] = wr[perm[i]];
96: for (i=0;i<n;i++) wr[i] = A[i];
97: for (i=0;i<n;i++) A[i] = wi[perm[i]];
98: for (i=0;i<n;i++) wi[i] = A[i];
99: DSPermuteColumnsTwo_Private(ds,0,n,ds->n,DS_MAT_X,DS_MAT_Y,perm);
100: return(0);
101: }
103: PetscErrorCode DSSolve_PEP_QZ(DS ds,PetscScalar *wr,PetscScalar *wi)
104: {
106: DS_PEP *ctx = (DS_PEP*)ds->data;
107: PetscInt i,j,k,off;
108: PetscScalar *A,*B,*W,*X,*U,*Y,*E,*work,*beta;
109: PetscReal *ca,*cb,*cg,norm,done=1.0;
110: PetscBLASInt info,n,ld,ldd,nd,lrwork=0,lwork,one=1,zero=0,cols;
111: #if defined(PETSC_USE_COMPLEX)
112: PetscReal *rwork;
113: #endif
116: if (!ds->mat[DS_MAT_A]) {
117: DSAllocateMat_Private(ds,DS_MAT_A);
118: }
119: if (!ds->mat[DS_MAT_B]) {
120: DSAllocateMat_Private(ds,DS_MAT_B);
121: }
122: if (!ds->mat[DS_MAT_W]) {
123: DSAllocateMat_Private(ds,DS_MAT_W);
124: }
125: if (!ds->mat[DS_MAT_U]) {
126: DSAllocateMat_Private(ds,DS_MAT_U);
127: }
128: PetscBLASIntCast(ds->n*ctx->d,&nd);
129: PetscBLASIntCast(ds->n,&n);
130: PetscBLASIntCast(ds->ld,&ld);
131: PetscBLASIntCast(ds->ld*ctx->d,&ldd);
132: #if defined(PETSC_USE_COMPLEX)
133: PetscBLASIntCast(nd+2*nd,&lwork);
134: PetscBLASIntCast(8*nd,&lrwork);
135: #else
136: PetscBLASIntCast(nd+8*nd,&lwork);
137: #endif
138: DSAllocateWork_Private(ds,lwork,lrwork,0);
139: beta = ds->work;
140: work = ds->work + nd;
141: lwork -= nd;
142: A = ds->mat[DS_MAT_A];
143: B = ds->mat[DS_MAT_B];
144: W = ds->mat[DS_MAT_W];
145: U = ds->mat[DS_MAT_U];
146: X = ds->mat[DS_MAT_X];
147: Y = ds->mat[DS_MAT_Y];
148: E = ds->mat[DSMatExtra[ctx->d]];
150: /* build matrices A and B of the linearization */
151: PetscArrayzero(A,ldd*ldd);
152: if (!ctx->pbc) { /* monomial basis */
153: for (i=0;i<nd-ds->n;i++) A[i+(i+ds->n)*ldd] = 1.0;
154: for (i=0;i<ctx->d;i++) {
155: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
156: for (j=0;j<ds->n;j++) {
157: PetscArraycpy(A+off+j*ldd,ds->mat[DSMatExtra[i]]+j*ds->ld,ds->n);
158: }
159: }
160: } else {
161: ca = ctx->pbc;
162: cb = ca+ctx->d+1;
163: cg = cb+ctx->d+1;
164: for (i=0;i<ds->n;i++) {
165: A[i+(i+ds->n)*ldd] = ca[0];
166: A[i+i*ldd] = cb[0];
167: }
168: for (;i<nd-ds->n;i++) {
169: j = i/ds->n;
170: A[i+(i+ds->n)*ldd] = ca[j];
171: A[i+i*ldd] = cb[j];
172: A[i+(i-ds->n)*ldd] = cg[j];
173: }
174: for (i=0;i<ctx->d-2;i++) {
175: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
176: for (j=0;j<ds->n;j++)
177: for (k=0;k<ds->n;k++)
178: *(A+off+j*ldd+k) = *(ds->mat[DSMatExtra[i]]+j*ds->ld+k)*ca[ctx->d-1];
179: }
180: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
181: for (j=0;j<ds->n;j++)
182: for (k=0;k<ds->n;k++)
183: *(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];
184: off = (++i)*ds->n*ldd+(ctx->d-1)*ds->n;
185: for (j=0;j<ds->n;j++)
186: for (k=0;k<ds->n;k++)
187: *(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];
188: }
189: PetscArrayzero(B,ldd*ldd);
190: for (i=0;i<nd-ds->n;i++) B[i+i*ldd] = 1.0;
191: off = (ctx->d-1)*ds->n*(ldd+1);
192: for (j=0;j<ds->n;j++) {
193: for (i=0;i<ds->n;i++) B[off+i+j*ldd] = -E[i+j*ds->ld];
194: }
196: /* solve generalized eigenproblem */
197: #if defined(PETSC_USE_COMPLEX)
198: rwork = ds->rwork;
199: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,rwork,&info));
200: #else
201: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
202: #endif
203: SlepcCheckLapackInfo("ggev",info);
205: /* copy eigenvalues */
206: for (i=0;i<nd;i++) {
207: if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
208: else wr[i] /= beta[i];
209: #if !defined(PETSC_USE_COMPLEX)
210: if (beta[i]==0.0) wi[i] = 0.0;
211: else wi[i] /= beta[i];
212: #else
213: if (wi) wi[i] = 0.0;
214: #endif
215: }
217: /* copy and normalize eigenvectors */
218: for (j=0;j<nd;j++) {
219: PetscArraycpy(X+j*ds->ld,W+j*ldd,ds->n);
220: PetscArraycpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n);
221: }
222: for (j=0;j<nd;j++) {
223: cols = 1;
224: norm = BLASnrm2_(&n,X+j*ds->ld,&one);
225: #if !defined(PETSC_USE_COMPLEX)
226: if (wi[j] != 0.0) {
227: norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,X+(j+1)*ds->ld,&one));
228: cols = 2;
229: }
230: #endif
231: PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,X+j*ds->ld,&ld,&info));
232: SlepcCheckLapackInfo("lascl",info);
233: norm = BLASnrm2_(&n,Y+j*ds->ld,&one);
234: #if !defined(PETSC_USE_COMPLEX)
235: if (wi[j] != 0.0) norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,Y+(j+1)*ds->ld,&one));
236: #endif
237: PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Y+j*ds->ld,&ld,&info));
238: SlepcCheckLapackInfo("lascl",info);
239: #if !defined(PETSC_USE_COMPLEX)
240: if (wi[j] != 0.0) j++;
241: #endif
242: }
243: return(0);
244: }
246: PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
247: {
249: DS_PEP *ctx = (DS_PEP*)ds->data;
250: PetscInt ld=ds->ld,k=0;
251: PetscMPIInt ldnd,rank,off=0,size,dn;
254: if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
255: if (eigr) k += ctx->d*ds->n;
256: if (eigi) k += ctx->d*ds->n;
257: DSAllocateWork_Private(ds,k,0,0);
258: PetscMPIIntCast(k*sizeof(PetscScalar),&size);
259: PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd);
260: PetscMPIIntCast(ctx->d*ds->n,&dn);
261: MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
262: if (!rank) {
263: if (ds->state>=DS_STATE_CONDENSED) {
264: MPI_Pack(ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
265: MPI_Pack(ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
266: }
267: if (eigr) {
268: MPI_Pack(eigr,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
269: }
270: #if !defined(PETSC_USE_COMPLEX)
271: if (eigi) {
272: MPI_Pack(eigi,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
273: }
274: #endif
275: }
276: MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
277: if (rank) {
278: if (ds->state>=DS_STATE_CONDENSED) {
279: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
280: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
281: }
282: if (eigr) {
283: MPI_Unpack(ds->work,size,&off,eigr,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
284: }
285: #if !defined(PETSC_USE_COMPLEX)
286: if (eigi) {
287: MPI_Unpack(ds->work,size,&off,eigi,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
288: }
289: #endif
290: }
291: return(0);
292: }
294: static PetscErrorCode DSPEPSetDegree_PEP(DS ds,PetscInt d)
295: {
296: DS_PEP *ctx = (DS_PEP*)ds->data;
299: if (d<0) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The degree must be a non-negative integer");
300: 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);
301: ctx->d = d;
302: return(0);
303: }
305: /*@
306: DSPEPSetDegree - Sets the polynomial degree for a DSPEP.
308: Logically Collective on ds
310: Input Parameters:
311: + ds - the direct solver context
312: - d - the degree
314: Level: intermediate
316: .seealso: DSPEPGetDegree()
317: @*/
318: PetscErrorCode DSPEPSetDegree(DS ds,PetscInt d)
319: {
325: PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
326: return(0);
327: }
329: static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
330: {
331: DS_PEP *ctx = (DS_PEP*)ds->data;
334: *d = ctx->d;
335: return(0);
336: }
338: /*@
339: DSPEPGetDegree - Returns the polynomial degree for a DSPEP.
341: Not collective
343: Input Parameter:
344: . ds - the direct solver context
346: Output Parameters:
347: . d - the degree
349: Level: intermediate
351: .seealso: DSPEPSetDegree()
352: @*/
353: PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
354: {
360: PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
361: return(0);
362: }
364: static PetscErrorCode DSPEPSetCoefficients_PEP(DS ds,PetscReal *pbc)
365: {
367: DS_PEP *ctx = (DS_PEP*)ds->data;
368: PetscInt i;
371: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
372: if (ctx->pbc) { PetscFree(ctx->pbc); }
373: PetscMalloc1(3*(ctx->d+1),&ctx->pbc);
374: for (i=0;i<3*(ctx->d+1);i++) ctx->pbc[i] = pbc[i];
375: ds->state = DS_STATE_RAW;
376: return(0);
377: }
379: /*@C
380: DSPEPSetCoefficients - Sets the polynomial basis coefficients for a DSPEP.
382: Logically Collective on ds
384: Input Parameters:
385: + ds - the direct solver context
386: - pbc - the polynomial basis coefficients
388: Notes:
389: This function is required only in the case of a polynomial specified in a
390: non-monomial basis, to provide the coefficients that will be used
391: during the linearization, multiplying the identity blocks on the three main
392: diagonal blocks. Depending on the polynomial basis (Chebyshev, Legendre, ...)
393: the coefficients must be different.
395: There must be a total of 3*(d+1) coefficients, where d is the degree of the
396: polynomial. The coefficients are arranged in three groups: alpha, beta, and
397: gamma, according to the definition of the three-term recurrence. In the case
398: of the monomial basis, alpha=1 and beta=gamma=0, in which case it is not
399: necessary to invoke this function.
401: Level: advanced
403: .seealso: DSPEPGetCoefficients(), DSPEPSetDegree()
404: @*/
405: PetscErrorCode DSPEPSetCoefficients(DS ds,PetscReal *pbc)
406: {
411: PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
412: return(0);
413: }
415: static PetscErrorCode DSPEPGetCoefficients_PEP(DS ds,PetscReal **pbc)
416: {
418: DS_PEP *ctx = (DS_PEP*)ds->data;
419: PetscInt i;
422: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
423: PetscCalloc1(3*(ctx->d+1),pbc);
424: if (ctx->pbc) for (i=0;i<3*(ctx->d+1);i++) (*pbc)[i] = ctx->pbc[i];
425: else for (i=0;i<ctx->d+1;i++) (*pbc)[i] = 1.0;
426: return(0);
427: }
429: /*@C
430: DSPEPGetCoefficients - Returns the polynomial basis coefficients for a DSPEP.
432: Not collective
434: Input Parameter:
435: . ds - the direct solver context
437: Output Parameters:
438: . pbc - the polynomial basis coefficients
440: Note:
441: The returned array has length 3*(d+1) and should be freed by the user.
443: Fortran Note:
444: The calling sequence from Fortran is
445: .vb
446: DSPEPGetCoefficients(eps,pbc,ierr)
447: double precision pbc(d+1) output
448: .ve
450: Level: advanced
452: .seealso: DSPEPSetCoefficients()
453: @*/
454: PetscErrorCode DSPEPGetCoefficients(DS ds,PetscReal **pbc)
455: {
461: PetscUseMethod(ds,"DSPEPGetCoefficients_C",(DS,PetscReal**),(ds,pbc));
462: return(0);
463: }
465: PetscErrorCode DSDestroy_PEP(DS ds)
466: {
468: DS_PEP *ctx = (DS_PEP*)ds->data;
471: if (ctx->pbc) { PetscFree(ctx->pbc); }
472: PetscFree(ds->data);
473: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL);
474: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL);
475: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL);
476: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL);
477: return(0);
478: }
480: PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
481: {
482: DS_PEP *ctx = (DS_PEP*)ds->data;
485: if (!ctx->d) SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
486: *rows = ds->n;
487: if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
488: *cols = ds->n;
489: 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;
490: return(0);
491: }
493: /*MC
494: DSPEP - Dense Polynomial Eigenvalue Problem.
496: Level: beginner
498: Notes:
499: The problem is expressed as P(lambda)*x = 0, where P(.) is a matrix
500: polynomial of degree d. The eigenvalues lambda are the arguments
501: returned by DSSolve().
503: The degree of the polynomial, d, can be set with DSPEPSetDegree(), with
504: the first d+1 extra matrices of the DS defining the matrix polynomial. By
505: default, the polynomial is expressed in the monomial basis, but a
506: different basis can be used by setting the corresponding coefficients
507: via DSPEPSetCoefficients().
509: The problem is solved via linearization, by building a pencil (A,B) of
510: size p*n and solving the corresponding GNHEP.
512: Used DS matrices:
513: + DS_MAT_Ex - coefficients of the matrix polynomial
514: . DS_MAT_A - (workspace) first matrix of the linearization
515: . DS_MAT_B - (workspace) second matrix of the linearization
516: . DS_MAT_W - (workspace) right eigenvectors of the linearization
517: - DS_MAT_U - (workspace) left eigenvectors of the linearization
519: Implemented methods:
520: . 0 - QZ iteration on the linearization (_ggev)
522: .seealso: DSCreate(), DSSetType(), DSType, DSPEPSetDegree(), DSPEPSetCoefficients()
523: M*/
524: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
525: {
526: DS_PEP *ctx;
530: PetscNewLog(ds,&ctx);
531: ds->data = (void*)ctx;
533: ds->ops->allocate = DSAllocate_PEP;
534: ds->ops->view = DSView_PEP;
535: ds->ops->vectors = DSVectors_PEP;
536: ds->ops->solve[0] = DSSolve_PEP_QZ;
537: ds->ops->sort = DSSort_PEP;
538: ds->ops->synchronize = DSSynchronize_PEP;
539: ds->ops->destroy = DSDestroy_PEP;
540: ds->ops->matgetsize = DSMatGetSize_PEP;
541: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP);
542: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP);
543: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP);
544: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP);
545: return(0);
546: }