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             :    Explicit linearization for polynomial eigenproblems
      12             : */
      13             : 
      14             : #include <slepc/private/pepimpl.h>         /*I "slepcpep.h" I*/
      15             : #include "linear.h"
      16             : 
      17         801 : static PetscErrorCode MatMult_Linear_Shift(Mat M,Vec x,Vec y)
      18             : {
      19         801 :   PEP_LINEAR        *ctx;
      20         801 :   PEP               pep;
      21         801 :   const PetscScalar *px;
      22         801 :   PetscScalar       *py,a,sigma=0.0;
      23         801 :   PetscInt          nmat,deg,i,m;
      24         801 :   Vec               x1,x2,x3,y1,aux;
      25         801 :   PetscReal         *ca,*cb,*cg;
      26         801 :   PetscBool         flg;
      27             : 
      28         801 :   PetscFunctionBegin;
      29         801 :   PetscCall(MatShellGetContext(M,&ctx));
      30         801 :   pep = ctx->pep;
      31         801 :   PetscCall(STGetTransform(pep->st,&flg));
      32         801 :   if (!flg) PetscCall(STGetShift(pep->st,&sigma));
      33         801 :   nmat = pep->nmat;
      34         801 :   deg = nmat-1;
      35         801 :   m = pep->nloc;
      36         801 :   ca = pep->pbc;
      37         801 :   cb = pep->pbc+nmat;
      38         801 :   cg = pep->pbc+2*nmat;
      39         801 :   x1=ctx->w[0];x2=ctx->w[1];x3=ctx->w[2];y1=ctx->w[3];aux=ctx->w[4];
      40             : 
      41         801 :   PetscCall(VecSet(y,0.0));
      42         801 :   PetscCall(VecGetArrayRead(x,&px));
      43         801 :   PetscCall(VecGetArray(y,&py));
      44         801 :   a = 1.0;
      45             : 
      46             :   /* first block */
      47         801 :   PetscCall(VecPlaceArray(x2,px));
      48         801 :   PetscCall(VecPlaceArray(x3,px+m));
      49         801 :   PetscCall(VecPlaceArray(y1,py));
      50         801 :   PetscCall(VecAXPY(y1,cb[0]-sigma,x2));
      51         801 :   PetscCall(VecAXPY(y1,ca[0],x3));
      52         801 :   PetscCall(VecResetArray(x2));
      53         801 :   PetscCall(VecResetArray(x3));
      54         801 :   PetscCall(VecResetArray(y1));
      55             : 
      56             :   /* inner blocks */
      57         893 :   for (i=1;i<deg-1;i++) {
      58          92 :     PetscCall(VecPlaceArray(x1,px+(i-1)*m));
      59          92 :     PetscCall(VecPlaceArray(x2,px+i*m));
      60          92 :     PetscCall(VecPlaceArray(x3,px+(i+1)*m));
      61          92 :     PetscCall(VecPlaceArray(y1,py+i*m));
      62          92 :     PetscCall(VecAXPY(y1,cg[i],x1));
      63          92 :     PetscCall(VecAXPY(y1,cb[i]-sigma,x2));
      64          92 :     PetscCall(VecAXPY(y1,ca[i],x3));
      65          92 :     PetscCall(VecResetArray(x1));
      66          92 :     PetscCall(VecResetArray(x2));
      67          92 :     PetscCall(VecResetArray(x3));
      68          92 :     PetscCall(VecResetArray(y1));
      69             :   }
      70             : 
      71             :   /* last block */
      72         801 :   PetscCall(VecPlaceArray(y1,py+(deg-1)*m));
      73        2495 :   for (i=0;i<deg;i++) {
      74        1694 :     PetscCall(VecPlaceArray(x1,px+i*m));
      75        1694 :     PetscCall(STMatMult(pep->st,i,x1,aux));
      76        1694 :     PetscCall(VecAXPY(y1,a,aux));
      77        1694 :     PetscCall(VecResetArray(x1));
      78        1694 :     a *= pep->sfactor;
      79             :   }
      80         801 :   PetscCall(VecCopy(y1,aux));
      81         801 :   PetscCall(STMatSolve(pep->st,aux,y1));
      82         801 :   PetscCall(VecScale(y1,-ca[deg-1]/a));
      83         801 :   PetscCall(VecPlaceArray(x1,px+(deg-2)*m));
      84         801 :   PetscCall(VecPlaceArray(x2,px+(deg-1)*m));
      85         801 :   PetscCall(VecAXPY(y1,cg[deg-1],x1));
      86         801 :   PetscCall(VecAXPY(y1,cb[deg-1]-sigma,x2));
      87         801 :   PetscCall(VecResetArray(x1));
      88         801 :   PetscCall(VecResetArray(x2));
      89         801 :   PetscCall(VecResetArray(y1));
      90             : 
      91         801 :   PetscCall(VecRestoreArrayRead(x,&px));
      92         801 :   PetscCall(VecRestoreArray(y,&py));
      93         801 :   PetscFunctionReturn(PETSC_SUCCESS);
      94             : }
      95             : 
      96        1284 : static PetscErrorCode MatMult_Linear_Sinvert(Mat M,Vec x,Vec y)
      97             : {
      98        1284 :   PEP_LINEAR        *ctx;
      99        1284 :   PEP               pep;
     100        1284 :   const PetscScalar *px;
     101        1284 :   PetscScalar       *py,a,sigma,t=1.0,tp=0.0,tt;
     102        1284 :   PetscInt          nmat,deg,i,m;
     103        1284 :   Vec               x1,y1,y2,y3,aux,aux2;
     104        1284 :   PetscReal         *ca,*cb,*cg;
     105             : 
     106        1284 :   PetscFunctionBegin;
     107        1284 :   PetscCall(MatShellGetContext(M,&ctx));
     108        1284 :   pep = ctx->pep;
     109        1284 :   nmat = pep->nmat;
     110        1284 :   deg = nmat-1;
     111        1284 :   m = pep->nloc;
     112        1284 :   ca = pep->pbc;
     113        1284 :   cb = pep->pbc+nmat;
     114        1284 :   cg = pep->pbc+2*nmat;
     115        1284 :   x1=ctx->w[0];y1=ctx->w[1];y2=ctx->w[2];y3=ctx->w[3];aux=ctx->w[4];aux2=ctx->w[5];
     116        1284 :   PetscCall(EPSGetTarget(ctx->eps,&sigma));
     117        1284 :   PetscCall(VecSet(y,0.0));
     118        1284 :   PetscCall(VecGetArrayRead(x,&px));
     119        1284 :   PetscCall(VecGetArray(y,&py));
     120        1284 :   a = pep->sfactor;
     121             : 
     122             :   /* first block */
     123        1284 :   PetscCall(VecPlaceArray(x1,px));
     124        1284 :   PetscCall(VecPlaceArray(y1,py+m));
     125        1284 :   PetscCall(VecCopy(x1,y1));
     126        1284 :   PetscCall(VecScale(y1,1.0/ca[0]));
     127        1284 :   PetscCall(VecResetArray(x1));
     128        1284 :   PetscCall(VecResetArray(y1));
     129             : 
     130             :   /* second block */
     131        1284 :   if (deg>2) {
     132         452 :     PetscCall(VecPlaceArray(x1,px+m));
     133         452 :     PetscCall(VecPlaceArray(y1,py+m));
     134         452 :     PetscCall(VecPlaceArray(y2,py+2*m));
     135         452 :     PetscCall(VecCopy(x1,y2));
     136         452 :     PetscCall(VecAXPY(y2,sigma-cb[1],y1));
     137         452 :     PetscCall(VecScale(y2,1.0/ca[1]));
     138         452 :     PetscCall(VecResetArray(x1));
     139         452 :     PetscCall(VecResetArray(y1));
     140         452 :     PetscCall(VecResetArray(y2));
     141             :   }
     142             : 
     143             :   /* inner blocks */
     144        1736 :   for (i=2;i<deg-1;i++) {
     145         452 :     PetscCall(VecPlaceArray(x1,px+i*m));
     146         452 :     PetscCall(VecPlaceArray(y1,py+(i-1)*m));
     147         452 :     PetscCall(VecPlaceArray(y2,py+i*m));
     148         452 :     PetscCall(VecPlaceArray(y3,py+(i+1)*m));
     149         452 :     PetscCall(VecCopy(x1,y3));
     150         452 :     PetscCall(VecAXPY(y3,sigma-cb[i],y2));
     151         452 :     PetscCall(VecAXPY(y3,-cg[i],y1));
     152         452 :     PetscCall(VecScale(y3,1.0/ca[i]));
     153         452 :     PetscCall(VecResetArray(x1));
     154         452 :     PetscCall(VecResetArray(y1));
     155         452 :     PetscCall(VecResetArray(y2));
     156         452 :     PetscCall(VecResetArray(y3));
     157             :   }
     158             : 
     159             :   /* last block */
     160        1284 :   PetscCall(VecPlaceArray(y1,py));
     161        2188 :   for (i=0;i<deg-2;i++) {
     162         904 :     PetscCall(VecPlaceArray(y2,py+(i+1)*m));
     163         904 :     PetscCall(STMatMult(pep->st,i+1,y2,aux));
     164         904 :     PetscCall(VecAXPY(y1,a,aux));
     165         904 :     PetscCall(VecResetArray(y2));
     166         904 :     a *= pep->sfactor;
     167             :   }
     168        1284 :   i = deg-2;
     169        1284 :   PetscCall(VecPlaceArray(y2,py+(i+1)*m));
     170        1284 :   PetscCall(VecPlaceArray(y3,py+i*m));
     171        1284 :   PetscCall(VecCopy(y2,aux2));
     172        1284 :   PetscCall(VecAXPY(aux2,cg[i+1]/ca[i+1],y3));
     173        1284 :   PetscCall(STMatMult(pep->st,i+1,aux2,aux));
     174        1284 :   PetscCall(VecAXPY(y1,a,aux));
     175        1284 :   PetscCall(VecResetArray(y2));
     176        1284 :   PetscCall(VecResetArray(y3));
     177        1284 :   a *= pep->sfactor;
     178        1284 :   i = deg-1;
     179        1284 :   PetscCall(VecPlaceArray(x1,px+i*m));
     180        1284 :   PetscCall(VecPlaceArray(y3,py+i*m));
     181        1284 :   PetscCall(VecCopy(x1,aux2));
     182        1284 :   PetscCall(VecAXPY(aux2,sigma-cb[i],y3));
     183        1284 :   PetscCall(VecScale(aux2,1.0/ca[i]));
     184        1284 :   PetscCall(STMatMult(pep->st,i+1,aux2,aux));
     185        1284 :   PetscCall(VecAXPY(y1,a,aux));
     186        1284 :   PetscCall(VecResetArray(x1));
     187        1284 :   PetscCall(VecResetArray(y3));
     188             : 
     189        1284 :   PetscCall(VecCopy(y1,aux));
     190        1284 :   PetscCall(STMatSolve(pep->st,aux,y1));
     191        1284 :   PetscCall(VecScale(y1,-1.0));
     192             : 
     193             :   /* final update */
     194        3472 :   for (i=1;i<deg;i++) {
     195        2188 :     PetscCall(VecPlaceArray(y2,py+i*m));
     196        2188 :     tt = t;
     197        2188 :     t = ((sigma-cb[i-1])*t-cg[i-1]*tp)/ca[i-1]; /* i-th basis polynomial */
     198        2188 :     tp = tt;
     199        2188 :     PetscCall(VecAXPY(y2,t,y1));
     200        2188 :     PetscCall(VecResetArray(y2));
     201             :   }
     202        1284 :   PetscCall(VecResetArray(y1));
     203             : 
     204        1284 :   PetscCall(VecRestoreArrayRead(x,&px));
     205        1284 :   PetscCall(VecRestoreArray(y,&py));
     206        1284 :   PetscFunctionReturn(PETSC_SUCCESS);
     207             : }
     208             : 
     209       18080 : static PetscErrorCode BackTransform_Linear(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
     210             : {
     211       18080 :   PEP_LINEAR     *ctx;
     212       18080 :   ST             stctx;
     213             : 
     214       18080 :   PetscFunctionBegin;
     215       18080 :   PetscCall(STShellGetContext(st,&ctx));
     216       18080 :   PetscCall(PEPGetST(ctx->pep,&stctx));
     217       18080 :   PetscCall(STBackTransform(stctx,n,eigr,eigi));
     218       18080 :   PetscFunctionReturn(PETSC_SUCCESS);
     219             : }
     220             : 
     221             : /*
     222             :    Dummy backtransform operation
     223             :  */
     224          46 : static PetscErrorCode BackTransform_Skip(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
     225             : {
     226          46 :   PetscFunctionBegin;
     227          46 :   PetscFunctionReturn(PETSC_SUCCESS);
     228             : }
     229             : 
     230        2085 : static PetscErrorCode Apply_Linear(ST st,Vec x,Vec y)
     231             : {
     232        2085 :   PEP_LINEAR     *ctx;
     233             : 
     234        2085 :   PetscFunctionBegin;
     235        2085 :   PetscCall(STShellGetContext(st,&ctx));
     236        2085 :   PetscCall(MatMult(ctx->A,x,y));
     237        2085 :   PetscFunctionReturn(PETSC_SUCCESS);
     238             : }
     239             : 
     240          34 : static PetscErrorCode PEPSetUp_Linear(PEP pep)
     241             : {
     242          34 :   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
     243          34 :   ST             st;
     244          34 :   PetscInt       i=0,deg=pep->nmat-1;
     245          34 :   EPSWhich       which = EPS_LARGEST_MAGNITUDE;
     246          34 :   EPSProblemType ptype;
     247          34 :   PetscBool      trackall,istrivial,transf,sinv,ks;
     248          34 :   PetscScalar    sigma,*epsarray,*peparray;
     249          34 :   Vec            veps,w=NULL;
     250             :   /* function tables */
     251          34 :   PetscErrorCode (*fcreate[][2])(MPI_Comm,PEP_LINEAR*,Mat*) = {
     252             :     { MatCreateExplicit_Linear_NA, MatCreateExplicit_Linear_NB },
     253             :     { MatCreateExplicit_Linear_SA, MatCreateExplicit_Linear_SB },
     254             :     { MatCreateExplicit_Linear_HA, MatCreateExplicit_Linear_HB },
     255             :   };
     256             : 
     257          34 :   PetscFunctionBegin;
     258          34 :   PEPCheckShiftSinvert(pep);
     259          34 :   PEPCheckUnsupported(pep,PEP_FEATURE_STOPPING);
     260          34 :   PEPCheckIgnored(pep,PEP_FEATURE_CONVERGENCE);
     261          34 :   PetscCall(STGetTransform(pep->st,&transf));
     262          34 :   PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
     263          34 :   if (!pep->which) PetscCall(PEPSetWhichEigenpairs_Default(pep));
     264          34 :   PetscCheck(pep->which!=PEP_ALL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues");
     265          34 :   PetscCall(STSetUp(pep->st));
     266          34 :   if (!ctx->eps) PetscCall(PEPLinearGetEPS(pep,&ctx->eps));
     267          34 :   PetscCall(EPSGetST(ctx->eps,&st));
     268          34 :   if (!transf && !ctx->usereps) PetscCall(EPSSetTarget(ctx->eps,pep->target));
     269          34 :   if (sinv && !transf && !ctx->usereps) PetscCall(STSetDefaultShift(st,pep->target));
     270             :   /* compute scale factor if not set by user */
     271          34 :   PetscCall(PEPComputeScaleFactor(pep));
     272             : 
     273          34 :   if (ctx->explicitmatrix) {
     274          11 :     PEPCheckQuadraticCondition(pep,PETSC_TRUE," (with explicit matrix)");
     275          11 :     PEPCheckUnsupportedCondition(pep,PEP_FEATURE_NONMONOMIAL,PETSC_TRUE," (with explicit matrix)");
     276          11 :     PetscCheck(!transf,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Explicit matrix option is not implemented with st-transform flag active");
     277          11 :     PetscCheck(pep->scale!=PEP_SCALE_DIAGONAL && pep->scale!=PEP_SCALE_BOTH,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Diagonal scaling not allowed in PEPLINEAR with explicit matrices");
     278          11 :     if (sinv && !transf) PetscCall(STSetType(st,STSINVERT));
     279          11 :     PetscCall(RGPushScale(pep->rg,1.0/pep->sfactor));
     280          11 :     PetscCall(STGetMatrixTransformed(pep->st,0,&ctx->K));
     281          11 :     PetscCall(STGetMatrixTransformed(pep->st,1,&ctx->C));
     282          11 :     PetscCall(STGetMatrixTransformed(pep->st,2,&ctx->M));
     283          11 :     ctx->sfactor = pep->sfactor;
     284          11 :     ctx->dsfactor = pep->dsfactor;
     285             : 
     286          11 :     PetscCall(MatDestroy(&ctx->A));
     287          11 :     PetscCall(MatDestroy(&ctx->B));
     288          11 :     PetscCall(VecDestroy(&ctx->w[0]));
     289          11 :     PetscCall(VecDestroy(&ctx->w[1]));
     290          11 :     PetscCall(VecDestroy(&ctx->w[2]));
     291          11 :     PetscCall(VecDestroy(&ctx->w[3]));
     292             : 
     293          11 :     switch (pep->problem_type) {
     294             :       case PEP_GENERAL:    i = 0; break;
     295           3 :       case PEP_HERMITIAN:
     296           3 :       case PEP_HYPERBOLIC: i = 1; break;
     297           2 :       case PEP_GYROSCOPIC: i = 2; break;
     298             :     }
     299             : 
     300          11 :     PetscCall((*fcreate[i][0])(PetscObjectComm((PetscObject)pep),ctx,&ctx->A));
     301          11 :     PetscCall((*fcreate[i][1])(PetscObjectComm((PetscObject)pep),ctx,&ctx->B));
     302             : 
     303             :   } else {   /* implicit matrix */
     304          23 :     PetscCheck(pep->problem_type==PEP_GENERAL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Must use the explicit matrix option if problem type is not general");
     305          23 :     if (!((PetscObject)ctx->eps)->type_name) PetscCall(EPSSetType(ctx->eps,EPSKRYLOVSCHUR));
     306             :     else {
     307          22 :       PetscCall(PetscObjectTypeCompare((PetscObject)ctx->eps,EPSKRYLOVSCHUR,&ks));
     308          22 :       PetscCheck(ks,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option only implemented for Krylov-Schur");
     309             :     }
     310          23 :     PetscCheck(ctx->alpha==1.0 && ctx->beta==0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Implicit matrix option does not support setting alpha,beta parameters of the linearization");
     311          23 :     PetscCall(STSetType(st,STSHELL));
     312          23 :     PetscCall(STShellSetContext(st,ctx));
     313          23 :     if (!transf) PetscCall(STShellSetBackTransform(st,BackTransform_Linear));
     314           1 :     else PetscCall(STShellSetBackTransform(st,BackTransform_Skip));
     315          23 :     PetscCall(MatCreateVecsEmpty(pep->A[0],&ctx->w[0],&ctx->w[1]));
     316          23 :     PetscCall(MatCreateVecsEmpty(pep->A[0],&ctx->w[2],&ctx->w[3]));
     317          23 :     PetscCall(MatCreateVecs(pep->A[0],&ctx->w[4],&ctx->w[5]));
     318          23 :     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pep),deg*pep->nloc,deg*pep->nloc,deg*pep->n,deg*pep->n,ctx,&ctx->A));
     319          23 :     if (sinv && !transf) PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Sinvert));
     320          13 :     else PetscCall(MatShellSetOperation(ctx->A,MATOP_MULT,(void(*)(void))MatMult_Linear_Shift));
     321          23 :     PetscCall(STShellSetApply(st,Apply_Linear));
     322          23 :     ctx->pep = pep;
     323             : 
     324          23 :     PetscCall(PEPBasisCoefficients(pep,pep->pbc));
     325          23 :     if (!transf) {
     326          22 :       PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
     327          22 :       if (sinv) PetscCall(PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL));
     328             :       else {
     329          38 :         for (i=0;i<deg;i++) pep->solvematcoeffs[i] = 0.0;
     330          12 :         pep->solvematcoeffs[deg] = 1.0;
     331             :       }
     332          22 :       PetscCall(STScaleShift(pep->st,1.0/pep->sfactor));
     333          22 :       PetscCall(RGPushScale(pep->rg,1.0/pep->sfactor));
     334             :     }
     335          23 :     if (pep->sfactor!=1.0) {
     336           4 :       for (i=0;i<pep->nmat;i++) {
     337           3 :         pep->pbc[pep->nmat+i] /= pep->sfactor;
     338           3 :         pep->pbc[2*pep->nmat+i] /= pep->sfactor*pep->sfactor;
     339             :       }
     340             :     }
     341             :   }
     342             : 
     343          34 :   PetscCall(EPSSetOperators(ctx->eps,ctx->A,ctx->B));
     344          34 :   PetscCall(EPSGetProblemType(ctx->eps,&ptype));
     345          34 :   if (!ptype) {
     346          32 :     if (ctx->explicitmatrix) PetscCall(EPSSetProblemType(ctx->eps,EPS_GNHEP));
     347          22 :     else PetscCall(EPSSetProblemType(ctx->eps,EPS_NHEP));
     348             :   }
     349          34 :   if (!ctx->usereps) {
     350          33 :     if (transf) which = EPS_LARGEST_MAGNITUDE;
     351             :     else {
     352          32 :       switch (pep->which) {
     353             :         case PEP_LARGEST_MAGNITUDE:  which = EPS_LARGEST_MAGNITUDE; break;
     354           0 :         case PEP_SMALLEST_MAGNITUDE: which = EPS_SMALLEST_MAGNITUDE; break;
     355           0 :         case PEP_LARGEST_REAL:       which = EPS_LARGEST_REAL; break;
     356           0 :         case PEP_SMALLEST_REAL:      which = EPS_SMALLEST_REAL; break;
     357           0 :         case PEP_LARGEST_IMAGINARY:  which = EPS_LARGEST_IMAGINARY; break;
     358           0 :         case PEP_SMALLEST_IMAGINARY: which = EPS_SMALLEST_IMAGINARY; break;
     359          14 :         case PEP_TARGET_MAGNITUDE:   which = EPS_TARGET_MAGNITUDE; break;
     360           0 :         case PEP_TARGET_REAL:        which = EPS_TARGET_REAL; break;
     361           0 :         case PEP_TARGET_IMAGINARY:   which = EPS_TARGET_IMAGINARY; break;
     362           0 :         case PEP_ALL:                which = EPS_ALL; break;
     363           0 :         case PEP_WHICH_USER:         which = EPS_WHICH_USER;
     364           0 :           PetscCall(EPSSetEigenvalueComparison(ctx->eps,pep->sc->comparison,pep->sc->comparisonctx));
     365             :           break;
     366             :       }
     367          33 :     }
     368          33 :     PetscCall(EPSSetWhichEigenpairs(ctx->eps,which));
     369             : 
     370          33 :     PetscCall(EPSSetDimensions(ctx->eps,pep->nev,pep->ncv,pep->mpd));
     371          47 :     PetscCall(EPSSetTolerances(ctx->eps,SlepcDefaultTol(pep->tol),pep->max_it));
     372             :   }
     373          34 :   PetscCall(RGIsTrivial(pep->rg,&istrivial));
     374          34 :   if (!istrivial) {
     375           1 :     PetscCheck(!transf,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"PEPLINEAR does not support a nontrivial region with st-transform");
     376           1 :     PetscCall(EPSSetRG(ctx->eps,pep->rg));
     377             :   }
     378             :   /* Transfer the trackall option from pep to eps */
     379          34 :   PetscCall(PEPGetTrackAll(pep,&trackall));
     380          34 :   PetscCall(EPSSetTrackAll(ctx->eps,trackall));
     381             : 
     382             :   /* temporary change of target */
     383          34 :   if (pep->sfactor!=1.0) {
     384           2 :     PetscCall(EPSGetTarget(ctx->eps,&sigma));
     385           2 :     PetscCall(EPSSetTarget(ctx->eps,sigma/pep->sfactor));
     386             :   }
     387             : 
     388             :   /* process initial vector */
     389          34 :   if (pep->nini<0) {
     390           2 :     PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)ctx->eps),deg*pep->nloc,deg*pep->n,&veps));
     391           2 :     PetscCall(VecGetArray(veps,&epsarray));
     392           6 :     for (i=0;i<deg;i++) {
     393           4 :       if (i<-pep->nini) {
     394           4 :         PetscCall(VecGetArray(pep->IS[i],&peparray));
     395           4 :         PetscCall(PetscArraycpy(epsarray+i*pep->nloc,peparray,pep->nloc));
     396           4 :         PetscCall(VecRestoreArray(pep->IS[i],&peparray));
     397             :       } else {
     398           0 :         if (!w) PetscCall(VecDuplicate(pep->IS[0],&w));
     399           0 :         PetscCall(VecSetRandom(w,NULL));
     400           0 :         PetscCall(VecGetArray(w,&peparray));
     401           0 :         PetscCall(PetscArraycpy(epsarray+i*pep->nloc,peparray,pep->nloc));
     402           4 :         PetscCall(VecRestoreArray(w,&peparray));
     403             :       }
     404             :     }
     405           2 :     PetscCall(VecRestoreArray(veps,&epsarray));
     406           2 :     PetscCall(EPSSetInitialSpace(ctx->eps,1,&veps));
     407           2 :     PetscCall(VecDestroy(&veps));
     408           2 :     PetscCall(VecDestroy(&w));
     409           2 :     PetscCall(SlepcBasisDestroy_Private(&pep->nini,&pep->IS));
     410             :   }
     411             : 
     412          34 :   PetscCall(EPSSetUp(ctx->eps));
     413          34 :   PetscCall(EPSGetDimensions(ctx->eps,NULL,&pep->ncv,&pep->mpd));
     414          34 :   PetscCall(EPSGetTolerances(ctx->eps,NULL,&pep->max_it));
     415          34 :   PetscCall(PEPAllocateSolution(pep,0));
     416          34 :   PetscFunctionReturn(PETSC_SUCCESS);
     417             : }
     418             : 
     419             : /*
     420             :    PEPLinearExtract_Residual - Auxiliary routine that copies the solution of the
     421             :    linear eigenproblem to the PEP object. The eigenvector of the generalized
     422             :    problem is supposed to be
     423             :                                z = [  x  ]
     424             :                                    [ l*x ]
     425             :    The eigenvector is taken from z(1:n) or z(n+1:2*n) depending on the explicitly
     426             :    computed residual norm.
     427             :    Finally, x is normalized so that ||x||_2 = 1.
     428             : */
     429           2 : static PetscErrorCode PEPLinearExtract_Residual(PEP pep,EPS eps)
     430             : {
     431           2 :   PetscInt          i,k;
     432           2 :   const PetscScalar *px;
     433           2 :   PetscScalar       *er=pep->eigr,*ei=pep->eigi;
     434           2 :   PetscReal         rn1,rn2;
     435           2 :   Vec               xr,xi=NULL,wr;
     436           2 :   Mat               A;
     437             : #if !defined(PETSC_USE_COMPLEX)
     438           2 :   Vec               wi;
     439           2 :   const PetscScalar *py;
     440             : #endif
     441             : 
     442           2 :   PetscFunctionBegin;
     443             : #if defined(PETSC_USE_COMPLEX)
     444             :   PetscCall(PEPSetWorkVecs(pep,2));
     445             : #else
     446           2 :   PetscCall(PEPSetWorkVecs(pep,4));
     447             : #endif
     448           2 :   PetscCall(EPSGetOperators(eps,&A,NULL));
     449           2 :   PetscCall(MatCreateVecs(A,&xr,NULL));
     450           2 :   PetscCall(MatCreateVecsEmpty(pep->A[0],&wr,NULL));
     451             : #if !defined(PETSC_USE_COMPLEX)
     452           2 :   PetscCall(VecDuplicate(xr,&xi));
     453           2 :   PetscCall(VecDuplicateEmpty(wr,&wi));
     454             : #endif
     455          18 :   for (i=0;i<pep->nconv;i++) {
     456          16 :     PetscCall(EPSGetEigenpair(eps,i,NULL,NULL,xr,xi));
     457             : #if !defined(PETSC_USE_COMPLEX)
     458          16 :     if (ei[i]!=0.0) {   /* complex conjugate pair */
     459           2 :       PetscCall(VecGetArrayRead(xr,&px));
     460           2 :       PetscCall(VecGetArrayRead(xi,&py));
     461           2 :       PetscCall(VecPlaceArray(wr,px));
     462           2 :       PetscCall(VecPlaceArray(wi,py));
     463           2 :       PetscCall(VecNormalizeComplex(wr,wi,PETSC_TRUE,NULL));
     464           2 :       PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,wi,pep->work,&rn1));
     465           2 :       PetscCall(BVInsertVec(pep->V,i,wr));
     466           2 :       PetscCall(BVInsertVec(pep->V,i+1,wi));
     467           4 :       for (k=1;k<pep->nmat-1;k++) {
     468           2 :         PetscCall(VecResetArray(wr));
     469           2 :         PetscCall(VecResetArray(wi));
     470           2 :         PetscCall(VecPlaceArray(wr,px+k*pep->nloc));
     471           2 :         PetscCall(VecPlaceArray(wi,py+k*pep->nloc));
     472           2 :         PetscCall(VecNormalizeComplex(wr,wi,PETSC_TRUE,NULL));
     473           2 :         PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,wi,pep->work,&rn2));
     474           2 :         if (rn1>rn2) {
     475           0 :           PetscCall(BVInsertVec(pep->V,i,wr));
     476           0 :           PetscCall(BVInsertVec(pep->V,i+1,wi));
     477           0 :           rn1 = rn2;
     478             :         }
     479             :       }
     480           2 :       PetscCall(VecResetArray(wr));
     481           2 :       PetscCall(VecResetArray(wi));
     482           2 :       PetscCall(VecRestoreArrayRead(xr,&px));
     483           2 :       PetscCall(VecRestoreArrayRead(xi,&py));
     484             :       i++;
     485             :     } else   /* real eigenvalue */
     486             : #endif
     487             :     {
     488          14 :       PetscCall(VecGetArrayRead(xr,&px));
     489          14 :       PetscCall(VecPlaceArray(wr,px));
     490          14 :       PetscCall(VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL));
     491          14 :       PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn1));
     492          14 :       PetscCall(BVInsertVec(pep->V,i,wr));
     493          28 :       for (k=1;k<pep->nmat-1;k++) {
     494          14 :         PetscCall(VecResetArray(wr));
     495          14 :         PetscCall(VecPlaceArray(wr,px+k*pep->nloc));
     496          14 :         PetscCall(VecNormalizeComplex(wr,NULL,PETSC_FALSE,NULL));
     497          14 :         PetscCall(PEPComputeResidualNorm_Private(pep,er[i],ei[i],wr,NULL,pep->work,&rn2));
     498          14 :         if (rn1>rn2) {
     499          14 :           PetscCall(BVInsertVec(pep->V,i,wr));
     500          14 :           rn1 = rn2;
     501             :         }
     502             :       }
     503          14 :       PetscCall(VecResetArray(wr));
     504          16 :       PetscCall(VecRestoreArrayRead(xr,&px));
     505             :     }
     506             :   }
     507           2 :   PetscCall(VecDestroy(&wr));
     508           2 :   PetscCall(VecDestroy(&xr));
     509             : #if !defined(PETSC_USE_COMPLEX)
     510           2 :   PetscCall(VecDestroy(&wi));
     511           2 :   PetscCall(VecDestroy(&xi));
     512             : #endif
     513           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     514             : }
     515             : 
     516             : /*
     517             :    PEPLinearExtract_None - Same as PEPLinearExtract_Norm but always takes
     518             :    the first block.
     519             : */
     520           2 : static PetscErrorCode PEPLinearExtract_None(PEP pep,EPS eps)
     521             : {
     522           2 :   PetscInt          i;
     523           2 :   const PetscScalar *px;
     524           2 :   Mat               A;
     525           2 :   Vec               xr,xi=NULL,w;
     526             : #if !defined(PETSC_USE_COMPLEX)
     527           2 :   PetscScalar       *ei=pep->eigi;
     528             : #endif
     529             : 
     530           2 :   PetscFunctionBegin;
     531           2 :   PetscCall(EPSGetOperators(eps,&A,NULL));
     532           2 :   PetscCall(MatCreateVecs(A,&xr,NULL));
     533             : #if !defined(PETSC_USE_COMPLEX)
     534           2 :   PetscCall(VecDuplicate(xr,&xi));
     535             : #endif
     536           2 :   PetscCall(MatCreateVecsEmpty(pep->A[0],&w,NULL));
     537          18 :   for (i=0;i<pep->nconv;i++) {
     538          16 :     PetscCall(EPSGetEigenvector(eps,i,xr,xi));
     539             : #if !defined(PETSC_USE_COMPLEX)
     540          16 :     if (ei[i]!=0.0) {   /* complex conjugate pair */
     541           2 :       PetscCall(VecGetArrayRead(xr,&px));
     542           2 :       PetscCall(VecPlaceArray(w,px));
     543           2 :       PetscCall(BVInsertVec(pep->V,i,w));
     544           2 :       PetscCall(VecResetArray(w));
     545           2 :       PetscCall(VecRestoreArrayRead(xr,&px));
     546           2 :       PetscCall(VecGetArrayRead(xi,&px));
     547           2 :       PetscCall(VecPlaceArray(w,px));
     548           2 :       PetscCall(BVInsertVec(pep->V,i+1,w));
     549           2 :       PetscCall(VecResetArray(w));
     550           2 :       PetscCall(VecRestoreArrayRead(xi,&px));
     551             :       i++;
     552             :     } else   /* real eigenvalue */
     553             : #endif
     554             :     {
     555          14 :       PetscCall(VecGetArrayRead(xr,&px));
     556          14 :       PetscCall(VecPlaceArray(w,px));
     557          14 :       PetscCall(BVInsertVec(pep->V,i,w));
     558          14 :       PetscCall(VecResetArray(w));
     559          16 :       PetscCall(VecRestoreArrayRead(xr,&px));
     560             :     }
     561             :   }
     562           2 :   PetscCall(VecDestroy(&w));
     563           2 :   PetscCall(VecDestroy(&xr));
     564             : #if !defined(PETSC_USE_COMPLEX)
     565           2 :   PetscCall(VecDestroy(&xi));
     566             : #endif
     567           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     568             : }
     569             : 
     570             : /*
     571             :    PEPLinearExtract_Norm - Auxiliary routine that copies the solution of the
     572             :    linear eigenproblem to the PEP object. The eigenvector of the generalized
     573             :    problem is supposed to be
     574             :                                z = [  x  ]
     575             :                                    [ l*x ]
     576             :    If |l|<1.0, the eigenvector is taken from z(1:n), otherwise from z(n+1:2*n).
     577             :    Finally, x is normalized so that ||x||_2 = 1.
     578             : */
     579          30 : static PetscErrorCode PEPLinearExtract_Norm(PEP pep,EPS eps)
     580             : {
     581          30 :   PetscInt          i,offset;
     582          30 :   const PetscScalar *px;
     583          30 :   PetscScalar       *er=pep->eigr;
     584          30 :   Mat               A;
     585          30 :   Vec               xr,xi=NULL,w;
     586             : #if !defined(PETSC_USE_COMPLEX)
     587          30 :   PetscScalar       *ei=pep->eigi;
     588             : #endif
     589             : 
     590          30 :   PetscFunctionBegin;
     591          30 :   PetscCall(EPSGetOperators(eps,&A,NULL));
     592          30 :   PetscCall(MatCreateVecs(A,&xr,NULL));
     593             : #if !defined(PETSC_USE_COMPLEX)
     594          30 :   PetscCall(VecDuplicate(xr,&xi));
     595             : #endif
     596          30 :   PetscCall(MatCreateVecsEmpty(pep->A[0],&w,NULL));
     597         160 :   for (i=0;i<pep->nconv;i++) {
     598         130 :     PetscCall(EPSGetEigenpair(eps,i,NULL,NULL,xr,xi));
     599         130 :     if (SlepcAbsEigenvalue(er[i],ei[i])>1.0) offset = (pep->nmat-2)*pep->nloc;
     600             :     else offset = 0;
     601             : #if !defined(PETSC_USE_COMPLEX)
     602         130 :     if (ei[i]!=0.0) {   /* complex conjugate pair */
     603          36 :       PetscCall(VecGetArrayRead(xr,&px));
     604          36 :       PetscCall(VecPlaceArray(w,px+offset));
     605          36 :       PetscCall(BVInsertVec(pep->V,i,w));
     606          36 :       PetscCall(VecResetArray(w));
     607          36 :       PetscCall(VecRestoreArrayRead(xr,&px));
     608          36 :       PetscCall(VecGetArrayRead(xi,&px));
     609          36 :       PetscCall(VecPlaceArray(w,px+offset));
     610          36 :       PetscCall(BVInsertVec(pep->V,i+1,w));
     611          36 :       PetscCall(VecResetArray(w));
     612          36 :       PetscCall(VecRestoreArrayRead(xi,&px));
     613             :       i++;
     614             :     } else /* real eigenvalue */
     615             : #endif
     616             :     {
     617          94 :       PetscCall(VecGetArrayRead(xr,&px));
     618          94 :       PetscCall(VecPlaceArray(w,px+offset));
     619          94 :       PetscCall(BVInsertVec(pep->V,i,w));
     620          94 :       PetscCall(VecResetArray(w));
     621         130 :       PetscCall(VecRestoreArrayRead(xr,&px));
     622             :     }
     623             :   }
     624          30 :   PetscCall(VecDestroy(&w));
     625          30 :   PetscCall(VecDestroy(&xr));
     626             : #if !defined(PETSC_USE_COMPLEX)
     627          30 :   PetscCall(VecDestroy(&xi));
     628             : #endif
     629          30 :   PetscFunctionReturn(PETSC_SUCCESS);
     630             : }
     631             : 
     632          34 : static PetscErrorCode PEPExtractVectors_Linear(PEP pep)
     633             : {
     634          34 :   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
     635             : 
     636          34 :   PetscFunctionBegin;
     637          34 :   switch (pep->extract) {
     638           2 :   case PEP_EXTRACT_NONE:
     639           2 :     PetscCall(PEPLinearExtract_None(pep,ctx->eps));
     640             :     break;
     641          30 :   case PEP_EXTRACT_NORM:
     642          30 :     PetscCall(PEPLinearExtract_Norm(pep,ctx->eps));
     643             :     break;
     644           2 :   case PEP_EXTRACT_RESIDUAL:
     645           2 :     PetscCall(PEPLinearExtract_Residual(pep,ctx->eps));
     646             :     break;
     647           0 :   case PEP_EXTRACT_STRUCTURED:
     648           0 :     SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Extraction not implemented in this solver");
     649             :   }
     650          34 :   PetscFunctionReturn(PETSC_SUCCESS);
     651             : }
     652             : 
     653          34 : static PetscErrorCode PEPSolve_Linear(PEP pep)
     654             : {
     655          34 :   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
     656          34 :   PetscScalar    sigma;
     657          34 :   PetscBool      flg;
     658          34 :   PetscInt       i;
     659             : 
     660          34 :   PetscFunctionBegin;
     661          34 :   PetscCall(EPSSolve(ctx->eps));
     662          34 :   PetscCall(EPSGetConverged(ctx->eps,&pep->nconv));
     663          34 :   PetscCall(EPSGetIterationNumber(ctx->eps,&pep->its));
     664          34 :   PetscCall(EPSGetConvergedReason(ctx->eps,(EPSConvergedReason*)&pep->reason));
     665             : 
     666             :   /* recover eigenvalues */
     667         236 :   for (i=0;i<pep->nconv;i++) {
     668         202 :     PetscCall(EPSGetEigenpair(ctx->eps,i,&pep->eigr[i],&pep->eigi[i],NULL,NULL));
     669         202 :     pep->eigr[i] *= pep->sfactor;
     670         202 :     pep->eigi[i] *= pep->sfactor;
     671             :   }
     672             : 
     673             :   /* restore target */
     674          34 :   PetscCall(EPSGetTarget(ctx->eps,&sigma));
     675          34 :   PetscCall(EPSSetTarget(ctx->eps,sigma*pep->sfactor));
     676             : 
     677          34 :   PetscCall(STGetTransform(pep->st,&flg));
     678          34 :   if (flg) PetscTryTypeMethod(pep,backtransform);
     679          34 :   if (pep->sfactor!=1.0) {
     680             :     /* Restore original values */
     681           8 :     for (i=0;i<pep->nmat;i++) {
     682           6 :       pep->pbc[pep->nmat+i] *= pep->sfactor;
     683           6 :       pep->pbc[2*pep->nmat+i] *= pep->sfactor*pep->sfactor;
     684             :     }
     685           2 :     if (!flg && !ctx->explicitmatrix) PetscCall(STScaleShift(pep->st,pep->sfactor));
     686             :   }
     687          34 :   if (ctx->explicitmatrix || !flg) PetscCall(RGPopScale(pep->rg));
     688          34 :   PetscFunctionReturn(PETSC_SUCCESS);
     689             : }
     690             : 
     691         413 : static PetscErrorCode EPSMonitor_Linear(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
     692             : {
     693         413 :   PEP            pep = (PEP)ctx;
     694             : 
     695         413 :   PetscFunctionBegin;
     696         413 :   PetscCall(PEPMonitor(pep,its,nconv,eigr,eigi,errest,nest));
     697         413 :   PetscFunctionReturn(PETSC_SUCCESS);
     698             : }
     699             : 
     700          31 : static PetscErrorCode PEPSetFromOptions_Linear(PEP pep,PetscOptionItems *PetscOptionsObject)
     701             : {
     702          31 :   PetscBool      set,val;
     703          31 :   PetscInt       k;
     704          31 :   PetscReal      array[2]={0,0};
     705          31 :   PetscBool      flg;
     706          31 :   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
     707             : 
     708          31 :   PetscFunctionBegin;
     709          31 :   PetscOptionsHeadBegin(PetscOptionsObject,"PEP Linear Options");
     710             : 
     711          31 :     k = 2;
     712          31 :     PetscCall(PetscOptionsRealArray("-pep_linear_linearization","Parameters of the linearization","PEPLinearSetLinearization",array,&k,&flg));
     713          31 :     if (flg) PetscCall(PEPLinearSetLinearization(pep,array[0],array[1]));
     714             : 
     715          31 :     PetscCall(PetscOptionsBool("-pep_linear_explicitmatrix","Use explicit matrix in linearization","PEPLinearSetExplicitMatrix",ctx->explicitmatrix,&val,&set));
     716          31 :     if (set) PetscCall(PEPLinearSetExplicitMatrix(pep,val));
     717             : 
     718          31 :   PetscOptionsHeadEnd();
     719             : 
     720          31 :   if (!ctx->eps) PetscCall(PEPLinearGetEPS(pep,&ctx->eps));
     721          31 :   PetscCall(EPSSetFromOptions(ctx->eps));
     722          31 :   PetscFunctionReturn(PETSC_SUCCESS);
     723             : }
     724             : 
     725           4 : static PetscErrorCode PEPLinearSetLinearization_Linear(PEP pep,PetscReal alpha,PetscReal beta)
     726             : {
     727           4 :   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
     728             : 
     729           4 :   PetscFunctionBegin;
     730           4 :   PetscCheck(beta!=0.0 || alpha!=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
     731           4 :   ctx->alpha = alpha;
     732           4 :   ctx->beta  = beta;
     733           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     734             : }
     735             : 
     736             : /*@
     737             :    PEPLinearSetLinearization - Set the coefficients that define
     738             :    the linearization of a quadratic eigenproblem.
     739             : 
     740             :    Logically Collective
     741             : 
     742             :    Input Parameters:
     743             : +  pep   - polynomial eigenvalue solver
     744             : .  alpha - first parameter of the linearization
     745             : -  beta  - second parameter of the linearization
     746             : 
     747             :    Options Database Key:
     748             : .  -pep_linear_linearization <alpha,beta> - Sets the coefficients
     749             : 
     750             :    Notes:
     751             :    Cannot pass zero for both alpha and beta. The default values are
     752             :    alpha=1 and beta=0.
     753             : 
     754             :    Level: advanced
     755             : 
     756             : .seealso: PEPLinearGetLinearization()
     757             : @*/
     758           4 : PetscErrorCode PEPLinearSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
     759             : {
     760           4 :   PetscFunctionBegin;
     761           4 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     762          12 :   PetscValidLogicalCollectiveReal(pep,alpha,2);
     763          12 :   PetscValidLogicalCollectiveReal(pep,beta,3);
     764           4 :   PetscTryMethod(pep,"PEPLinearSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
     765           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     766             : }
     767             : 
     768           1 : static PetscErrorCode PEPLinearGetLinearization_Linear(PEP pep,PetscReal *alpha,PetscReal *beta)
     769             : {
     770           1 :   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
     771             : 
     772           1 :   PetscFunctionBegin;
     773           1 :   if (alpha) *alpha = ctx->alpha;
     774           1 :   if (beta)  *beta  = ctx->beta;
     775           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     776             : }
     777             : 
     778             : /*@
     779             :    PEPLinearGetLinearization - Returns the coefficients that define
     780             :    the linearization of a quadratic eigenproblem.
     781             : 
     782             :    Not Collective
     783             : 
     784             :    Input Parameter:
     785             : .  pep  - polynomial eigenvalue solver
     786             : 
     787             :    Output Parameters:
     788             : +  alpha - the first parameter of the linearization
     789             : -  beta  - the second parameter of the linearization
     790             : 
     791             :    Level: advanced
     792             : 
     793             : .seealso: PEPLinearSetLinearization()
     794             : @*/
     795           1 : PetscErrorCode PEPLinearGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
     796             : {
     797           1 :   PetscFunctionBegin;
     798           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     799           1 :   PetscUseMethod(pep,"PEPLinearGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
     800           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     801             : }
     802             : 
     803          11 : static PetscErrorCode PEPLinearSetExplicitMatrix_Linear(PEP pep,PetscBool explicitmatrix)
     804             : {
     805          11 :   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
     806             : 
     807          11 :   PetscFunctionBegin;
     808          11 :   if (ctx->explicitmatrix != explicitmatrix) {
     809          11 :     ctx->explicitmatrix = explicitmatrix;
     810          11 :     pep->state = PEP_STATE_INITIAL;
     811             :   }
     812          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     813             : }
     814             : 
     815             : /*@
     816             :    PEPLinearSetExplicitMatrix - Indicate if the matrices A and B for the
     817             :    linearization of the problem must be built explicitly.
     818             : 
     819             :    Logically Collective
     820             : 
     821             :    Input Parameters:
     822             : +  pep         - polynomial eigenvalue solver
     823             : -  explicitmat - boolean flag indicating if the matrices are built explicitly
     824             : 
     825             :    Options Database Key:
     826             : .  -pep_linear_explicitmatrix <boolean> - Indicates the boolean flag
     827             : 
     828             :    Level: advanced
     829             : 
     830             : .seealso: PEPLinearGetExplicitMatrix()
     831             : @*/
     832          11 : PetscErrorCode PEPLinearSetExplicitMatrix(PEP pep,PetscBool explicitmat)
     833             : {
     834          11 :   PetscFunctionBegin;
     835          11 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     836          33 :   PetscValidLogicalCollectiveBool(pep,explicitmat,2);
     837          11 :   PetscTryMethod(pep,"PEPLinearSetExplicitMatrix_C",(PEP,PetscBool),(pep,explicitmat));
     838          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     839             : }
     840             : 
     841           1 : static PetscErrorCode PEPLinearGetExplicitMatrix_Linear(PEP pep,PetscBool *explicitmat)
     842             : {
     843           1 :   PEP_LINEAR *ctx = (PEP_LINEAR*)pep->data;
     844             : 
     845           1 :   PetscFunctionBegin;
     846           1 :   *explicitmat = ctx->explicitmatrix;
     847           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     848             : }
     849             : 
     850             : /*@
     851             :    PEPLinearGetExplicitMatrix - Returns the flag indicating if the matrices
     852             :    A and B for the linearization are built explicitly.
     853             : 
     854             :    Not Collective
     855             : 
     856             :    Input Parameter:
     857             : .  pep  - polynomial eigenvalue solver
     858             : 
     859             :    Output Parameter:
     860             : .  explicitmat - the mode flag
     861             : 
     862             :    Level: advanced
     863             : 
     864             : .seealso: PEPLinearSetExplicitMatrix()
     865             : @*/
     866           1 : PetscErrorCode PEPLinearGetExplicitMatrix(PEP pep,PetscBool *explicitmat)
     867             : {
     868           1 :   PetscFunctionBegin;
     869           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     870           1 :   PetscAssertPointer(explicitmat,2);
     871           1 :   PetscUseMethod(pep,"PEPLinearGetExplicitMatrix_C",(PEP,PetscBool*),(pep,explicitmat));
     872           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     873             : }
     874             : 
     875           1 : static PetscErrorCode PEPLinearSetEPS_Linear(PEP pep,EPS eps)
     876             : {
     877           1 :   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
     878             : 
     879           1 :   PetscFunctionBegin;
     880           1 :   PetscCall(PetscObjectReference((PetscObject)eps));
     881           1 :   PetscCall(EPSDestroy(&ctx->eps));
     882           1 :   ctx->eps     = eps;
     883           1 :   ctx->usereps = PETSC_TRUE;
     884           1 :   pep->state   = PEP_STATE_INITIAL;
     885           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     886             : }
     887             : 
     888             : /*@
     889             :    PEPLinearSetEPS - Associate an eigensolver object (EPS) to the
     890             :    polynomial eigenvalue solver.
     891             : 
     892             :    Collective
     893             : 
     894             :    Input Parameters:
     895             : +  pep - polynomial eigenvalue solver
     896             : -  eps - the eigensolver object
     897             : 
     898             :    Level: advanced
     899             : 
     900             : .seealso: PEPLinearGetEPS()
     901             : @*/
     902           1 : PetscErrorCode PEPLinearSetEPS(PEP pep,EPS eps)
     903             : {
     904           1 :   PetscFunctionBegin;
     905           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     906           1 :   PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
     907           1 :   PetscCheckSameComm(pep,1,eps,2);
     908           1 :   PetscTryMethod(pep,"PEPLinearSetEPS_C",(PEP,EPS),(pep,eps));
     909           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     910             : }
     911             : 
     912          32 : static PetscErrorCode PEPLinearGetEPS_Linear(PEP pep,EPS *eps)
     913             : {
     914          32 :   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
     915             : 
     916          32 :   PetscFunctionBegin;
     917          32 :   if (!ctx->eps) {
     918          32 :     PetscCall(EPSCreate(PetscObjectComm((PetscObject)pep),&ctx->eps));
     919          32 :     PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->eps,(PetscObject)pep,1));
     920          32 :     PetscCall(EPSSetOptionsPrefix(ctx->eps,((PetscObject)pep)->prefix));
     921          32 :     PetscCall(EPSAppendOptionsPrefix(ctx->eps,"pep_linear_"));
     922          32 :     PetscCall(PetscObjectSetOptions((PetscObject)ctx->eps,((PetscObject)pep)->options));
     923          32 :     PetscCall(EPSMonitorSet(ctx->eps,EPSMonitor_Linear,pep,NULL));
     924             :   }
     925          32 :   *eps = ctx->eps;
     926          32 :   PetscFunctionReturn(PETSC_SUCCESS);
     927             : }
     928             : 
     929             : /*@
     930             :    PEPLinearGetEPS - Retrieve the eigensolver object (EPS) associated
     931             :    to the polynomial eigenvalue solver.
     932             : 
     933             :    Collective
     934             : 
     935             :    Input Parameter:
     936             : .  pep - polynomial eigenvalue solver
     937             : 
     938             :    Output Parameter:
     939             : .  eps - the eigensolver object
     940             : 
     941             :    Level: advanced
     942             : 
     943             : .seealso: PEPLinearSetEPS()
     944             : @*/
     945          32 : PetscErrorCode PEPLinearGetEPS(PEP pep,EPS *eps)
     946             : {
     947          32 :   PetscFunctionBegin;
     948          32 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     949          32 :   PetscAssertPointer(eps,2);
     950          32 :   PetscUseMethod(pep,"PEPLinearGetEPS_C",(PEP,EPS*),(pep,eps));
     951          32 :   PetscFunctionReturn(PETSC_SUCCESS);
     952             : }
     953             : 
     954           0 : static PetscErrorCode PEPView_Linear(PEP pep,PetscViewer viewer)
     955             : {
     956           0 :   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
     957           0 :   PetscBool      isascii;
     958             : 
     959           0 :   PetscFunctionBegin;
     960           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     961           0 :   if (isascii) {
     962           0 :     if (!ctx->eps) PetscCall(PEPLinearGetEPS(pep,&ctx->eps));
     963           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  %s matrices\n",ctx->explicitmatrix? "explicit": "implicit"));
     964           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta));
     965           0 :     PetscCall(PetscViewerASCIIPushTab(viewer));
     966           0 :     PetscCall(EPSView(ctx->eps,viewer));
     967           0 :     PetscCall(PetscViewerASCIIPopTab(viewer));
     968             :   }
     969           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     970             : }
     971             : 
     972          34 : static PetscErrorCode PEPReset_Linear(PEP pep)
     973             : {
     974          34 :   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
     975             : 
     976          34 :   PetscFunctionBegin;
     977          34 :   if (!ctx->eps) PetscCall(EPSReset(ctx->eps));
     978          34 :   PetscCall(MatDestroy(&ctx->A));
     979          34 :   PetscCall(MatDestroy(&ctx->B));
     980          34 :   PetscCall(VecDestroy(&ctx->w[0]));
     981          34 :   PetscCall(VecDestroy(&ctx->w[1]));
     982          34 :   PetscCall(VecDestroy(&ctx->w[2]));
     983          34 :   PetscCall(VecDestroy(&ctx->w[3]));
     984          34 :   PetscCall(VecDestroy(&ctx->w[4]));
     985          34 :   PetscCall(VecDestroy(&ctx->w[5]));
     986          34 :   PetscFunctionReturn(PETSC_SUCCESS);
     987             : }
     988             : 
     989          33 : static PetscErrorCode PEPDestroy_Linear(PEP pep)
     990             : {
     991          33 :   PEP_LINEAR     *ctx = (PEP_LINEAR*)pep->data;
     992             : 
     993          33 :   PetscFunctionBegin;
     994          33 :   PetscCall(EPSDestroy(&ctx->eps));
     995          33 :   PetscCall(PetscFree(pep->data));
     996          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",NULL));
     997          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",NULL));
     998          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",NULL));
     999          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",NULL));
    1000          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",NULL));
    1001          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",NULL));
    1002          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1003             : }
    1004             : 
    1005          33 : SLEPC_EXTERN PetscErrorCode PEPCreate_Linear(PEP pep)
    1006             : {
    1007          33 :   PEP_LINEAR     *ctx;
    1008             : 
    1009          33 :   PetscFunctionBegin;
    1010          33 :   PetscCall(PetscNew(&ctx));
    1011          33 :   pep->data = (void*)ctx;
    1012             : 
    1013          33 :   pep->lineariz       = PETSC_TRUE;
    1014          33 :   ctx->explicitmatrix = PETSC_FALSE;
    1015          33 :   ctx->alpha          = 1.0;
    1016          33 :   ctx->beta           = 0.0;
    1017             : 
    1018          33 :   pep->ops->solve          = PEPSolve_Linear;
    1019          33 :   pep->ops->setup          = PEPSetUp_Linear;
    1020          33 :   pep->ops->setfromoptions = PEPSetFromOptions_Linear;
    1021          33 :   pep->ops->destroy        = PEPDestroy_Linear;
    1022          33 :   pep->ops->reset          = PEPReset_Linear;
    1023          33 :   pep->ops->view           = PEPView_Linear;
    1024          33 :   pep->ops->backtransform  = PEPBackTransform_Default;
    1025          33 :   pep->ops->computevectors = PEPComputeVectors_Default;
    1026          33 :   pep->ops->extractvectors = PEPExtractVectors_Linear;
    1027             : 
    1028          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetLinearization_C",PEPLinearSetLinearization_Linear));
    1029          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetLinearization_C",PEPLinearGetLinearization_Linear));
    1030          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetEPS_C",PEPLinearSetEPS_Linear));
    1031          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetEPS_C",PEPLinearGetEPS_Linear));
    1032          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearSetExplicitMatrix_C",PEPLinearSetExplicitMatrix_Linear));
    1033          33 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPLinearGetExplicitMatrix_C",PEPLinearGetExplicitMatrix_Linear));
    1034          33 :   PetscFunctionReturn(PETSC_SUCCESS);
    1035             : }