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         314 : static PetscErrorCode DSSort_PEP(DS ds,PetscScalar *wr,PetscScalar *wi,PetscScalar *rr,PetscScalar *ri,PetscInt *kout)
      68             : {
      69         314 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
      70         314 :   PetscInt       n,i,*perm,told;
      71         314 :   PetscScalar    *A;
      72             : 
      73         314 :   PetscFunctionBegin;
      74         314 :   if (!ds->sc) PetscFunctionReturn(PETSC_SUCCESS);
      75         314 :   n = ds->n*ctx->d;
      76         314 :   perm = ds->perm;
      77        7963 :   for (i=0;i<n;i++) perm[i] = i;
      78         314 :   told = ds->t;
      79         314 :   ds->t = n;  /* force the sorting routines to consider d*n eigenvalues */
      80         314 :   if (rr) PetscCall(DSSortEigenvalues_Private(ds,rr,ri,perm,PETSC_FALSE));
      81         304 :   else PetscCall(DSSortEigenvalues_Private(ds,wr,wi,perm,PETSC_FALSE));
      82         314 :   ds->t = told;  /* restore value of t */
      83         314 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
      84        7963 :   for (i=0;i<n;i++) A[i]  = wr[perm[i]];
      85        7963 :   for (i=0;i<n;i++) wr[i] = A[i];
      86        7963 :   for (i=0;i<n;i++) A[i]  = wi[perm[i]];
      87        7963 :   for (i=0;i<n;i++) wi[i] = A[i];
      88         314 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
      89         314 :   PetscCall(DSPermuteColumnsTwo_Private(ds,0,n,ds->n,DS_MAT_X,DS_MAT_Y,perm));
      90         314 :   PetscFunctionReturn(PETSC_SUCCESS);
      91             : }
      92             : 
      93         314 : static PetscErrorCode DSSolve_PEP_QZ(DS ds,PetscScalar *wr,PetscScalar *wi)
      94             : {
      95         314 :   DS_PEP            *ctx = (DS_PEP*)ds->data;
      96         314 :   PetscInt          i,j,k,off;
      97         314 :   PetscScalar       *A,*B,*W,*X,*U,*Y,*work,*beta,a;
      98         314 :   const PetscScalar *Ed,*Ei;
      99         314 :   PetscReal         *ca,*cb,*cg,norm,done=1.0;
     100         314 :   PetscBLASInt      info,n,ld,ldd,nd,lwork,one=1,zero=0,cols;
     101         314 :   PetscBool         useggev3=(ds->method==1)?PETSC_TRUE:PETSC_FALSE;
     102             : 
     103         314 :   PetscFunctionBegin;
     104         314 :   PetscCall(PetscBLASIntCast(ds->n*ctx->d,&nd));
     105         314 :   PetscCall(PetscBLASIntCast(ds->n,&n));
     106         314 :   PetscCall(PetscBLASIntCast(ds->ld,&ld));
     107         314 :   PetscCall(PetscBLASIntCast(ds->ld*ctx->d,&ldd));
     108         314 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_A));
     109         314 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_B));
     110         314 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_W));
     111         314 :   PetscCall(DSAllocateMat_Private(ds,DS_MAT_U));
     112         314 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_A],&A));
     113         314 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_B],&B));
     114             : 
     115             :   /* build matrices A and B of the linearization */
     116         314 :   PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[ctx->d]],&Ed));
     117         314 :   PetscCall(PetscArrayzero(A,ldd*ldd));
     118         314 :   if (!ctx->pbc) { /* monomial basis */
     119        3073 :     for (i=0;i<nd-ds->n;i++) A[i+(i+ds->n)*ldd] = 1.0;
     120         815 :     for (i=0;i<ctx->d;i++) {
     121         544 :       PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     122         544 :       off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     123        6130 :       for (j=0;j<ds->n;j++) PetscCall(PetscArraycpy(A+off+j*ldd,Ei+j*ds->ld,ds->n));
     124         544 :       PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     125             :     }
     126             :   } else {
     127          43 :     ca = ctx->pbc;
     128          43 :     cb = ca+ctx->d+1;
     129          43 :     cg = cb+ctx->d+1;
     130         776 :     for (i=0;i<ds->n;i++) {
     131         733 :       A[i+(i+ds->n)*ldd] = ca[0];
     132         733 :       A[i+i*ldd] = cb[0];
     133             :     }
     134         640 :     for (;i<nd-ds->n;i++) {
     135         597 :       j = i/ds->n;
     136         597 :       A[i+(i+ds->n)*ldd] = ca[j];
     137         597 :       A[i+i*ldd] = cb[j];
     138         597 :       A[i+(i-ds->n)*ldd] = cg[j];
     139             :     }
     140          86 :     for (i=0;i<ctx->d-2;i++) {
     141          43 :       PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     142          43 :       off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     143         640 :       for (j=0;j<ds->n;j++)
     144       16284 :         for (k=0;k<ds->n;k++)
     145       15687 :           A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1];
     146          43 :       PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     147             :     }
     148          43 :     PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     149          43 :     off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     150         776 :     for (j=0;j<ds->n;j++)
     151       32494 :       for (k=0;k<ds->n;k++)
     152       31761 :         A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1]-Ed[j*ds->ld+k]*cg[ctx->d-1];
     153          43 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     154          43 :     i++;
     155          43 :     PetscCall(MatDenseGetArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     156          43 :     off = i*ds->n*ldd+(ctx->d-1)*ds->n;
     157         776 :     for (j=0;j<ds->n;j++)
     158       32494 :       for (k=0;k<ds->n;k++)
     159       31761 :         A[off+j*ldd+k] = Ei[j*ds->ld+k]*ca[ctx->d-1]-Ed[j*ds->ld+k]*cb[ctx->d-1];
     160          43 :     PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[i]],&Ei));
     161             :   }
     162         314 :   PetscCall(PetscArrayzero(B,ldd*ldd));
     163        4446 :   for (i=0;i<nd-ds->n;i++) B[i+i*ldd] = 1.0;
     164         314 :   off = (ctx->d-1)*ds->n*(ldd+1);
     165        3831 :   for (j=0;j<ds->n;j++) {
     166       71658 :     for (i=0;i<ds->n;i++) B[off+i+j*ldd] = -Ed[i+j*ds->ld];
     167             :   }
     168         314 :   PetscCall(MatDenseRestoreArrayRead(ds->omat[DSMatExtra[ctx->d]],&Ed));
     169             : 
     170             :   /* solve generalized eigenproblem */
     171         314 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_W],&W));
     172         314 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_U],&U));
     173         314 :   lwork = -1;
     174             : #if defined(PETSC_USE_COMPLEX)
     175         314 :   if (useggev3) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,NULL,U,&ldd,W,&ldd,&a,&lwork,NULL,&info));
     176         314 :   else PetscCallBLAS("LAPACKggev",LAPACKggev_("V","V",&nd,A,&ldd,B,&ldd,wr,NULL,U,&ldd,W,&ldd,&a,&lwork,NULL,&info));
     177         314 :   PetscCall(PetscBLASIntCast((PetscInt)PetscRealPart(a),&lwork));
     178         314 :   PetscCall(DSAllocateWork_Private(ds,lwork+nd,8*nd,0));
     179         314 :   beta  = ds->work;
     180         314 :   work  = ds->work + nd;
     181         314 :   if (useggev3) PetscCallBLAS("LAPACKggev3",LAPACKggev3_("V","V",&nd,A,&ldd,B,&ldd,wr,beta,U,&ldd,W,&ldd,work,&lwork,ds->rwork,&info));
     182         314 :   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         314 :   SlepcCheckLapackInfo(useggev3?"ggev3":"ggev",info);
     194         314 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_A],&A));
     195         314 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_B],&B));
     196             : 
     197             :   /* copy eigenvalues */
     198        7963 :   for (i=0;i<nd;i++) {
     199        7649 :     if (beta[i]==0.0) wr[i] = (PetscRealPart(wr[i])>0.0)? PETSC_MAX_REAL: PETSC_MIN_REAL;
     200        7649 :     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        7649 :     if (wi) wi[i] = 0.0;
     206             : #endif
     207             :   }
     208             : 
     209             :   /* copy and normalize eigenvectors */
     210         314 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
     211         314 :   PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
     212        7963 :   for (j=0;j<nd;j++) {
     213        7649 :     PetscCall(PetscArraycpy(X+j*ds->ld,W+j*ldd,ds->n));
     214        7649 :     PetscCall(PetscArraycpy(Y+j*ds->ld,U+ds->n*(ctx->d-1)+j*ldd,ds->n));
     215             :   }
     216         314 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_W],&W));
     217         314 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_U],&U));
     218        7963 :   for (j=0;j<nd;j++) {
     219        7649 :     cols = 1;
     220        7649 :     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        7649 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,X+j*ds->ld,&ld,&info));
     228        7649 :     SlepcCheckLapackInfo("lascl",info);
     229        7649 :     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        7649 :     PetscCallBLAS("LAPACKlascl",LAPACKlascl_("G",&zero,&zero,&norm,&done,&n,&cols,Y+j*ds->ld,&ld,&info));
     234        7649 :     SlepcCheckLapackInfo("lascl",info);
     235             : #if !defined(PETSC_USE_COMPLEX)
     236             :     if (wi[j] != 0.0) j++;
     237             : #endif
     238             :   }
     239         314 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
     240         314 :   PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
     241         314 :   PetscFunctionReturn(PETSC_SUCCESS);
     242             : }
     243             : 
     244             : #if !defined(PETSC_HAVE_MPIUNI)
     245          54 : static PetscErrorCode DSSynchronize_PEP(DS ds,PetscScalar eigr[],PetscScalar eigi[])
     246             : {
     247          54 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     248          54 :   PetscInt       ld=ds->ld,k=0;
     249          54 :   PetscMPIInt    ldnd,rank,off=0,size,dn;
     250          54 :   PetscScalar    *X,*Y;
     251             : 
     252          54 :   PetscFunctionBegin;
     253          54 :   if (ds->state>=DS_STATE_CONDENSED) k += 2*ctx->d*ds->n*ld;
     254          54 :   if (eigr) k += ctx->d*ds->n;
     255          54 :   if (eigi) k += ctx->d*ds->n;
     256          54 :   PetscCall(DSAllocateWork_Private(ds,k,0,0));
     257          54 :   PetscCall(PetscMPIIntCast(k*sizeof(PetscScalar),&size));
     258          54 :   PetscCall(PetscMPIIntCast(ds->n*ctx->d*ld,&ldnd));
     259          54 :   PetscCall(PetscMPIIntCast(ctx->d*ds->n,&dn));
     260          54 :   if (ds->state>=DS_STATE_CONDENSED) {
     261          54 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_X],&X));
     262          54 :     PetscCall(MatDenseGetArray(ds->omat[DS_MAT_Y],&Y));
     263             :   }
     264          54 :   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)ds),&rank));
     265          54 :   if (!rank) {
     266          27 :     if (ds->state>=DS_STATE_CONDENSED) {
     267          27 :       PetscCallMPI(MPI_Pack(X,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     268          27 :       PetscCallMPI(MPI_Pack(Y,ldnd,MPIU_SCALAR,ds->work,size,&off,PetscObjectComm((PetscObject)ds)));
     269             :     }
     270          27 :     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         108 :   PetscCallMPI(MPI_Bcast(ds->work,size,MPI_BYTE,0,PetscObjectComm((PetscObject)ds)));
     276          54 :   if (rank) {
     277          27 :     if (ds->state>=DS_STATE_CONDENSED) {
     278          27 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,X,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     279          27 :       PetscCallMPI(MPI_Unpack(ds->work,size,&off,Y,ldnd,MPIU_SCALAR,PetscObjectComm((PetscObject)ds)));
     280             :     }
     281          27 :     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          54 :   if (ds->state>=DS_STATE_CONDENSED) {
     287          54 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_X],&X));
     288          54 :     PetscCall(MatDenseRestoreArray(ds->omat[DS_MAT_Y],&Y));
     289             :   }
     290          54 :   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          75 :   PetscValidLogicalCollectiveInt(ds,d,2);
     323          25 :   PetscTryMethod(ds,"DSPEPSetDegree_C",(DS,PetscInt),(ds,d));
     324          25 :   PetscFunctionReturn(PETSC_SUCCESS);
     325             : }
     326             : 
     327        1305 : static PetscErrorCode DSPEPGetDegree_PEP(DS ds,PetscInt *d)
     328             : {
     329        1305 :   DS_PEP *ctx = (DS_PEP*)ds->data;
     330             : 
     331        1305 :   PetscFunctionBegin;
     332        1305 :   *d = ctx->d;
     333        1305 :   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        1305 : PetscErrorCode DSPEPGetDegree(DS ds,PetscInt *d)
     352             : {
     353        1305 :   PetscFunctionBegin;
     354        1305 :   PetscValidHeaderSpecific(ds,DS_CLASSID,1);
     355        1305 :   PetscAssertPointer(d,2);
     356        1305 :   PetscUseMethod(ds,"DSPEPGetDegree_C",(DS,PetscInt*),(ds,d));
     357        1305 :   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             : /*@
     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 :   PetscAssertPointer(pbc,2);
     405          12 :   PetscTryMethod(ds,"DSPEPSetCoefficients_C",(DS,PetscReal*),(ds,pbc));
     406          12 :   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          24 : static PetscErrorCode DSDestroy_PEP(DS ds)
     457             : {
     458          24 :   DS_PEP         *ctx = (DS_PEP*)ds->data;
     459             : 
     460          24 :   PetscFunctionBegin;
     461          24 :   if (ctx->pbc) PetscCall(PetscFree(ctx->pbc));
     462          24 :   PetscCall(PetscFree(ds->data));
     463          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",NULL));
     464          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",NULL));
     465          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",NULL));
     466          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",NULL));
     467          24 :   PetscFunctionReturn(PETSC_SUCCESS);
     468             : }
     469             : 
     470        1027 : static PetscErrorCode DSMatGetSize_PEP(DS ds,DSMatType t,PetscInt *rows,PetscInt *cols)
     471             : {
     472        1027 :   DS_PEP *ctx = (DS_PEP*)ds->data;
     473             : 
     474        1027 :   PetscFunctionBegin;
     475        1027 :   PetscCheck(ctx->d,PetscObjectComm((PetscObject)ds),PETSC_ERR_ARG_WRONGSTATE,"DSPEP requires specifying the polynomial degree via DSPEPSetDegree()");
     476        1027 :   *rows = ds->n;
     477        1027 :   if (t==DS_MAT_A || t==DS_MAT_B || t==DS_MAT_W || t==DS_MAT_U) *rows *= ctx->d;
     478        1027 :   *cols = ds->n;
     479        1027 :   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        1027 :   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          24 : SLEPC_EXTERN PetscErrorCode DSCreate_PEP(DS ds)
     517             : {
     518          24 :   DS_PEP         *ctx;
     519             : 
     520          24 :   PetscFunctionBegin;
     521          24 :   PetscCall(PetscNew(&ctx));
     522          24 :   ds->data = (void*)ctx;
     523             : 
     524          24 :   ds->ops->allocate      = DSAllocate_PEP;
     525          24 :   ds->ops->view          = DSView_PEP;
     526          24 :   ds->ops->vectors       = DSVectors_PEP;
     527          24 :   ds->ops->solve[0]      = DSSolve_PEP_QZ;
     528             : #if !defined(SLEPC_MISSING_LAPACK_GGES3)
     529          24 :   ds->ops->solve[1]      = DSSolve_PEP_QZ;
     530             : #endif
     531          24 :   ds->ops->sort          = DSSort_PEP;
     532             : #if !defined(PETSC_HAVE_MPIUNI)
     533          24 :   ds->ops->synchronize   = DSSynchronize_PEP;
     534             : #endif
     535          24 :   ds->ops->destroy       = DSDestroy_PEP;
     536          24 :   ds->ops->matgetsize    = DSMatGetSize_PEP;
     537          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetDegree_C",DSPEPSetDegree_PEP));
     538          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetDegree_C",DSPEPGetDegree_PEP));
     539          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPSetCoefficients_C",DSPEPSetCoefficients_PEP));
     540          24 :   PetscCall(PetscObjectComposeFunction((PetscObject)ds,"DSPEPGetCoefficients_C",DSPEPGetCoefficients_PEP));
     541          24 :   PetscFunctionReturn(PETSC_SUCCESS);
     542             : }