Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : 
      11             : #include <slepc/private/dsimpl.h>       /*I "slepcds.h" I*/
      12             : #include <slepcblaslapack.h>
      13             : 
      14             : typedef struct {
      15             :   PetscInt  d;              /* polynomial degree */
      16             :   PetscReal *pbc;           /* polynomial basis coefficients */
      17             : } DS_PEP;
      18             : 
      19          24 : static PetscErrorCode DSAllocate_PEP(DS ds,PetscInt ld)
      20             : {
      21          24 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
      22          24 :   PetscInt       i;
      23             : 
      24          24 :   PetscFunctionBegin;
      25          24 :   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
      26          24 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
      27          24 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_Y));
      28         109 :   for (i=0;i<=ctx->d;i++) PetscCall(DSAllocateMat_Private(ds,DSMatExtra[i]));
      29          24 :   PetscCall(PetscFree(ds->perm));
      30          24 :   PetscCall(PetscMalloc1(ld*ctx->d,&ds->perm));
      31          24 :   PetscFunctionReturn(PETSC_SUCCESS);
      32             : }
      33             : 
      34           2 : static PetscErrorCode DSView_PEP(DS ds,PetscViewer viewer)
      35             : {
      36           2 :   DS_PEP            *ctx = (DS_PEP*)ds->data;
      37           2 :   PetscViewerFormat format;
      38           2 :   PetscInt          i;
      39             : 
      40           2 :   PetscFunctionBegin;
      41           2 :   PetscCall(PetscViewerGetFormat(viewer,&format));
      42           2 :   if (format == PETSC_VIEWER_ASCII_INFO) PetscFunctionReturn(PETSC_SUCCESS);
      43           2 :   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
      44           2 :     PetscCall(PetscViewerASCIIPrintf(viewer,"polynomial degree: %" PetscInt_FMT "\n",ctx->d));
      45           2 :     PetscFunctionReturn(PETSC_SUCCESS);
      46             :   }
      47           0 :   for (i=0;i<=ctx->d;i++) PetscCall(DSViewMat(ds,viewer,DSMatExtra[i]));
      48           0 :   if (ds->state>DS_STATE_INTERMEDIATE) PetscCall(DSViewMat(ds,viewer,DS_MAT_X));
      49           0 :   PetscFunctionReturn(PETSC_SUCCESS);
      50             : }
      51             : 
      52          23 : static PetscErrorCode DSVectors_PEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
      53             : {
      54          23 :   PetscFunctionBegin;
      55          23 :   PetscCheck(!rnorm,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
      56          23 :   switch (mat) {
      57             :     case DS_MAT_X:
      58             :       break;
      59             :     case DS_MAT_Y:
      60             :       break;
      61           0 :     default:
      62           0 :       SETERRQ(PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"Invalid mat parameter");
      63             :   }
      64          23 :   PetscFunctionReturn(PETSC_SUCCESS);
      65             : }
      66             : 
      67         322 : static PetscErrorCode DSSort_PEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *kout)
      68             : {
      69         322 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
      70         322 :   PetscInt       n,i,*perm,told;
      71         322 :   PetscScalar    *A;
      72             : 
      73         322 :   PetscFunctionBegin;
      74         322 :   if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
      75         322 :   n = ds->n*ctx->d;
      76         322 :   perm = ds->perm;
      77        8270 :   for (i=0;i<n;i++) perm[i] = i;
      78         322 :   told = ds->t;
      79         322 :   ds->t = n;  /* force the sorting routines to consider d*n eigenvalues */
      80         322 :   if (rr) PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
      81         312 :   else PetscCall(DSSortEigenvalues_Private(ds,wr,wi,perm,PETSC_FALSE));
      82         322 :   ds->t = told;  /* restore value of t */
      83         322 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
      84        8270 :   for (i=0;i<n;i++) A[i]  = wr[perm[i]];
      85        8270 :   for (i=0;i<n;i++) wr[i] = A[i];
      86        8270 :   for (i=0;i<n;i++) A[i]  = wi[perm[i]];
      87        8270 :   for (i=0;i<n;i++) wi[i] = A[i];
      88         322 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
      89         322 :   PetscCall(DSPermuteColumnsTwo_Private(ds,0,n,ds->n,DS_MAT_X,DS_MAT_Y,perm));
      90         322 :   PetscFunctionReturn(PETSC_SUCCESS);
      91             : }
      92             : 
      93         322 : static PetscErrorCode DSSolve_PEP_QZ(DS ds,PetscScalar *wr,PetscScalar *wi)
      94             : {
      95         322 :   DS_PEP            *ctx = (DS_PEP*)ds->data;
      96         322 :   PetscInt          i,j,k,off;
      97         322 :   PetscScalar       *A,*B,*W,*X,*U,*Y,*work,*beta,a;
      98         322 :   const PetscScalar *Ed,*Ei;
      99         322 :   PetscReal         *ca,*cb,*cg,norm,done=1.0;
     100         322 :   PetscBLASInt      info,n,ld,ldd,nd,lwork,one=1,zero=0,cols;
     101         322 :   PetscBool         useggev3=(ds->method==1)?PETSC_TRUE:PETSC_FALSE;
     102             : 
     103         322 :   PetscFunctionBegin;
     104         322 :   PetscCall(PetscBLASIntCast(ds->n*ctx->d,&nd));
     105         322 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     106         322 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     107         322 :   PetscCall(PetscBLASIntCast(ds->ld*ctx->d,&ldd));
     108         322 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
     109         322 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
     110         322 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
     111         322 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
     112         322 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     113         322 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
     114             : 
     115             :   /* build matrices A and B of the linearization */
     116         322 :   PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[ctx->d]],&Ed));
     117         322 :   PetscCall(PetscArrayzero(A,ldd*ldd));
     118         322 :   if (!ctx->pbc) { /* monomial basis */
     119        3274 :     for (i=0;i<nd-ds->n;i++) A[i+(i+ds->n)*ldd] = 1.0;
     120         848 :     for (i=0;i<ctx->d;i++) {
     121         566 :       PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     122         566 :       off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     123        6532 :       for (j=0;j<ds->n;j++) PetscCall(PetscArraycpy(A+off+j*ldd,Ei+j*ds->ld,ds->n));
     124         566 :       PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     125             :     }
     126             :   } else {
     127          40 :     ca = ctx->pbc;
     128          40 :     cb = ca+ctx->d+1;
     129          40 :     cg = cb+ctx->d+1;
     130         746 :     for (i=0;i<ds->n;i++) {
     131         706 :       A[i+(i+ds->n)*ldd] = ca[0];
     132         706 :       A[i+i*ldd] = cb[0];
     133             :     }
     134         610 :     for (;i<nd-ds->n;i++) {
     135         570 :       j = i/ds->n;
     136         570 :       A[i+(i+ds->n)*ldd] = ca[j];
     137         570 :       A[i+i*ldd] = cb[j];
     138         570 :       A[i+(i-ds->n)*ldd] = cg[j];
     139             :     }
     140          80 :     for (i=0;i<ctx->d-2;i++) {
     141          40 :       PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     142          40 :       off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     143         610 :       for (j=0;j<ds->n;j++)
     144       16012 :         for (k=0;k<ds->n;k++)
     145       15442 :           A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1];
     146          40 :       PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     147             :     }
     148          40 :     PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     149          40 :     off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     150         746 :     for (j=0;j<ds->n;j++)
     151       32222 :       for (k=0;k<ds->n;k++)
     152       31516 :         A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1]-Ed[j*ds->ld+k]*cg[ctx->d-1];
     153          40 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     154          40 :     i++;
     155          40 :     PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     156          40 :     off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     157         746 :     for (j=0;j<ds->n;j++)
     158       32222 :       for (k=0;k<ds->n;k++)
     159       31516 :         A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1]-Ed[j*ds->ld+k]*cb[ctx->d-1];
     160          40 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     161             :   }
     162         322 :   PetscCall(PetscArrayzero(B,ldd*ldd));
     163        4590 :   for (i=0;i<nd-ds->n;i++) B[i+i*ldd] = 1.0;
     164         322 :   off = (ctx->d-1)*ds->n*(ldd+1);
     165        4002 :   for (j=0;j<ds->n;j++) {
     166       75000 :     for (i=0;i<ds->n;i++) B[off+i+j*ldd] = -Ed[i+j*ds->ld];
     167             :   }
     168         322 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[ctx->d]],&Ed));
     169             : 
     170             :   /* solve generalized eigenproblem */
     171         322 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
     172         322 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
     173         322 :   lwork = -1;
     174             : #if defined(PETSC_USE_COMPLEX)
     175         322 :   if (useggev3) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,NULL,U,&ldd,W,&ldd,&a,&lwork,NULL,&info));
     176         322 :   else PetscCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,NULL,U,&ldd,W,&ldd,&a,&lwork,NULL,&info));
     177         322 :   PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
     178         322 :   PetscCall(DSAllocateWork_Private(ds,lwork+nd,8*nd,0));
     179         322 :   beta  = ds->work;
     180         322 :   work  = ds->work + nd;
     181         322 :   if (useggev3) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,ds->rwork,&info));
     182         322 :   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         322 :   SlepcCheckLapackInfo(useggev3?"ggev3":"ggev",info);
     194         322 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     195         322 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
     196             : 
     197             :   /* copy eigenvalues */
     198        8270 :   for (i=0;i<nd;i++) {
     199        7948 :     if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
     200        7947 :     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        7948 :     if (wi) wi[i] = 0.0;
     206             : #endif
     207             :   }
     208             : 
     209             :   /* copy and normalize eigenvectors */
     210         322 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
     211         322 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
     212        8270 :   for (j=0;j<nd;j++) {
     213        7948 :     PetscCall(PetscArraycpy(X+j*ds->ld,W+j*ldd,ds->n));
     214        7948 :     PetscCall(PetscArraycpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n));
     215             :   }
     216         322 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
     217         322 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
     218        8270 :   for (j=0;j<nd;j++) {
     219        7948 :     cols = 1;
     220        7948 :     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        7948 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,X+j*ds->ld,&ld,&info));
     228        7948 :     SlepcCheckLapackInfo("lascl",info);
     229        7948 :     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        7948 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Y+j*ds->ld,&ld,&info));
     234        7948 :     SlepcCheckLapackInfo("lascl",info);
     235             : #if !defined(PETSC_USE_COMPLEX)
     236             :     if (wi[j] != 0.0) j++;
     237             : #endif
     238             :   }
     239         322 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
     240         322 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
     241         322 :   PetscFunctionReturn(PETSC_SUCCESS);
     242             : }
     243             : 
     244             : #if !defined(PETSC_HAVE_MPIUNI)
     245          52 : static PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
     246             : {
     247          52 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     248          52 :   PetscInt       ld=ds->ld,k=0;
     249          52 :   PetscMPIInt    ldnd,rank,off=0,size,dn;
     250          52 :   PetscScalar    *X,*Y;
     251             : 
     252          52 :   PetscFunctionBegin;
     253          52 :   if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
     254          52 :   if (eigr) k += ctx->d*ds->n;
     255          52 :   if (eigi) k += ctx->d*ds->n;
     256          52 :   PetscCall(DSAllocateWork_Private(ds,k,0,0));
     257          52 :   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
     258          52 :   PetscCall(PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd));
     259          52 :   PetscCall(PetscMPIIntCast(ctx->d*ds->n,&dn));
     260          52 :   if (ds->state>=DS_STATE_CONDENSED) {
     261          52 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
     262          52 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
     263             :   }
     264          52 :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
     265          52 :   if (!rank) {
     266          26 :     if (ds->state>=DS_STATE_CONDENSED) {
     267          26 :       PetscCallMPI(MPI_Pack(X,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     268          26 :       PetscCallMPI(MPI_Pack(Y,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     269             :     }
     270          26 :     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         104 :   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
     276          52 :   if (rank) {
     277          26 :     if (ds->state>=DS_STATE_CONDENSED) {
     278          26 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,X,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     279          26 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,Y,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     280             :     }
     281          26 :     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          52 :   if (ds->state>=DS_STATE_CONDENSED) {
     287          52 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
     288          52 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
     289             :   }
     290          52 :   PetscFunctionReturn(PETSC_SUCCESS);
     291             : }
     292             : #endif
     293             : 
     294          25 : static PetscErrorCode DSPEPSetDegree_PEP(DS ds,PetscInt d)
     295             : {
     296          25 :   DS_PEP *ctx = (DS_PEP*)ds->data;
     297             : 
     298          25 :   PetscFunctionBegin;
     299          25 :   PetscCheck(d>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The degree must be a non-negative integer");
     300          25 :   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          25 :   ctx->d = d;
     302          25 :   PetscFunctionReturn(PETSC_SUCCESS);
     303             : }
     304             : 
     305             : /*@
     306             :    DSPEPSetDegree - Sets the polynomial degree for a DSPEP.
     307             : 
     308             :    Logically Collective
     309             : 
     310             :    Input Parameters:
     311             : +  ds - the direct solver context
     312             : -  d  - the degree
     313             : 
     314             :    Level: intermediate
     315             : 
     316             : .seealso: DSPEPGetDegree()
     317             : @*/
     318          25 : PetscErrorCode DSPEPSetDegree(DS ds,PetscInt d)
     319             : {
     320          25 :   PetscFunctionBegin;
     321          25 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     322         100 :   PetscValidLogicalCollectiveInt(ds,d,2);
     323          25 :   PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
     324          25 :   PetscFunctionReturn(PETSC_SUCCESS);
     325             : }
     326             : 
     327        1337 : static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
     328             : {
     329        1337 :   DS_PEP *ctx = (DS_PEP*)ds->data;
     330             : 
     331        1337 :   PetscFunctionBegin;
     332        1337 :   *d = ctx->d;
     333        1337 :   PetscFunctionReturn(PETSC_SUCCESS);
     334             : }
     335             : 
     336             : /*@
     337             :    DSPEPGetDegree - Returns the polynomial degree for a DSPEP.
     338             : 
     339             :    Not Collective
     340             : 
     341             :    Input Parameter:
     342             : .  ds - the direct solver context
     343             : 
     344             :    Output Parameters:
     345             : .  d - the degree
     346             : 
     347             :    Level: intermediate
     348             : 
     349             : .seealso: DSPEPSetDegree()
     350             : @*/
     351        1337 : PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
     352             : {
     353        1337 :   PetscFunctionBegin;
     354        1337 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     355        1337 :   PetscAssertPointer(d,2);
     356        1337 :   PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
     357        1337 :   PetscFunctionReturn(PETSC_SUCCESS);
     358             : }
     359             : 
     360          12 : static PetscErrorCode DSPEPSetCoefficients_PEP(DS ds,PetscReal *pbc)
     361             : {
     362          12 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     363          12 :   PetscInt       i;
     364             : 
     365          12 :   PetscFunctionBegin;
     366          12 :   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
     367          12 :   if (ctx->pbc) PetscCall(PetscFree(ctx->pbc));
     368          12 :   PetscCall(PetscMalloc1(3*(ctx->d+1),&ctx->pbc));
     369         156 :   for (i=0;i<3*(ctx->d+1);i++) ctx->pbc[i] = pbc[i];
     370          12 :   ds->state = DS_STATE_RAW;
     371          12 :   PetscFunctionReturn(PETSC_SUCCESS);
     372             : }
     373             : 
     374             : /*@C
     375             :    DSPEPSetCoefficients - Sets the polynomial basis coefficients for a DSPEP.
     376             : 
     377             :    Logically Collective
     378             : 
     379             :    Input Parameters:
     380             : +  ds  - the direct solver context
     381             : -  pbc - the polynomial basis coefficients
     382             : 
     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.
     389             : 
     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.
     395             : 
     396             :    Level: advanced
     397             : 
     398             : .seealso: DSPEPGetCoefficients(), DSPEPSetDegree()
     399             : @*/
     400          12 : PetscErrorCode DSPEPSetCoefficients(DS ds,PetscReal *pbc)
     401             : {
     402          12 :   PetscFunctionBegin;
     403          12 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     404          12 :   PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
     405          12 :   PetscFunctionReturn(PETSC_SUCCESS);
     406             : }
     407             : 
     408           1 : static PetscErrorCode DSPEPGetCoefficients_PEP(DS ds,PetscReal **pbc)
     409             : {
     410           1 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     411           1 :   PetscInt       i;
     412             : 
     413           1 :   PetscFunctionBegin;
     414           1 :   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
     415           1 :   PetscCall(PetscCalloc1(3*(ctx->d+1),pbc));
     416           1 :   if (ctx->pbc) for (i=0;i<3*(ctx->d+1);i++) (*pbc)[i] = ctx->pbc[i];
     417           4 :   else for (i=0;i<ctx->d+1;i++) (*pbc)[i] = 1.0;
     418           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     419             : }
     420             : 
     421             : /*@C
     422             :    DSPEPGetCoefficients - Returns the polynomial basis coefficients for a DSPEP.
     423             : 
     424             :    Not Collective
     425             : 
     426             :    Input Parameter:
     427             : .  ds - the direct solver context
     428             : 
     429             :    Output Parameters:
     430             : .  pbc - the polynomial basis coefficients
     431             : 
     432             :    Note:
     433             :    The returned array has length 3*(d+1) and should be freed by the user.
     434             : 
     435             :    Fortran Notes:
     436             :    The calling sequence from Fortran is
     437             : .vb
     438             :    DSPEPGetCoefficients(eps,pbc,ierr)
     439             :    double precision pbc(d+1) output
     440             : .ve
     441             : 
     442             :    Level: advanced
     443             : 
     444             : .seealso: DSPEPSetCoefficients()
     445             : @*/
     446           1 : PetscErrorCode DSPEPGetCoefficients(DS ds,PetscReal **pbc)
     447             : {
     448           1 :   PetscFunctionBegin;
     449           1 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     450           1 :   PetscAssertPointer(pbc,2);
     451           1 :   PetscUseMethod(ds,"DSPEPGetCoefficients_C",(DS,PetscReal**),(ds,pbc));
     452           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     453             : }
     454             : 
     455          24 : static PetscErrorCode DSDestroy_PEP(DS ds)
     456             : {
     457          24 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     458             : 
     459          24 :   PetscFunctionBegin;
     460          24 :   if (ctx->pbc) PetscCall(PetscFree(ctx->pbc));
     461          24 :   PetscCall(PetscFree(ds->data));
     462          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL));
     463          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL));
     464          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL));
     465          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL));
     466          24 :   PetscFunctionReturn(PETSC_SUCCESS);
     467             : }
     468             : 
     469        1050 : static PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
     470             : {
     471        1050 :   DS_PEP *ctx = (DS_PEP*)ds->data;
     472             : 
     473        1050 :   PetscFunctionBegin;
     474        1050 :   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
     475        1050 :   *rows = ds->n;
     476        1050 :   if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
     477        1050 :   *cols = ds->n;
     478        1050 :   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;
     479        1050 :   PetscFunctionReturn(PETSC_SUCCESS);
     480             : }
     481             : 
     482             : /*MC
     483             :    DSPEP - Dense Polynomial Eigenvalue Problem.
     484             : 
     485             :    Level: beginner
     486             : 
     487             :    Notes:
     488             :    The problem is expressed as P(lambda)*x = 0, where P(.) is a matrix
     489             :    polynomial of degree d. The eigenvalues lambda are the arguments
     490             :    returned by DSSolve().
     491             : 
     492             :    The degree of the polynomial, d, can be set with DSPEPSetDegree(), with
     493             :    the first d+1 extra matrices of the DS defining the matrix polynomial. By
     494             :    default, the polynomial is expressed in the monomial basis, but a
     495             :    different basis can be used by setting the corresponding coefficients
     496             :    via DSPEPSetCoefficients().
     497             : 
     498             :    The problem is solved via linearization, by building a pencil (A,B) of
     499             :    size p*n and solving the corresponding GNHEP.
     500             : 
     501             :    Used DS matrices:
     502             : +  DS_MAT_Ex - coefficients of the matrix polynomial
     503             : .  DS_MAT_A  - (workspace) first matrix of the linearization
     504             : .  DS_MAT_B  - (workspace) second matrix of the linearization
     505             : .  DS_MAT_W  - (workspace) right eigenvectors of the linearization
     506             : -  DS_MAT_U  - (workspace) left eigenvectors of the linearization
     507             : 
     508             :    Implemented methods:
     509             : .  0 - QZ iteration on the linearization (_ggev)
     510             : 
     511             : .seealso: DSCreate(), DSSetType(), DSType, DSPEPSetDegree(), DSPEPSetCoefficients()
     512             : M*/
     513          24 : SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
     514             : {
     515          24 :   DS_PEP         *ctx;
     516             : 
     517          24 :   PetscFunctionBegin;
     518          24 :   PetscCall(PetscNew(&ctx));
     519          24 :   ds->data = (void*)ctx;
     520             : 
     521          24 :   ds->ops->allocate      = DSAllocate_PEP;
     522          24 :   ds->ops->view          = DSView_PEP;
     523          24 :   ds->ops->vectors       = DSVectors_PEP;
     524          24 :   ds->ops->solve[0]      = DSSolve_PEP_QZ;
     525             : #if !defined(SLEPC_MISSING_LAPACK_GGES3)
     526          24 :   ds->ops->solve[1]      = DSSolve_PEP_QZ;
     527             : #endif
     528          24 :   ds->ops->sort          = DSSort_PEP;
     529             : #if !defined(PETSC_HAVE_MPIUNI)
     530          24 :   ds->ops->synchronize   = DSSynchronize_PEP;
     531             : #endif
     532          24 :   ds->ops->destroy       = DSDestroy_PEP;
     533          24 :   ds->ops->matgetsize    = DSMatGetSize_PEP;
     534          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP));
     535          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP));
     536          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP));
     537          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP));
     538          24 :   PetscFunctionReturn(PETSC_SUCCESS);
     539             : }