Actual source code: dspep.c

slepc-3.22.1 2024-10-28
Report Typos and Errors
  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      info,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 defined(PETSC_USE_COMPLEX)
175:   if (useggev3) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,NULL,U,&ldd,W,&ldd,&a,&lwork,NULL,&info));
176:   else PetscCallBLAS("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) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,ds->rwork,&info));
182:   else PetscCallBLAS("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) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,NULL,U,&ldd,W,&ldd,&a,&lwork,&info));
185:   else PetscCallBLAS("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) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
191:   else PetscCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
192: #endif
193:   SlepcCheckLapackInfo(useggev3?"ggev3":"ggev",info);
194:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
195:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));

197:   /* copy eigenvalues */
198:   for (i=0;i<nd;i++) {
199:     if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
200:     else wr[i] /= beta[i];
201: #if !defined(PETSC_USE_COMPLEX)
202:     if (beta[i]==0.0) wi[i] = 0.0;
203:     else wi[i] /= beta[i];
204: #else
205:     if (wi) wi[i] = 0.0;
206: #endif
207:   }

209:   /* copy and normalize eigenvectors */
210:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
211:   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
212:   for (j=0;j<nd;j++) {
213:     PetscCall(PetscArraycpy(X+j*ds->ld,W+j*ldd,ds->n));
214:     PetscCall(PetscArraycpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n));
215:   }
216:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
217:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
218:   for (j=0;j<nd;j++) {
219:     cols = 1;
220:     norm = BLASnrm2_(&n,X+j*ds->ld,&one);
221: #if !defined(PETSC_USE_COMPLEX)
222:     if (wi[j] != 0.0) {
223:       norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,X+(j+1)*ds->ld,&one));
224:       cols = 2;
225:     }
226: #endif
227:     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,X+j*ds->ld,&ld,&info));
228:     SlepcCheckLapackInfo("lascl",info);
229:     norm = BLASnrm2_(&n,Y+j*ds->ld,&one);
230: #if !defined(PETSC_USE_COMPLEX)
231:     if (wi[j] != 0.0) norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,Y+(j+1)*ds->ld,&one));
232: #endif
233:     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Y+j*ds->ld,&ld,&info));
234:     SlepcCheckLapackInfo("lascl",info);
235: #if !defined(PETSC_USE_COMPLEX)
236:     if (wi[j] != 0.0) j++;
237: #endif
238:   }
239:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
240:   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
241:   PetscFunctionReturn(PETSC_SUCCESS);
242: }

