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          13 : static PetscErrorCode DSAllocate_PEP(DS ds,PetscInt ld)
      20             : {
      21          13 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
      22          13 :   PetscInt       i;
      23             : 
      24          13 :   PetscFunctionBegin;
      25          13 :   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
      26          13 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_X));
      27          13 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_Y));
      28          55 :   for (i=0;i<=ctx->d;i++) PetscCall(DSAllocateMat_Private(ds,DSMatExtra[i]));
      29          13 :   PetscCall(PetscFree(ds->perm));
      30          13 :   PetscCall(PetscMalloc1(ld*ctx->d,&ds->perm));
      31          13 :   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           3 : static PetscErrorCode DSVectors_PEP(DS ds,DSMatType mat,PetscInt *j,PetscReal *rnorm)
      53             : {
      54           3 :   PetscFunctionBegin;
      55           3 :   PetscCheck(!rnorm,PetscObjectComm((PetscObject)ds),PETSC_ERR_SUP,"Not implemented yet");
      56           3 :   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           3 :   PetscFunctionReturn(PETSC_SUCCESS);
      65             : }
      66             : 
      67         267 : static PetscErrorCode DSSort_PEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *kout)
      68             : {
      69         267 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
      70         267 :   PetscInt       n,i,*perm,told;
      71         267 :   PetscScalar    *A;
      72             : 
      73         267 :   PetscFunctionBegin;
      74         267 :   if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
      75         267 :   n = ds->n*ctx->d;
      76         267 :   perm = ds->perm;
      77        5813 :   for (i=0;i<n;i++) perm[i] = i;
      78         267 :   told = ds->t;
      79         267 :   ds->t = n;  /* force the sorting routines to consider d*n eigenvalues */
      80         267 :   if (rr) PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
      81         267 :   else PetscCall(DSSortEigenvalues_Private(ds,wr,wi,perm,PETSC_FALSE));
      82         267 :   ds->t = told;  /* restore value of t */
      83         267 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
      84        5813 :   for (i=0;i<n;i++) A[i]  = wr[perm[i]];
      85        5813 :   for (i=0;i<n;i++) wr[i] = A[i];
      86        5813 :   for (i=0;i<n;i++) A[i]  = wi[perm[i]];
      87        5813 :   for (i=0;i<n;i++) wi[i] = A[i];
      88         267 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
      89         267 :   PetscCall(DSPermuteColumnsTwo_Private(ds,0,n,ds->n,DS_MAT_X,DS_MAT_Y,perm));
      90         267 :   PetscFunctionReturn(PETSC_SUCCESS);
      91             : }
      92             : 
      93         267 : static PetscErrorCode DSSolve_PEP_QZ(DS ds,PetscScalar *wr,PetscScalar *wi)
      94             : {
      95         267 :   DS_PEP            *ctx = (DS_PEP*)ds->data;
      96         267 :   PetscInt          i,j,k,off;
      97         267 :   PetscScalar       *A,*B,*W,*X,*U,*Y,*work,*beta,a;
      98         267 :   const PetscScalar *Ed,*Ei;
      99         267 :   PetscReal         *ca,*cb,*cg,norm,done=1.0;
     100         267 :   PetscBLASInt      info,n,ld,ldd,nd,lwork,one=1,zero=0,cols;
     101         267 :   PetscBool         useggev3=(ds->method==1)?PETSC_TRUE:PETSC_FALSE;
     102             : 
     103         267 :   PetscFunctionBegin;
     104         267 :   PetscCall(PetscBLASIntCast(ds->n*ctx->d,&nd));
     105         267 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     106         267 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     107         267 :   PetscCall(PetscBLASIntCast(ds->ld*ctx->d,&ldd));
     108         267 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
     109         267 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
     110         267 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
     111         267 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
     112         267 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     113         267 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
     114             : 
     115             :   /* build matrices A and B of the linearization */
     116         267 :   PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[ctx->d]],&Ed));
     117         267 :   PetscCall(PetscArrayzero(A,ldd*ldd));
     118         267 :   if (!ctx->pbc) { /* monomial basis */
     119        2724 :     for (i=0;i<nd-ds->n;i++) A[i+(i+ds->n)*ldd] = 1.0;
     120         719 :     for (i=0;i<ctx->d;i++) {
     121         480 :       PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     122         480 :       off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     123        5432 :       for (j=0;j<ds->n;j++) PetscCall(PetscArraycpy(A+off+j*ldd,Ei+j*ds->ld,ds->n));
     124         480 :       PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     125             :     }
     126             :   } else {
     127          28 :     ca = ctx->pbc;
     128          28 :     cb = ca+ctx->d+1;
     129          28 :     cg = cb+ctx->d+1;
     130         226 :     for (i=0;i<ds->n;i++) {
     131         198 :       A[i+(i+ds->n)*ldd] = ca[0];
     132         198 :       A[i+i*ldd] = cb[0];
     133             :     }
     134         226 :     for (;i<nd-ds->n;i++) {
     135         198 :       j = i/ds->n;
     136         198 :       A[i+(i+ds->n)*ldd] = ca[j];
     137         198 :       A[i+i*ldd] = cb[j];
     138         198 :       A[i+(i-ds->n)*ldd] = cg[j];
     139             :     }
     140          56 :     for (i=0;i<ctx->d-2;i++) {
     141          28 :       PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     142          28 :       off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     143         226 :       for (j=0;j<ds->n;j++)
     144        1976 :         for (k=0;k<ds->n;k++)
     145        1778 :           A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1];
     146          28 :       PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     147             :     }
     148          28 :     PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     149          28 :     off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     150         226 :     for (j=0;j<ds->n;j++)
     151        1976 :       for (k=0;k<ds->n;k++)
     152        1778 :         A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1]-Ed[j*ds->ld+k]*cg[ctx->d-1];
     153          28 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     154          28 :     i++;
     155          28 :     PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     156          28 :     off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     157         226 :     for (j=0;j<ds->n;j++)
     158        1976 :       for (k=0;k<ds->n;k++)
     159        1778 :         A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1]-Ed[j*ds->ld+k]*cb[ctx->d-1];
     160          28 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     161             :   }
     162         267 :   PetscCall(PetscArrayzero(B,ldd*ldd));
     163        3148 :   for (i=0;i<nd-ds->n;i++) B[i+i*ldd] = 1.0;
     164         267 :   off = (ctx->d-1)*ds->n*(ldd+1);
     165        2932 :   for (j=0;j<ds->n;j++) {
     166       36796 :     for (i=0;i<ds->n;i++) B[off+i+j*ldd] = -Ed[i+j*ds->ld];
     167             :   }
     168         267 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[ctx->d]],&Ed));
     169             : 
     170             :   /* solve generalized eigenproblem */
     171         267 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
     172         267 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
     173         267 :   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         267 :   if (useggev3) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,NULL,U,&ldd,W,&ldd,&a,&lwork,&info));
     185         267 :   else PetscCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,NULL,U,&ldd,W,&ldd,&a,&lwork,&info));
     186         267 :   PetscCall(PetscBLASIntCast((PetscInt)a,&lwork));
     187         267 :   PetscCall(DSAllocateWork_Private(ds,lwork+nd,0,0));
     188         267 :   beta = ds->work;
     189         267 :   work = ds->work + nd;
     190         267 :   if (useggev3) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
     191         267 :   else PetscCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,wi,beta,U,&ldd,W,&ldd,work,&lwork,&info));
     192             : #endif
     193         267 :   SlepcCheckLapackInfo(useggev3?"ggev3":"ggev",info);
     194         267 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     195         267 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
     196             : 
     197             :   /* copy eigenvalues */
     198        5813 :   for (i=0;i<nd;i++) {
     199        5546 :     if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
     200        5546 :     else wr[i] /= beta[i];
     201             : #if !defined(PETSC_USE_COMPLEX)
     202        5546 :     if (beta[i]==0.0) wi[i] = 0.0;
     203        5546 :     else wi[i] /= beta[i];
     204             : #else
     205             :     if (wi) wi[i] = 0.0;
     206             : #endif
     207             :   }
     208             : 
     209             :   /* copy and normalize eigenvectors */
     210         267 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
     211         267 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
     212        5813 :   for (j=0;j<nd;j++) {
     213        5546 :     PetscCall(PetscArraycpy(X+j*ds->ld,W+j*ldd,ds->n));
     214        5546 :     PetscCall(PetscArraycpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n));
     215             :   }
     216         267 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
     217         267 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
     218        5350 :   for (j=0;j<nd;j++) {
     219        5083 :     cols = 1;
     220        5083 :     norm = BLASnrm2_(&n,X+j*ds->ld,&one);
     221             : #if !defined(PETSC_USE_COMPLEX)
     222        5083 :     if (wi[j] != 0.0) {
     223         463 :       norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,X+(j+1)*ds->ld,&one));
     224         463 :       cols = 2;
     225             :     }
     226             : #endif
     227        5083 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,X+j*ds->ld,&ld,&info));
     228        5083 :     SlepcCheckLapackInfo("lascl",info);
     229        5083 :     norm = BLASnrm2_(&n,Y+j*ds->ld,&one);
     230             : #if !defined(PETSC_USE_COMPLEX)
     231        5083 :     if (wi[j] != 0.0) norm = SlepcAbsEigenvalue(norm,BLASnrm2_(&n,Y+(j+1)*ds->ld,&one));
     232             : #endif
     233        5083 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Y+j*ds->ld,&ld,&info));
     234        5083 :     SlepcCheckLapackInfo("lascl",info);
     235             : #if !defined(PETSC_USE_COMPLEX)
     236        5083 :     if (wi[j] != 0.0) j++;
     237             : #endif
     238             :   }
     239         267 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
     240         267 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
     241         267 :   PetscFunctionReturn(PETSC_SUCCESS);
     242             : }
     243             : 
     244             : #if !defined(PETSC_HAVE_MPIUNI)
     245          58 : static PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
     246             : {
     247          58 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     248          58 :   PetscInt       ld=ds->ld,k=0;
     249          58 :   PetscMPIInt    ldnd,rank,off=0,size,dn;
     250          58 :   PetscScalar    *X,*Y;
     251             : 
     252          58 :   PetscFunctionBegin;
     253          58 :   if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
     254          58 :   if (eigr) k += ctx->d*ds->n;
     255          58 :   if (eigi) k += ctx->d*ds->n;
     256          58 :   PetscCall(DSAllocateWork_Private(ds,k,0,0));
     257          58 :   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
     258          58 :   PetscCall(PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd));
     259          58 :   PetscCall(PetscMPIIntCast(ctx->d*ds->n,&dn));
     260          58 :   if (ds->state>=DS_STATE_CONDENSED) {
     261          58 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
     262          58 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
     263             :   }
     264          58 :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
     265          58 :   if (!rank) {
     266          29 :     if (ds->state>=DS_STATE_CONDENSED) {
     267          29 :       PetscCallMPI(MPI_Pack(X,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     268          29 :       PetscCallMPI(MPI_Pack(Y,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     269             :     }
     270          29 :     if (eigr) PetscCallMPI(MPI_Pack(eigr,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     271             : #if !defined(PETSC_USE_COMPLEX)
     272          29 :     if (eigi) PetscCallMPI(MPI_Pack(eigi,dn,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     273             : #endif
     274             :   }
     275         116 :   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
     276          58 :   if (rank) {
     277          29 :     if (ds->state>=DS_STATE_CONDENSED) {
     278          29 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,X,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     279          29 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,Y,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     280             :     }
     281          29 :     if (eigr) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigr,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     282             : #if !defined(PETSC_USE_COMPLEX)
     283          29 :     if (eigi) PetscCallMPI(MPI_Unpack(ds->work,size,&off,eigi,dn,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     284             : #endif
     285             :   }
     286          58 :   if (ds->state>=DS_STATE_CONDENSED) {
     287          58 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
     288          58 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
     289             :   }
     290          58 :   PetscFunctionReturn(PETSC_SUCCESS);
     291             : }
     292             : #endif
     293             : 
     294          14 : static PetscErrorCode DSPEPSetDegree_PEP(DS ds,PetscInt d)
     295             : {
     296          14 :   DS_PEP *ctx = (DS_PEP*)ds->data;
     297             : 
     298          14 :   PetscFunctionBegin;
     299          14 :   PetscCheck(d>=0,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_OUTOFRANGE,"The degree must be a non-negative integer");
     300          14 :   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          14 :   ctx->d = d;
     302          14 :   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          14 : PetscErrorCode DSPEPSetDegree(DS ds,PetscInt d)
     319             : {
     320          14 :   PetscFunctionBegin;
     321          14 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     322          42 :   PetscValidLogicalCollectiveInt(ds,d,2);
     323          14 :   PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
     324          14 :   PetscFunctionReturn(PETSC_SUCCESS);
     325             : }
     326             : 
     327        1095 : static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
     328             : {
     329        1095 :   DS_PEP *ctx = (DS_PEP*)ds->data;
     330             : 
     331        1095 :   PetscFunctionBegin;
     332        1095 :   *d = ctx->d;
     333        1095 :   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        1095 : PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
     352             : {
     353        1095 :   PetscFunctionBegin;
     354        1095 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     355        1095 :   PetscAssertPointer(d,2);
     356        1095 :   PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
     357        1095 :   PetscFunctionReturn(PETSC_SUCCESS);
     358             : }
     359             : 
     360           2 : static PetscErrorCode DSPEPSetCoefficients_PEP(DS ds,PetscReal *pbc)
     361             : {
     362           2 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     363           2 :   PetscInt       i;
     364             : 
     365           2 :   PetscFunctionBegin;
     366           2 :   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
     367           2 :   if (ctx->pbc) PetscCall(PetscFree(ctx->pbc));
     368           2 :   PetscCall(PetscMalloc1(3*(ctx->d+1),&ctx->pbc));
     369          26 :   for (i=0;i<3*(ctx->d+1);i++) ctx->pbc[i] = pbc[i];
     370           2 :   ds->state = DS_STATE_RAW;
     371           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     372             : }
     373             : 
     374             : /*@
     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           2 : PetscErrorCode DSPEPSetCoefficients(DS ds,PetscReal pbc[])
     401             : {
     402           2 :   PetscFunctionBegin;
     403           2 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     404           2 :   PetscAssertPointer(pbc,2);
     405           2 :   PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
     406           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     407             : }
     408             : 
     409           1 : static PetscErrorCode DSPEPGetCoefficients_PEP(DS ds,PetscReal **pbc)
     410             : {
     411           1 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     412           1 :   PetscInt       i;
     413             : 
     414           1 :   PetscFunctionBegin;
     415           1 :   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"Must first specify the polynomial degree via DSPEPSetDegree()");
     416           1 :   PetscCall(PetscCalloc1(3*(ctx->d+1),pbc));
     417           1 :   if (ctx->pbc) for (i=0;i<3*(ctx->d+1);i++) (*pbc)[i] = ctx->pbc[i];
     418           4 :   else for (i=0;i<ctx->d+1;i++) (*pbc)[i] = 1.0;
     419           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     420             : }
     421             : 
     422             : /*@C
     423             :    DSPEPGetCoefficients - Returns the polynomial basis coefficients for a DSPEP.
     424             : 
     425             :    Not Collective
     426             : 
     427             :    Input Parameter:
     428             : .  ds - the direct solver context
     429             : 
     430             :    Output Parameters:
     431             : .  pbc - the polynomial basis coefficients
     432             : 
     433             :    Note:
     434             :    The returned array has length 3*(d+1) and should be freed by the user.
     435             : 
     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
     442             : 
     443             :    Level: advanced
     444             : 
     445             : .seealso: DSPEPSetCoefficients()
     446             : @*/
     447           1 : PetscErrorCode DSPEPGetCoefficients(DS ds,PetscReal **pbc)
     448             : {
     449           1 :   PetscFunctionBegin;
     450           1 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     451           1 :   PetscAssertPointer(pbc,2);
     452           1 :   PetscUseMethod(ds,"DSPEPGetCoefficients_C",(DS,PetscReal**),(ds,pbc));
     453           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     454             : }
     455             : 
     456          13 : static PetscErrorCode DSDestroy_PEP(DS ds)
     457             : {
     458          13 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     459             : 
     460          13 :   PetscFunctionBegin;
     461          13 :   if (ctx->pbc) PetscCall(PetscFree(ctx->pbc));
     462          13 :   PetscCall(PetscFree(ds->data));
     463          13 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL));
     464          13 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL));
     465          13 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL));
     466          13 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL));
     467          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     468             : }
     469             : 
     470         850 : static PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
     471             : {
     472         850 :   DS_PEP *ctx = (DS_PEP*)ds->data;
     473             : 
     474         850 :   PetscFunctionBegin;
     475         850 :   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
     476         850 :   *rows = ds->n;
     477         850 :   if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
     478         850 :   *cols = ds->n;
     479         850 :   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         850 :   PetscFunctionReturn(PETSC_SUCCESS);
     481             : }
     482             : 
     483             : /*MC
     484             :    DSPEP - Dense Polynomial Eigenvalue Problem.
     485             : 
     486             :    Level: beginner
     487             : 
     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().
     492             : 
     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().
     498             : 
     499             :    The problem is solved via linearization, by building a pencil (A,B) of
     500             :    size p*n and solving the corresponding GNHEP.
     501             : 
     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
     510             : 
     511             :    Implemented methods:
     512             : .  0 - QZ iteration on the linearization (_ggev)
     513             : 
     514             : .seealso: DSCreate(), DSSetType(), DSType, DSPEPSetDegree(), DSPEPSetCoefficients()
     515             : M*/
     516          13 : SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
     517             : {
     518          13 :   DS_PEP         *ctx;
     519             : 
     520          13 :   PetscFunctionBegin;
     521          13 :   PetscCall(PetscNew(&ctx));
     522          13 :   ds->data = (void*)ctx;
     523             : 
     524          13 :   ds->ops->allocate      = DSAllocate_PEP;
     525          13 :   ds->ops->view          = DSView_PEP;
     526          13 :   ds->ops->vectors       = DSVectors_PEP;
     527          13 :   ds->ops->solve[0]      = DSSolve_PEP_QZ;
     528             : #if !defined(SLEPC_MISSING_LAPACK_GGES3)
     529          13 :   ds->ops->solve[1]      = DSSolve_PEP_QZ;
     530             : #endif
     531          13 :   ds->ops->sort          = DSSort_PEP;
     532             : #if !defined(PETSC_HAVE_MPIUNI)
     533          13 :   ds->ops->synchronize   = DSSynchronize_PEP;
     534             : #endif
     535          13 :   ds->ops->destroy       = DSDestroy_PEP;
     536          13 :   ds->ops->matgetsize    = DSMatGetSize_PEP;
     537          13 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP));
     538          13 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP));
     539          13 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP));
     540          13 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP));
     541          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     542             : }