Actual source code: dspep.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 <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: {
21: DS_PEP *ctx = (DS_PEP*)ds->data;
22: PetscInt i;
25: DSAllocateMat_Private(ds,DS_MAT_X);
26: DSAllocateMat_Private(ds,DS_MAT_Y);
27: for (i=0;i<=ctx->d;i++) DSAllocateMat_Private(ds,DSMatExtra[i]);
28: PetscFree(ds->perm);
29: PetscMalloc1(ld*ctx->d,&ds->perm);
30: PetscLogObjectMemory((PetscObject)ds,ld*ctx->d*sizeof(PetscInt));
31: PetscFunctionReturn(0);
32: }
34: PetscErrorCode DSView_PEP(DS ds,PetscViewer viewer)
35: {
36: DS_PEP *ctx = (DS_PEP*)ds->data;
37: PetscViewerFormat format;
38: PetscInt i;
40: PetscViewerGetFormat(viewer,&format);
41: if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(0);
42: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
43: PetscViewerASCIIPrintf(viewer,"polynomial degree: %" PetscInt_FMT "\n",ctx->d);
44: PetscFunctionReturn(0);
45: }
46: for (i=0;i<=ctx->d;i++) DSViewMat(ds,viewer,DSMatExtra[i]);
47: if (ds->state>DS_STATE_INTERMEDIATE) DSViewMat(ds,viewer,DS_MAT_X);
48: PetscFunctionReturn(0);
49: }
51: PetscErrorCode DSVectors_PEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
52: {
54: switch (mat) {
55: case DS_MAT_X:
56: break;
57: case DS_MAT_Y:
58: break;
59: default:
60: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
61: }
62: PetscFunctionReturn(0);
63: }
65: PetscErrorCode DSSort_PEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *kout)
66: {
67: DS_PEP *ctx = (DS_PEP*)ds->data;
68: PetscInt n,i,*perm,told;
69: PetscScalar *A;
71: if (!ds->sc) PetscFunctionReturn(0);
72: n = ds->n*ctx->d;
73: A = ds->mat[DS_MAT_A];
74: perm = ds->perm;
75: for (i=0;i<n;i++) perm[i] = i;
76: told = ds->t;
77: ds->t = n; /* force the sorting routines to consider d*n eigenvalues */
78: if (rr) DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE);
79: else DSSortEigenvalues_Private(ds,wr,wi,perm,PETSC_FALSE);
80: ds->t = told; /* restore value of t */
81: for (i=0;i<n;i++) A[i] = wr[perm[i]];
82: for (i=0;i<n;i++) wr[i] = A[i];
83: for (i=0;i<n;i++) A[i] = wi[perm[i]];
84: for (i=0;i<n;i++) wi[i] = A[i];
85: DSPermuteColumnsTwo_Private(ds,0,n,ds->n,DS_MAT_X,DS_MAT_Y,perm);
86: PetscFunctionReturn(0);
87: }
89: PetscErrorCode DSSolve_PEP_QZ(DS ds,PetscScalar *wr,PetscScalar *wi)
90: {
91: DS_PEP *ctx = (DS_PEP*)ds->data;
92: PetscInt i,j,k,off;
93: PetscScalar *A,*B,*W,*X,*U,*Y,*E,*work,*beta;
94: PetscReal *ca,*cb,*cg,norm,done=1.0;
95: PetscBLASInt info,n,ld,ldd,nd,lrwork=0,lwork,one=1,zero=0,cols;
96: #if defined(PETSC_USE_COMPLEX)
97: PetscReal *rwork;
98: #endif
100: if (!ds->mat[DS_MAT_A]) DSAllocateMat_Private(ds,DS_MAT_A);
101: if (!ds->mat[DS_MAT_B]) DSAllocateMat_Private(ds,DS_MAT_B);
102: if (!ds->mat[DS_MAT_W]) DSAllocateMat_Private(ds,DS_MAT_W);
103: if (!ds->mat[DS_MAT_U]) DSAllocateMat_Private(ds,DS_MAT_U);
104: PetscBLASIntCast(ds->n*ctx->d,&nd);
105: PetscBLASIntCast(ds->n,&n);
106: PetscBLASIntCast(ds->ld,&ld);
107: PetscBLASIntCast(ds->ld*ctx->d,&ldd);
108: #if defined(PETSC_USE_COMPLEX)
109: PetscBLASIntCast(nd+2*nd,&lwork);
110: PetscBLASIntCast(8*nd,&lrwork);
111: #else
112: PetscBLASIntCast(nd+8*nd,&lwork);
113: #endif
114: DSAllocateWork_Private(ds,lwork,lrwork,0);
115: beta = ds->work;
116: work = ds->work + nd;
117: lwork -= nd;
118: A = ds->mat[DS_MAT_A];
119: B = ds->mat[DS_MAT_B];
120: W = ds->mat[DS_MAT_W];
121: U = ds->mat[DS_MAT_U];
122: X = ds->mat[DS_MAT_X];
123: Y = ds->mat[DS_MAT_Y];
124: E = ds->mat[DSMatExtra[ctx->d]];
126: /* build matrices A and B of the linearization */
127: PetscArrayzero(A,ldd*ldd);
128: if (!ctx->pbc) { /* monomial basis */
129: for (i=0;i<nd-ds->n;i++) A[i+(i+ds->n)*ldd] = 1.0;
130: for (i=0;i<ctx->d;i++) {
131: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
132: for (j=0;j<ds->n;j++) PetscArraycpy(A+off+j*ldd,ds->mat[DSMatExtra[i]]+j*ds->ld,ds->n);
133: }
134: } else {
135: ca = ctx->pbc;
136: cb = ca+ctx->d+1;
137: cg = cb+ctx->d+1;
138: for (i=0;i<ds->n;i++) {
139: A[i+(i+ds->n)*ldd] = ca[0];
140: A[i+i*ldd] = cb[0];
141: }
142: for (;i<nd-ds->n;i++) {
143: j = i/ds->n;
144: A[i+(i+ds->n)*ldd] = ca[j];
145: A[i+i*ldd] = cb[j];
146: A[i+(i-ds->n)*ldd] = cg[j];
147: }
148: for (i=0;i<ctx->d-2;i++) {
149: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
150: for (j=0;j<ds->n;j++)
151: for (k=0;k<ds->n;k++)
152: *(A+off+j*ldd+k) = *(ds->mat[DSMatExtra[i]]+j*ds->ld+k)*ca[ctx->d-1];
153: }
154: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
155: for (j=0;j<ds->n;j++)
156: for (k=0;k<ds->n;k++)
157: *(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];
158: off = (++i)*ds->n*ldd+(ctx->d-1)*ds->n;
159: for (j=0;j<ds->n;j++)
160: for (k=0;k<ds->n;k++)
161: *(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];
162: }
163: PetscArrayzero(B,ldd*ldd);
164: for (i=0;i<nd-ds->n;i++) B[i+i*ldd] = 1.0;
165: off = (ctx->d-1)*ds->n*(ldd+1);
166: for (j=0;j<ds->n;j++) {
167: for (i=0;i<ds->n;i++) B[off+i+j*ldd] = -E[i+j*ds->ld];
168: }
170: /* solve generalized eigenproblem */
171: #if defined(PETSC_USE_COMPLEX)
172: rwork = ds->rwork;
173: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,rwork,&info));
174: #else
175: PetscStackCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
176: #endif
177: SlepcCheckLapackInfo("ggev",info);
179: /* copy eigenvalues */
180: for (i=0;i<nd;i++) {
181: if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
182: else wr[i] /= beta[i];
183: #if !defined(PETSC_USE_COMPLEX)
184: if (beta[i]==0.0) wi[i] = 0.0;
185: else wi[i] /= beta[i];
186: #else
187: if (wi) wi[i] = 0.0;
188: #endif
189: }
191: /* copy and normalize eigenvectors */
192: for (j=0;j<nd;j++) {
193: PetscArraycpy(X+j*ds->ld,W+j*ldd,ds->n);
194: PetscArraycpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n);
195: }
196: for (j=0;j<nd;j++) {
197: cols = 1;
198: norm = BLASnrm2_(&n,X+j*ds->ld,&one);
199: #if !defined(PETSC_USE_COMPLEX)
200: if (wi[j] != 0.0) {
201: norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,X+(j+1)*ds->ld,&one));
202: cols = 2;
203: }
204: #endif
205: PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,X+j*ds->ld,&ld,&info));
206: SlepcCheckLapackInfo("lascl",info);
207: norm = BLASnrm2_(&n,Y+j*ds->ld,&one);
208: #if !defined(PETSC_USE_COMPLEX)
209: if (wi[j] != 0.0) norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,Y+(j+1)*ds->ld,&one));
210: #endif
211: PetscStackCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Y+j*ds->ld,&ld,&info));
212: SlepcCheckLapackInfo("lascl",info);
213: #if !defined(PETSC_USE_COMPLEX)
214: if (wi[j] != 0.0) j++;
215: #endif
216: }
217: PetscFunctionReturn(0);
218: }
220: PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
221: {
222: DS_PEP *ctx = (DS_PEP*)ds->data;
223: PetscInt ld=ds->ld,k=0;
224: PetscMPIInt ldnd,rank,off=0,size,dn;
226: if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
227: if (eigr) k += ctx->d*ds->n;
228: if (eigi) k += ctx->d*ds->n;
229: DSAllocateWork_Private(ds,k,0,0);
230: PetscMPIIntCast(k*sizeof(PetscScalar),&size);
231: PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd);
232: PetscMPIIntCast(ctx->d*ds->n,&dn);
233: MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank);
234: if (!rank) {
235: if (ds->state>=DS_STATE_CONDENSED) {
236: MPI_Pack(ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
237: MPI_Pack(ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
238: }
239: if (eigr) MPI_Pack(eigr,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
240: #if !defined(PETSC_USE_COMPLEX)
241: if (eigi) MPI_Pack(eigi,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds));
242: #endif
243: }
244: MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds));
245: if (rank) {
246: if (ds->state>=DS_STATE_CONDENSED) {
247: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_X],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
248: MPI_Unpack(ds->work,size,&off,ds->mat[DS_MAT_Y],ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
249: }
250: if (eigr) MPI_Unpack(ds->work,size,&off,eigr,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
251: #if !defined(PETSC_USE_COMPLEX)
252: if (eigi) MPI_Unpack(ds->work,size,&off,eigi,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds));
253: #endif
254: }
255: PetscFunctionReturn(0);
256: }
258: static PetscErrorCode DSPEPSetDegree_PEP(DS ds,PetscInt d)
259: {
260: DS_PEP *ctx = (DS_PEP*)ds->data;
264: ctx->d = d;
265: PetscFunctionReturn(0);
266: }
268: /*@
269: DSPEPSetDegree - Sets the polynomial degree for a DSPEP.
271: Logically Collective on ds
273: Input Parameters:
274: + ds - the direct solver context
275: - d - the degree
277: Level: intermediate
279: .seealso: DSPEPGetDegree()
280: @*/
281: PetscErrorCode DSPEPSetDegree(DS ds,PetscInt d)
282: {
285: PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
286: PetscFunctionReturn(0);
287: }
289: static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
290: {
291: DS_PEP *ctx = (DS_PEP*)ds->data;
293: *d = ctx->d;
294: PetscFunctionReturn(0);
295: }
297: /*@
298: DSPEPGetDegree - Returns the polynomial degree for a DSPEP.
300: Not collective
302: Input Parameter:
303: . ds - the direct solver context
305: Output Parameters:
306: . d - the degree
308: Level: intermediate
310: .seealso: DSPEPSetDegree()
311: @*/
312: PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
313: {
316: PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
317: PetscFunctionReturn(0);
318: }
320: static PetscErrorCode DSPEPSetCoefficients_PEP(DS ds,PetscReal *pbc)
321: {
322: DS_PEP *ctx = (DS_PEP*)ds->data;
323: PetscInt i;
326: if (ctx->pbc) PetscFree(ctx->pbc);
327: PetscMalloc1(3*(ctx->d+1),&ctx->pbc);
328: for (i=0;i<3*(ctx->d+1);i++) ctx->pbc[i] = pbc[i];
329: ds->state = DS_STATE_RAW;
330: PetscFunctionReturn(0);
331: }
333: /*@C
334: DSPEPSetCoefficients - Sets the polynomial basis coefficients for a DSPEP.
336: Logically Collective on ds
338: Input Parameters:
339: + ds - the direct solver context
340: - pbc - the polynomial basis coefficients
342: Notes:
343: This function is required only in the case of a polynomial specified in a
344: non-monomial basis, to provide the coefficients that will be used
345: during the linearization, multiplying the identity blocks on the three main
346: diagonal blocks. Depending on the polynomial basis (Chebyshev, Legendre, ...)
347: the coefficients must be different.
349: There must be a total of 3*(d+1) coefficients, where d is the degree of the
350: polynomial. The coefficients are arranged in three groups, alpha, beta, and
351: gamma, according to the definition of the three-term recurrence. In the case
352: of the monomial basis, alpha=1 and beta=gamma=0, in which case it is not
353: necessary to invoke this function.
355: Level: advanced
357: .seealso: DSPEPGetCoefficients(), DSPEPSetDegree()
358: @*/
359: PetscErrorCode DSPEPSetCoefficients(DS ds,PetscReal *pbc)
360: {
362: PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
363: PetscFunctionReturn(0);
364: }
366: static PetscErrorCode DSPEPGetCoefficients_PEP(DS ds,PetscReal **pbc)
367: {
368: DS_PEP *ctx = (DS_PEP*)ds->data;
369: PetscInt i;
372: PetscCalloc1(3*(ctx->d+1),pbc);
373: if (ctx->pbc) for (i=0;i<3*(ctx->d+1);i++) (*pbc)[i] = ctx->pbc[i];
374: else for (i=0;i<ctx->d+1;i++) (*pbc)[i] = 1.0;
375: PetscFunctionReturn(0);
376: }
378: /*@C
379: DSPEPGetCoefficients - Returns the polynomial basis coefficients for a DSPEP.
381: Not collective
383: Input Parameter:
384: . ds - the direct solver context
386: Output Parameters:
387: . pbc - the polynomial basis coefficients
389: Note:
390: The returned array has length 3*(d+1) and should be freed by the user.
392: Fortran Notes:
393: The calling sequence from Fortran is
394: .vb
395: DSPEPGetCoefficients(eps,pbc,ierr)
396: double precision pbc(d+1) output
397: .ve
399: Level: advanced
401: .seealso: DSPEPSetCoefficients()
402: @*/
403: PetscErrorCode DSPEPGetCoefficients(DS ds,PetscReal **pbc)
404: {
407: PetscUseMethod(ds,"DSPEPGetCoefficients_C",(DS,PetscReal**),(ds,pbc));
408: PetscFunctionReturn(0);
409: }
411: PetscErrorCode DSDestroy_PEP(DS ds)
412: {
413: DS_PEP *ctx = (DS_PEP*)ds->data;
415: if (ctx->pbc) PetscFree(ctx->pbc);
416: PetscFree(ds->data);
417: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL);
418: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL);
419: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL);
420: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL);
421: PetscFunctionReturn(0);
422: }
424: PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
425: {
426: DS_PEP *ctx = (DS_PEP*)ds->data;
429: *rows = ds->n;
430: if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
431: *cols = ds->n;
432: 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;
433: PetscFunctionReturn(0);
434: }
436: /*MC
437: DSPEP - Dense Polynomial Eigenvalue Problem.
439: Level: beginner
441: Notes:
442: The problem is expressed as P(lambda)*x = 0, where P(.) is a matrix
443: polynomial of degree d. The eigenvalues lambda are the arguments
444: returned by DSSolve().
446: The degree of the polynomial, d, can be set with DSPEPSetDegree(), with
447: the first d+1 extra matrices of the DS defining the matrix polynomial. By
448: default, the polynomial is expressed in the monomial basis, but a
449: different basis can be used by setting the corresponding coefficients
450: via DSPEPSetCoefficients().
452: The problem is solved via linearization, by building a pencil (A,B) of
453: size p*n and solving the corresponding GNHEP.
455: Used DS matrices:
456: + DS_MAT_Ex - coefficients of the matrix polynomial
457: . DS_MAT_A - (workspace) first matrix of the linearization
458: . DS_MAT_B - (workspace) second matrix of the linearization
459: . DS_MAT_W - (workspace) right eigenvectors of the linearization
460: - DS_MAT_U - (workspace) left eigenvectors of the linearization
462: Implemented methods:
463: . 0 - QZ iteration on the linearization (_ggev)
465: .seealso: DSCreate(), DSSetType(), DSType, DSPEPSetDegree(), DSPEPSetCoefficients()
466: M*/
467: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
468: {
469: DS_PEP *ctx;
471: PetscNewLog(ds,&ctx);
472: ds->data = (void*)ctx;
474: ds->ops->allocate = DSAllocate_PEP;
475: ds->ops->view = DSView_PEP;
476: ds->ops->vectors = DSVectors_PEP;
477: ds->ops->solve[0] = DSSolve_PEP_QZ;
478: ds->ops->sort = DSSort_PEP;
479: ds->ops->synchronize = DSSynchronize_PEP;
480: ds->ops->destroy = DSDestroy_PEP;
481: ds->ops->matgetsize = DSMatGetSize_PEP;
482: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP);
483: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP);
484: PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP);
485: PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP);
486: PetscFunctionReturn(0);
487: }