244: #if !defined(PETSC_HAVE_MPIUNI)
245: static PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
246: {
247:   DS_PEP         *ctx = (DS_PEP*)ds->data;
248:   PetscInt       ld=ds->ld,k=0;
249:   PetscMPIInt    ldnd,rank,off=0,size,dn;
250:   PetscScalar    *X,*Y;

252:   PetscFunctionBegin;
253:   if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
254:   if (eigr) k += ctx->d*ds->n;
255:   if (eigi) k += ctx->d*ds->n;
256:   PetscCall(DSAllocateWork_Private(ds,k,0,0));
257:   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
258:   PetscCall(PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd));
259:   PetscCall(PetscMPIIntCast(ctx->d*ds->n,&dn));
260:   if (ds->state>=DS_STATE_CONDENSED) {
261:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
262:     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
263:   }
264:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
265:   if (!rank) {
266:     if (ds->state>=DS_STATE_CONDENSED) {
267:       PetscCallMPI(MPI_Pack(X,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
268:       PetscCallMPI(MPI_Pack(Y,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
269:     }
270:     if (eigr) PetscCallMPI(MPI_Pack(eigr,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
271: #if !defined(PETSC_USE_COMPLEX)
272:     if (eigi) PetscCallMPI(MPI_Pack(eigi,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
273: #endif
274:   }
275:   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
276:   if (rank) {
277:     if (ds->state>=DS_STATE_CONDENSED) {
278:       PetscCallMPI(MPI_Unpack(ds->work,size,&off,X,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
279:       PetscCallMPI(MPI_Unpack(ds->work,size,&off,Y,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
280:     }
281:     if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
282: #if !defined(PETSC_USE_COMPLEX)
283:     if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
284: #endif
285:   }
286:   if (ds->state>=DS_STATE_CONDENSED) {
287:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
288:     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
289:   }
290:   PetscFunctionReturn(PETSC_SUCCESS);
291: }
292: #endif

294: static PetscErrorCode DSPEPSetDegree_PEP(DS ds,PetscInt d)
295: {
296:   DS_PEP *ctx = (DS_PEP*)ds->data;

298:   PetscFunctionBegin;
299:   PetscCheck(d>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The degree must be a non-negative integer");
300:   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);
301:   ctx->d = d;
302:   PetscFunctionReturn(PETSC_SUCCESS);
303: }

305: /*@
306:    DSPEPSetDegree - Sets the polynomial degree for a DSPEP.

308:    Logically Collective

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: {
320:   PetscFunctionBegin;
323:   PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
324:   PetscFunctionReturn(PETSC_SUCCESS);
325: }

327: static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
328: {
329:   DS_PEP *ctx = (DS_PEP*)ds->data;

331:   PetscFunctionBegin;
332:   *d = ctx->d;
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }

336: /*@
337:    DSPEPGetDegree - Returns the polynomial degree for a DSPEP.

339:    Not Collective

341:    Input Parameter:
342: .  ds - the direct solver context

344:    Output Parameters:
345: .  d - the degree

347:    Level: intermediate

349: .seealso: DSPEPSetDegree()
350: @*/
351: PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
352: {
353:   PetscFunctionBegin;
355:   PetscAssertPointer(d,2);
356:   PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }

360: static PetscErrorCode DSPEPSetCoefficients_PEP(DS ds,PetscReal *pbc)
361: {
362:   DS_PEP         *ctx = (DS_PEP*)ds->data;
363:   PetscInt       i;

365:   PetscFunctionBegin;
366:   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
367:   if (ctx->pbc) PetscCall(PetscFree(ctx->pbc));
368:   PetscCall(PetscMalloc1(3*(ctx->d+1),&ctx->pbc));
369:   for (i=0;i<3*(ctx->d+1);i++) ctx->pbc[i] = pbc[i];
370:   ds->state = DS_STATE_RAW;
371:   PetscFunctionReturn(PETSC_SUCCESS);
372: }

374: /*@
375:    DSPEPSetCoefficients - Sets the polynomial basis coefficients for a DSPEP.

377:    Logically Collective

379:    Input Parameters:
380: +  ds  - the direct solver context
381: -  pbc - the polynomial basis coefficients

383:    Notes:
384:    This function is required only in the case of a polynomial specified in a
385:    non-monomial basis, to provide the coefficients that will be used
386:    during the linearization, multiplying the identity blocks on the three main
387:    diagonal blocks. Depending on the polynomial basis (Chebyshev, Legendre, ...)
388:    the coefficients must be different.

390:    There must be a total of 3*(d+1) coefficients, where d is the degree of the
391:    polynomial. The coefficients are arranged in three groups, alpha, beta, and
392:    gamma, according to the definition of the three-term recurrence. In the case
393:    of the monomial basis, alpha=1 and beta=gamma=0, in which case it is not
394:    necessary to invoke this function.

396:    Level: advanced

398: .seealso: DSPEPGetCoefficients(), DSPEPSetDegree()
399: @*/
400: PetscErrorCode DSPEPSetCoefficients(DS ds,PetscReal pbc[])
401: {
402:   PetscFunctionBegin;
404:   PetscAssertPointer(pbc,2);
405:   PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
406:   PetscFunctionReturn(PETSC_SUCCESS);
407: }

409: static PetscErrorCode DSPEPGetCoefficients_PEP(DS ds,PetscReal **pbc)
410: {
411:   DS_PEP         *ctx = (DS_PEP*)ds->data;
412:   PetscInt       i;

414:   PetscFunctionBegin;
415:   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
416:   PetscCall(PetscCalloc1(3*(ctx->d+1),pbc));
417:   if (ctx->pbc) for (i=0;i<3*(ctx->d+1);i++) (*pbc)[i] = ctx->pbc[i];
418:   else for (i=0;i<ctx->d+1;i++) (*pbc)[i] = 1.0;
419:   PetscFunctionReturn(PETSC_SUCCESS);
420: }

422: /*@C
423:    DSPEPGetCoefficients - Returns the polynomial basis coefficients for a DSPEP.

425:    Not Collective

427:    Input Parameter:
428: .  ds - the direct solver context

430:    Output Parameters:
431: .  pbc - the polynomial basis coefficients

433:    Note:
434:    The returned array has length 3*(d+1) and should be freed by the user.

436:    Fortran Notes:
437:    The calling sequence from Fortran is
438: .vb
439:    DSPEPGetCoefficients(eps,pbc,ierr)
440:    double precision pbc(d+1) output
441: .ve

443:    Level: advanced

445: .seealso: DSPEPSetCoefficients()
446: @*/
447: PetscErrorCode DSPEPGetCoefficients(DS ds,PetscReal **pbc)
448: {
449:   PetscFunctionBegin;
451:   PetscAssertPointer(pbc,2);
452:   PetscUseMethod(ds,"DSPEPGetCoefficients_C",(DS,PetscReal**),(ds,pbc));
453:   PetscFunctionReturn(PETSC_SUCCESS);
454: }

456: static PetscErrorCode DSDestroy_PEP(DS ds)
457: {
458:   DS_PEP         *ctx = (DS_PEP*)ds->data;

460:   PetscFunctionBegin;
461:   if (ctx->pbc) PetscCall(PetscFree(ctx->pbc));
462:   PetscCall(PetscFree(ds->data));
463:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL));
464:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL));
465:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL));
466:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL));
467:   PetscFunctionReturn(PETSC_SUCCESS);
468: }

470: static PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
471: {
472:   DS_PEP *ctx = (DS_PEP*)ds->data;

474:   PetscFunctionBegin;
475:   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
476:   *rows = ds->n;
477:   if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
478:   *cols = ds->n;
479:   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;
480:   PetscFunctionReturn(PETSC_SUCCESS);
481: }

483: /*MC
484:    DSPEP - Dense Polynomial Eigenvalue Problem.

486:    Level: beginner

488:    Notes:
489:    The problem is expressed as P(lambda)*x = 0, where P(.) is a matrix
490:    polynomial of degree d. The eigenvalues lambda are the arguments
491:    returned by DSSolve().

493:    The degree of the polynomial, d, can be set with DSPEPSetDegree(), with
494:    the first d+1 extra matrices of the DS defining the matrix polynomial. By
495:    default, the polynomial is expressed in the monomial basis, but a
496:    different basis can be used by setting the corresponding coefficients
497:    via DSPEPSetCoefficients().

499:    The problem is solved via linearization, by building a pencil (A,B) of
500:    size p*n and solving the corresponding GNHEP.

502:    Used DS matrices:
503: +  DS_MAT_Ex - coefficients of the matrix polynomial
504: .  DS_MAT_X  - right eigenvectors
505: .  DS_MAT_Y  - left eigenvectors
506: .  DS_MAT_A  - (workspace) first matrix of the linearization
507: .  DS_MAT_B  - (workspace) second matrix of the linearization
508: .  DS_MAT_W  - (workspace) right eigenvectors of the linearization
509: -  DS_MAT_U  - (workspace) left eigenvectors of the linearization

511:    Implemented methods:
512: .  0 - QZ iteration on the linearization (_ggev)

514: .seealso: DSCreate(), DSSetType(), DSType, DSPEPSetDegree(), DSPEPSetCoefficients()
515: M*/
516: SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
517: {
518:   DS_PEP         *ctx;

520:   PetscFunctionBegin;
521:   PetscCall(PetscNew(&ctx));
522:   ds->data = (void*)ctx;

524:   ds->ops->allocate      = DSAllocate_PEP;
525:   ds->ops->view          = DSView_PEP;
526:   ds->ops->vectors       = DSVectors_PEP;
527:   ds->ops->solve[0]      = DSSolve_PEP_QZ;
528: #if !defined(SLEPC_MISSING_LAPACK_GGES3)
529:   ds->ops->solve[1]      = DSSolve_PEP_QZ;
530: #endif
531:   ds->ops->sort          = DSSort_PEP;
532: #if !defined(PETSC_HAVE_MPIUNI)
533:   ds->ops->synchronize   = DSSynchronize_PEP;
534: #endif
535:   ds->ops->destroy       = DSDestroy_PEP;
536:   ds->ops->matgetsize    = DSMatGetSize_PEP;
537:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP));
538:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP));
539:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP));
540:   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP));
541:   PetscFunctionReturn(PETSC_SUCCESS);
542: }