Actual source code: dspep.c

slepc-3.17.0 2022-03-31
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: 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: }