Actual source code: dspep.c
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: static PetscErrorCode DSAllocate_PEP(DS ds,PetscInt ld)
20: {
21: DS_PEP *ctx = (DS_PEP*)ds->data;
22: PetscInt i;
24: PetscFunctionBegin;
25: PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
26: PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
27: PetscCall(DSAllocateMat_Private(ds,DS_MAT_Y));
28: for (i=0;i<=ctx->d;i++) PetscCall(DSAllocateMat_Private(ds,DSMatExtra[i]));
29: PetscCall(PetscFree(ds->perm));
30: PetscCall(PetscMalloc1(ld*ctx->d,&ds->perm));
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: static PetscErrorCode DSView_PEP(DS ds,PetscViewer viewer)
35: {
36: DS_PEP *ctx = (DS_PEP*)ds->data;
37: PetscViewerFormat format;
38: PetscInt i;
40: PetscFunctionBegin;
41: PetscCall(PetscViewerGetFormat(viewer,&format));
42: if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
43: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
44: PetscCall(PetscViewerASCIIPrintf(viewer,"polynomial degree: %" PetscInt_FMT "\n",ctx->d));
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
47: for (i=0;i<=ctx->d;i++) PetscCall(DSViewMat(ds,viewer,DSMatExtra[i]));
48: if (ds->state>DS_STATE_INTERMEDIATE) PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: static PetscErrorCode DSVectors_PEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
53: {
54: PetscFunctionBegin;
55: PetscCheck(!rnorm,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
56: switch (mat) {
57: case DS_MAT_X:
58: break;
59: case DS_MAT_Y:
60: break;
61: default:
62: SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
63: }
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static PetscErrorCode DSSort_PEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *kout)
68: {
69: DS_PEP *ctx = (DS_PEP*)ds->data;
70: PetscInt n,i,*perm,told;
71: PetscScalar *A;
73: PetscFunctionBegin;
74: if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
75: n = ds->n*ctx->d;
76: perm = ds->perm;
77: for (i=0;i<n;i++) perm[i] = i;
78: told = ds->t;
79: ds->t = n; /* force the sorting routines to consider d*n eigenvalues */
80: if (rr) PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
81: else PetscCall(DSSortEigenvalues_Private(ds,wr,wi,perm,PETSC_FALSE));
82: ds->t = told; /* restore value of t */
83: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
84: for (i=0;i<n;i++) A[i] = wr[perm[i]];
85: for (i=0;i<n;i++) wr[i] = A[i];
86: for (i=0;i<n;i++) A[i] = wi[perm[i]];
87: for (i=0;i<n;i++) wi[i] = A[i];
88: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
89: PetscCall(DSPermuteColumnsTwo_Private(ds,0,n,ds->n,DS_MAT_X,DS_MAT_Y,perm));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode DSSolve_PEP_QZ(DS ds,PetscScalar *wr,PetscScalar *wi)
94: {
95: DS_PEP *ctx = (DS_PEP*)ds->data;
96: PetscInt i,j,k,off;
97: PetscScalar *A,*B,*W,*X,*U,*Y,*work,*beta,a;
98: const PetscScalar *Ed,*Ei;
99: PetscReal *ca,*cb,*cg,norm,done=1.0;
100: PetscBLASInt n,ld,ldd,nd,lwork,one=1,zero=0,cols;
101: PetscBool useggev3=(ds->method==1)?PETSC_TRUE:PETSC_FALSE;
103: PetscFunctionBegin;
104: PetscCall(PetscBLASIntCast(ds->n*ctx->d,&nd));
105: PetscCall(PetscBLASIntCast(ds->n,&n));
106: PetscCall(PetscBLASIntCast(ds->ld,&ld));
107: PetscCall(PetscBLASIntCast(ds->ld*ctx->d,&ldd));
108: PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
109: PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
110: PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
111: PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
112: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
113: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
115: /* build matrices A and B of the linearization */
116: PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[ctx->d]],&Ed));
117: PetscCall(PetscArrayzero(A,ldd*ldd));
118: if (!ctx->pbc) { /* monomial basis */
119: for (i=0;i<nd-ds->n;i++) A[i+(i+ds->n)*ldd] = 1.0;
120: for (i=0;i<ctx->d;i++) {
121: PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
122: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
123: for (j=0;j<ds->n;j++) PetscCall(PetscArraycpy(A+off+j*ldd,Ei+j*ds->ld,ds->n));
124: PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
125: }
126: } else {
127: ca = ctx->pbc;
128: cb = ca+ctx->d+1;
129: cg = cb+ctx->d+1;
130: for (i=0;i<ds->n;i++) {
131: A[i+(i+ds->n)*ldd] = ca[0];
132: A[i+i*ldd] = cb[0];
133: }
134: for (;i<nd-ds->n;i++) {
135: j = i/ds->n;
136: A[i+(i+ds->n)*ldd] = ca[j];
137: A[i+i*ldd] = cb[j];
138: A[i+(i-ds->n)*ldd] = cg[j];
139: }
140: for (i=0;i<ctx->d-2;i++) {
141: PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
142: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
143: for (j=0;j<ds->n;j++)
144: for (k=0;k<ds->n;k++)
145: A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1];
146: PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
147: }
148: PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
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] = Ei[j*ds->ld+k]*ca[ctx->d-1]-Ed[j*ds->ld+k]*cg[ctx->d-1];
153: PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
154: i++;
155: PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
156: off = i*ds->n*ldd+(ctx->d-1)*ds->n;
157: for (j=0;j<ds->n;j++)
158: for (k=0;k<ds->n;k++)
159: A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1]-Ed[j*ds->ld+k]*cb[ctx->d-1];
160: PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
161: }
162: PetscCall(PetscArrayzero(B,ldd*ldd));
163: for (i=0;i<nd-ds->n;i++) B[i+i*ldd] = 1.0;
164: off = (ctx->d-1)*ds->n*(ldd+1);
165: for (j=0;j<ds->n;j++) {
166: for (i=0;i<ds->n;i++) B[off+i+j*ldd] = -Ed[i+j*ds->ld];
167: }
168: PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[ctx->d]],&Ed));
170: /* solve generalized eigenproblem */
171: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
172: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
173: lwork = -1;
174: #if PetscDefined(USE_COMPLEX)
175: if (useggev3) PetscCallLAPACKInfo("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,NULL,U,&ldd,W,&ldd,&a,&lwork,NULL,&info));
176: else PetscCallLAPACKInfo("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,NULL,U,&ldd,W,&ldd,&a,&lwork,NULL,&info));
177: PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
178: PetscCall(DSAllocateWork_Private(ds,lwork+nd,8*nd,0));
179: beta = ds->work;
180: work = ds->work + nd;
181: if (useggev3) PetscCallLAPACKInfo("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,ds->rwork,&info));
182: else PetscCallLAPACKInfo("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,ds->rwork,&info));
183: #else
184: if (useggev3) PetscCallLAPACKInfo("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,NULL,U,&ldd,W,&ldd,&a,&lwork,&info));
185: else PetscCallLAPACKInfo("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,NULL,U,&ldd,W,&ldd,&a,&lwork,&info));
186: PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
187: PetscCall(DSAllocateWork_Private(ds,lwork+nd,0,0));
188: beta = ds->work;
189: work = ds->work + nd;
190: if (useggev3) PetscCallLAPACKInfo("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
191: else PetscCallLAPACKInfo("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
192: #endif
193: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
194: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
196: /* copy eigenvalues */
197: for (i=0;i<nd;i++) {
198: if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
199: else wr[i] /= beta[i];
200: #if !PetscDefined(USE_COMPLEX)
201: if (beta[i]==0.0) wi[i] = 0.0;
202: else wi[i] /= beta[i];
203: #else
204: if (wi) wi[i] = 0.0;
205: #endif
206: }
208: /* copy and normalize eigenvectors */
209: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
210: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
211: for (j=0;j<nd;j++) {
212: PetscCall(PetscArraycpy(X+j*ds->ld,W+j*ldd,ds->n));
213: PetscCall(PetscArraycpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n));
214: }
215: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
216: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
217: for (j=0;j<nd;j++) {
218: cols = 1;
219: norm = BLASnrm2_(&n,X+j*ds->ld,&one);
220: #if !PetscDefined(USE_COMPLEX)
221: if (wi[j] != 0.0) {
222: norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,X+(j+1)*ds->ld,&one));
223: cols = 2;
224: }
225: #endif
226: PetscCallLAPACKInfo("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,X+j*ds->ld,&ld,&info));
227: norm = BLASnrm2_(&n,Y+j*ds->ld,&one);
228: #if !PetscDefined(USE_COMPLEX)
229: if (wi[j] != 0.0) norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,Y+(j+1)*ds->ld,&one));
230: #endif
231: PetscCallLAPACKInfo("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Y+j*ds->ld,&ld,&info));
232: #if !PetscDefined(USE_COMPLEX)
233: if (wi[j] != 0.0) j++;
234: #endif
235: }
236: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
237: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: #if !PetscDefined(HAVE_MPIUNI)
242: static PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
243: {
244: DS_PEP *ctx = (DS_PEP*)ds->data;
245: PetscInt ld=ds->ld,k=0;
246: PetscMPIInt ldnd,rank,off=0,size,dn;
247: PetscScalar *X,*Y;
249: PetscFunctionBegin;
250: if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
251: if (eigr) k += ctx->d*ds->n;
252: if (eigi) k += ctx->d*ds->n;
253: PetscCall(DSAllocateWork_Private(ds,k,0,0));
254: PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
255: PetscCall(PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd));
256: PetscCall(PetscMPIIntCast(ctx->d*ds->n,&dn));
257: if (ds->state>=DS_STATE_CONDENSED) {
258: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
259: PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
260: }
261: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
262: if (!rank) {
263: if (ds->state>=DS_STATE_CONDENSED) {
264: PetscCallMPI(MPI_Pack(X,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
265: PetscCallMPI(MPI_Pack(Y,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
266: }
267: if (eigr) PetscCallMPI(MPI_Pack(eigr,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
268: #if !PetscDefined(USE_COMPLEX)
269: if (eigi) PetscCallMPI(MPI_Pack(eigi,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
270: #endif
271: }
272: PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
273: if (rank) {
274: if (ds->state>=DS_STATE_CONDENSED) {
275: PetscCallMPI(MPI_Unpack(ds->work,size,&off,X,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
276: PetscCallMPI(MPI_Unpack(ds->work,size,&off,Y,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
277: }
278: if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
279: #if !PetscDefined(USE_COMPLEX)
280: if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
281: #endif
282: }
283: if (ds->state>=DS_STATE_CONDENSED) {
284: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
285: PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
286: }
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
289: #endif
291: static PetscErrorCode DSPEPSetDegree_PEP(DS ds,PetscInt d)
292: {
293: DS_PEP *ctx = (DS_PEP*)ds->data;
295: PetscFunctionBegin;
296: PetscCheck(d>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The degree must be a non-negative integer");
297: PetscCheck(d<DS_NUM_EXTRA,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Only implemented for polynomials of degree at most %d",DS_NUM_EXTRA-1);
298: ctx->d = d;
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: /*@
303: DSPEPSetDegree - Sets the polynomial degree for a `DSPEP`.
305: Logically Collective
307: Input Parameters:
308: + ds - the direct solver context
309: - d - the degree
311: Level: intermediate
313: .seealso: [](sec:ds), `DSPEP`, `DSPEPGetDegree()`
314: @*/
315: PetscErrorCode DSPEPSetDegree(DS ds,PetscInt d)
316: {
317: PetscFunctionBegin;
320: PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
325: {
326: DS_PEP *ctx = (DS_PEP*)ds->data;
328: PetscFunctionBegin;
329: *d = ctx->d;
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: /*@
334: DSPEPGetDegree - Returns the polynomial degree for a `DSPEP`.
336: Not Collective
338: Input Parameter:
339: . ds - the direct solver context
341: Output Parameter:
342: . d - the degree
344: Level: intermediate
346: .seealso: [](sec:ds), `DSPEP`, `DSPEPSetDegree()`
347: @*/
348: PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
349: {
350: PetscFunctionBegin;
352: PetscAssertPointer(d,2);
353: PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode DSPEPSetCoefficients_PEP(DS ds,PetscReal *pbc)
358: {
359: DS_PEP *ctx = (DS_PEP*)ds->data;
360: PetscInt i;
362: PetscFunctionBegin;
363: PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
364: PetscCall(PetscFree(ctx->pbc));
365: PetscCall(PetscMalloc1(3*(ctx->d+1),&ctx->pbc));
366: for (i=0;i<3*(ctx->d+1);i++) ctx->pbc[i] = pbc[i];
367: ds->state = DS_STATE_RAW;
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }
371: /*@
372: DSPEPSetCoefficients - Sets the polynomial basis coefficients for a `DSPEP`.
374: Logically Collective
376: Input Parameters:
377: + ds - the direct solver context
378: - pbc - the polynomial basis coefficients
380: Notes:
381: This function is required only in the case of a polynomial specified in a
382: non-monomial basis, to provide the coefficients that will be used
383: during the linearization, multiplying the identity blocks on the three main
384: diagonal blocks. Depending on the polynomial basis (Chebyshev, Legendre, ...)
385: the coefficients must be different.
387: There must be a total of $3(d+1)$ coefficients, where $d$ is the degree of the
388: polynomial. The coefficients are arranged in three groups, $\alpha_i$, $\beta_i$, and
389: $\gamma_i$, according to the definition of the three-term recurrence. In the case
390: of the monomial basis, $\alpha_i=1$ and $\beta_i=\gamma_i=0$, in which case it is not
391: necessary to invoke this function.
393: Level: advanced
395: .seealso: [](sec:ds), `DSPEP`, `DSPEPGetCoefficients()`, `DSPEPSetDegree()`
396: @*/
397: PetscErrorCode DSPEPSetCoefficients(DS ds,PetscReal pbc[])
398: {
399: PetscFunctionBegin;
401: PetscAssertPointer(pbc,2);
402: PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: static PetscErrorCode DSPEPGetCoefficients_PEP(DS ds,PetscReal *pbc[])
407: {
408: DS_PEP *ctx = (DS_PEP*)ds->data;
409: PetscInt i;
411: PetscFunctionBegin;
412: PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
413: PetscCall(PetscCalloc1(3*(ctx->d+1),pbc));
414: if (ctx->pbc) for (i=0;i<3*(ctx->d+1);i++) (*pbc)[i] = ctx->pbc[i];
415: else for (i=0;i<ctx->d+1;i++) (*pbc)[i] = 1.0;
416: PetscFunctionReturn(PETSC_SUCCESS);
417: }
419: /*@
420: DSPEPGetCoefficients - Returns the polynomial basis coefficients for a `DSPEP`.
422: Not Collective
424: Input Parameter:
425: . ds - the direct solver context
427: Output Parameter:
428: . pbc - the polynomial basis coefficients
430: Note:
431: The returned array has length $3(d+1)$, where $d$ is the degree of the
432: polynomial, and should be freed by the user.
434: Fortran Notes:
435: In Fortran the user must provide in argument `pbc` a sufficiently large array.
437: Level: advanced
439: .seealso: [](sec:ds), `DSPEP`, `DSPEPSetCoefficients()`
440: @*/
441: PetscErrorCode DSPEPGetCoefficients(DS ds,PetscReal *pbc[]) PeNS
442: {
443: PetscFunctionBegin;
445: PetscAssertPointer(pbc,2);
446: PetscUseMethod(ds,"DSPEPGetCoefficients_C",(DS,PetscReal**),(ds,pbc));
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: static PetscErrorCode DSDestroy_PEP(DS ds)
451: {
452: DS_PEP *ctx = (DS_PEP*)ds->data;
454: PetscFunctionBegin;
455: PetscCall(PetscFree(ctx->pbc));
456: PetscCall(PetscFree(ds->data));
457: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL));
458: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL));
459: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL));
460: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL));
461: PetscFunctionReturn(PETSC_SUCCESS);
462: }
464: static PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
465: {
466: DS_PEP *ctx = (DS_PEP*)ds->data;
468: PetscFunctionBegin;
469: PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
470: *rows = ds->n;
471: if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
472: *cols = ds->n;
473: 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;
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: /*MC
478: DSPEP - Dense Polynomial Eigenvalue Problem.
480: Notes:
481: The problem is expressed as $P(\lambda)x = 0$, where $P(\cdot)$ is a
482: polynomial matrix of degree $d$. The eigenvalues $\lambda$ are the arguments
483: returned by `DSSolve()`.
485: The degree of the polynomial, $d$, can be set with `DSPEPSetDegree()`, with
486: the first $d+1$ extra matrices of the `DS` storing the polynomial coefficient.
487: By default, the polynomial is expressed in the monomial basis, but a
488: different basis can be used by setting the corresponding coefficients
489: via `DSPEPSetCoefficients()`.
491: The problem is solved via linearization, by building a pencil $(A,B)$ of
492: size $d\cdot n$ and solving the corresponding `GNHEP`.
494: Used DS matrices:
495: + `DS_MAT_E0` to `DS_MAT_E9` - coefficients of the matrix polynomial
496: . `DS_MAT_X` - right eigenvectors
497: . `DS_MAT_Y` - left eigenvectors
498: . `DS_MAT_A` - (workspace) first matrix of the linearization
499: . `DS_MAT_B` - (workspace) second matrix of the linearization
500: . `DS_MAT_W` - (workspace) right eigenvectors of the linearization
501: - `DS_MAT_U` - (workspace) left eigenvectors of the linearization
503: Implemented methods:
504: . 0 - QZ iteration on the linearization (`_ggev`)
506: Level: beginner
508: .seealso: [](sec:ds), `DSCreate()`, `DSSetType()`, `DSType`, `DSPEPSetDegree()`, `DSPEPSetCoefficients()`
509: M*/
510: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
511: {
512: DS_PEP *ctx;
514: PetscFunctionBegin;
515: PetscCall(PetscNew(&ctx));
516: ds->data = (void*)ctx;
518: ds->ops->allocate = DSAllocate_PEP;
519: ds->ops->view = DSView_PEP;
520: ds->ops->vectors = DSVectors_PEP;
521: ds->ops->solve[0] = DSSolve_PEP_QZ;
522: #if !defined(SLEPC_MISSING_LAPACK_GGES3)
523: ds->ops->solve[1] = DSSolve_PEP_QZ;
524: #endif
525: ds->ops->sort = DSSort_PEP;
526: #if !PetscDefined(HAVE_MPIUNI)
527: ds->ops->synchronize = DSSynchronize_PEP;
528: #endif
529: ds->ops->destroy = DSDestroy_PEP;
530: ds->ops->matgetsize = DSMatGetSize_PEP;
531: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP));
532: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP));
533: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP));
534: PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP));
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }