LCOV - code coverage report
Current view: top level - pep/impls/krylov/stoar - stoar.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 517 549 94.2 %
Date: 2024-03-29 00:34:52 Functions: 32 33 97.0 %
Legend: Lines: hit not hit

          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             :    SLEPc polynomial eigensolver: "stoar"
      12             : 
      13             :    Method: S-TOAR
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Symmetric Two-Level Orthogonal Arnoldi.
      18             : 
      19             :    References:
      20             : 
      21             :        [1] C. Campos and J.E. Roman, "Restarted Q-Arnoldi-type methods
      22             :            exploiting symmetry in quadratic eigenvalue problems", BIT
      23             :            Numer. Math. 56(4):1213-1236, 2016.
      24             : */
      25             : 
      26             : #include <slepc/private/pepimpl.h>         /*I "slepcpep.h" I*/
      27             : #include "../src/pep/impls/krylov/pepkrylov.h"
      28             : #include <slepcblaslapack.h>
      29             : 
      30             : static PetscBool  cited = PETSC_FALSE;
      31             : static const char citation[] =
      32             :   "@Article{slepc-stoar,\n"
      33             :   "   author = \"C. Campos and J. E. Roman\",\n"
      34             :   "   title = \"Restarted {Q-Arnoldi-type} methods exploiting symmetry in quadratic eigenvalue problems\",\n"
      35             :   "   journal = \"{BIT} Numer. Math.\",\n"
      36             :   "   volume = \"56\",\n"
      37             :   "   number = \"4\",\n"
      38             :   "   pages = \"1213--1236\",\n"
      39             :   "   year = \"2016,\"\n"
      40             :   "   doi = \"https://doi.org/10.1007/s10543-016-0601-5\"\n"
      41             :   "}\n";
      42             : 
      43             : typedef struct {
      44             :   PetscReal   scal[2];
      45             :   Mat         A[2];
      46             :   Vec         t;
      47             : } PEP_STOAR_MATSHELL;
      48             : 
      49        9024 : static PetscErrorCode MatMult_STOAR(Mat A,Vec x,Vec y)
      50             : {
      51        9024 :   PEP_STOAR_MATSHELL *ctx;
      52             : 
      53        9024 :   PetscFunctionBegin;
      54        9024 :   PetscCall(MatShellGetContext(A,&ctx));
      55        9024 :   PetscCall(MatMult(ctx->A[0],x,y));
      56        9024 :   PetscCall(VecScale(y,ctx->scal[0]));
      57        9024 :   if (ctx->scal[1]) {
      58          28 :     PetscCall(MatMult(ctx->A[1],x,ctx->t));
      59          28 :     PetscCall(VecAXPY(y,ctx->scal[1],ctx->t));
      60             :   }
      61        9024 :   PetscFunctionReturn(PETSC_SUCCESS);
      62             : }
      63             : 
      64          60 : static PetscErrorCode MatDestroy_STOAR(Mat A)
      65             : {
      66          60 :   PEP_STOAR_MATSHELL *ctx;
      67             : 
      68          60 :   PetscFunctionBegin;
      69          60 :   PetscCall(MatShellGetContext(A,&ctx));
      70          60 :   PetscCall(VecDestroy(&ctx->t));
      71          60 :   PetscCall(PetscFree(ctx));
      72          60 :   PetscFunctionReturn(PETSC_SUCCESS);
      73             : }
      74             : 
      75          20 : PetscErrorCode PEPSTOARSetUpInnerMatrix(PEP pep,Mat *B)
      76             : {
      77          20 :   Mat                pB[4],Bs[3],D[3];
      78          20 :   PetscInt           i,j,n,m;
      79          20 :   PEP_STOAR_MATSHELL *ctxMat[3];
      80          20 :   PEP_STOAR          *ctx=(PEP_STOAR*)pep->data;
      81             : 
      82          20 :   PetscFunctionBegin;
      83          80 :   for (i=0;i<3;i++) {
      84          60 :     PetscCall(STGetMatrixTransformed(pep->st,i,&D[i])); /* D[2] = M */
      85             :   }
      86          20 :   PetscCall(MatGetLocalSize(D[2],&m,&n));
      87             : 
      88          80 :   for (j=0;j<3;j++) {
      89          60 :     PetscCall(PetscNew(ctxMat+j));
      90          60 :     PetscCall(MatCreateShell(PetscObjectComm((PetscObject)pep),m,n,PETSC_DETERMINE,PETSC_DETERMINE,ctxMat[j],&Bs[j]));
      91          60 :     PetscCall(MatShellSetOperation(Bs[j],MATOP_MULT,(void(*)(void))MatMult_STOAR));
      92          60 :     PetscCall(MatShellSetOperation(Bs[j],MATOP_DESTROY,(void(*)(void))MatDestroy_STOAR));
      93             :   }
      94         100 :   for (i=0;i<4;i++) pB[i] = NULL;
      95          20 :   if (ctx->alpha) {
      96          19 :     ctxMat[0]->A[0] = D[0]; ctxMat[0]->scal[0] = ctx->alpha; ctxMat[0]->scal[1] = 0.0;
      97          19 :     ctxMat[2]->A[0] = D[2]; ctxMat[2]->scal[0] = -ctx->alpha*pep->sfactor*pep->sfactor; ctxMat[2]->scal[1] = 0.0;
      98          19 :     pB[0] = Bs[0]; pB[3] = Bs[2];
      99             :   }
     100          20 :   if (ctx->beta) {
     101           2 :     i = (ctx->alpha)?1:0;
     102           2 :     ctxMat[0]->scal[1] = 0.0;
     103           2 :     ctxMat[0]->A[i] = D[1]; ctxMat[0]->scal[i] = -ctx->beta*pep->sfactor;
     104           2 :     ctxMat[1]->A[0] = D[2]; ctxMat[1]->scal[0] = -ctx->beta*pep->sfactor*pep->sfactor; ctxMat[1]->scal[1] = 0.0;
     105           2 :     pB[0] = Bs[0]; pB[1] = pB[2] = Bs[1];
     106             :   }
     107          20 :   PetscCall(BVCreateVec(pep->V,&ctxMat[0]->t));
     108          20 :   PetscCall(MatCreateNest(PetscObjectComm((PetscObject)pep),2,NULL,2,NULL,pB,B));
     109          80 :   for (j=0;j<3;j++) PetscCall(MatDestroy(&Bs[j]));
     110          20 :   PetscFunctionReturn(PETSC_SUCCESS);
     111             : }
     112             : 
     113          21 : static PetscErrorCode PEPSetUp_STOAR(PEP pep)
     114             : {
     115          21 :   PetscBool         sinv,flg;
     116          21 :   PEP_STOAR         *ctx = (PEP_STOAR*)pep->data;
     117          21 :   PetscInt          ld,i;
     118          21 :   PetscReal         eta;
     119          21 :   BVOrthogType      otype;
     120          21 :   BVOrthogBlockType obtype;
     121             : 
     122          21 :   PetscFunctionBegin;
     123          21 :   PEPCheckHermitian(pep);
     124          21 :   PEPCheckQuadratic(pep);
     125          21 :   PEPCheckShiftSinvert(pep);
     126             :   /* spectrum slicing requires special treatment of default values */
     127          21 :   if (pep->which==PEP_ALL) {
     128           8 :     pep->ops->solve = PEPSolve_STOAR_QSlice;
     129           8 :     pep->ops->extractvectors = NULL;
     130           8 :     pep->ops->setdefaultst   = NULL;
     131           8 :     PetscCall(PEPSetUp_STOAR_QSlice(pep));
     132             :   } else {
     133          13 :     PetscCall(PEPSetDimensions_Default(pep,pep->nev,&pep->ncv,&pep->mpd));
     134          13 :     PetscCheck(ctx->lock || pep->mpd>=pep->ncv,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
     135          13 :     if (pep->max_it==PETSC_DEFAULT) pep->max_it = PetscMax(100,2*(pep->nmat-1)*pep->n/pep->ncv);
     136          13 :     pep->ops->solve = PEPSolve_STOAR;
     137          13 :     ld   = pep->ncv+2;
     138          13 :     PetscCall(DSSetType(pep->ds,DSGHIEP));
     139          13 :     PetscCall(DSSetCompact(pep->ds,PETSC_TRUE));
     140          13 :     PetscCall(DSSetExtraRow(pep->ds,PETSC_TRUE));
     141          13 :     PetscCall(DSAllocate(pep->ds,ld));
     142          13 :     PetscCall(PEPBasisCoefficients(pep,pep->pbc));
     143          13 :     PetscCall(STGetTransform(pep->st,&flg));
     144          13 :     if (!flg) {
     145          11 :       PetscCall(PetscFree(pep->solvematcoeffs));
     146          11 :       PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
     147          11 :       PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
     148          11 :       if (sinv) PetscCall(PEPEvaluateBasis(pep,pep->target,0,pep->solvematcoeffs,NULL));
     149             :       else {
     150          15 :         for (i=0;i<pep->nmat-1;i++) pep->solvematcoeffs[i] = 0.0;
     151           5 :         pep->solvematcoeffs[pep->nmat-1] = 1.0;
     152             :       }
     153             :     }
     154             :   }
     155          21 :   if (!pep->which) PetscCall(PEPSetWhichEigenpairs_Default(pep));
     156          21 :   PEPCheckUnsupported(pep,PEP_FEATURE_NONMONOMIAL | PEP_FEATURE_REGION);
     157             : 
     158          21 :   PetscCall(PEPAllocateSolution(pep,2));
     159          21 :   PetscCall(PEPSetWorkVecs(pep,4));
     160          21 :   PetscCall(BVDestroy(&ctx->V));
     161          21 :   PetscCall(BVCreateTensor(pep->V,pep->nmat-1,&ctx->V));
     162          21 :   PetscCall(BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype));
     163          21 :   PetscCall(BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
     164          21 :   PetscFunctionReturn(PETSC_SUCCESS);
     165             : }
     166             : 
     167             : /*
     168             :   Compute a run of Lanczos iterations. dim(work)=(ctx->ld)*4
     169             : */
     170         112 : static PetscErrorCode PEPSTOARrun(PEP pep,PetscReal *a,PetscReal *b,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,Vec *t_)
     171             : {
     172         112 :   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
     173         112 :   PetscInt       i,j,m=*M,l,lock;
     174         112 :   PetscInt       lds,d,ld,offq,nqt,ldds;
     175         112 :   Vec            v=t_[0],t=t_[1],q=t_[2];
     176         112 :   PetscReal      norm,sym=0.0,fro=0.0,*f;
     177         112 :   PetscScalar    *y,*S,*x,sigma;
     178         112 :   PetscBLASInt   j_,one=1;
     179         112 :   PetscBool      lindep,flg,sinvert=PETSC_FALSE;
     180         112 :   Mat            MS;
     181             : 
     182         112 :   PetscFunctionBegin;
     183         112 :   PetscCall(PetscMalloc1(*M,&y));
     184         112 :   PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
     185         112 :   PetscCall(BVTensorGetDegree(ctx->V,&d));
     186         112 :   PetscCall(BVGetActiveColumns(pep->V,&lock,&nqt));
     187         112 :   lds = d*ld;
     188         112 :   offq = ld;
     189         112 :   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
     190         112 :   *breakdown = PETSC_FALSE; /* ----- */
     191         112 :   PetscCall(DSGetDimensions(pep->ds,NULL,&l,NULL,NULL));
     192         112 :   PetscCall(BVSetActiveColumns(ctx->V,0,m));
     193         112 :   PetscCall(BVSetActiveColumns(pep->V,0,nqt));
     194         112 :   PetscCall(STGetTransform(pep->st,&flg));
     195         112 :   if (!flg) {
     196             :     /* spectral transformation handled by the solver */
     197         106 :     PetscCall(PetscObjectTypeCompareAny((PetscObject)pep->st,&flg,STSINVERT,STSHIFT,""));
     198         106 :     PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"ST type not supported for TOAR without transforming matrices");
     199         106 :     PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinvert));
     200         106 :     PetscCall(STGetShift(pep->st,&sigma));
     201             :   }
     202        1236 :   for (j=k;j<m;j++) {
     203             :     /* apply operator */
     204        1124 :     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
     205        1124 :     PetscCall(MatDenseGetArray(MS,&S));
     206        1124 :     PetscCall(BVGetColumn(pep->V,nqt,&t));
     207        1124 :     PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+j*lds));
     208        1124 :     if (!sinvert) {
     209         848 :       PetscCall(STMatMult(pep->st,0,v,q));
     210         848 :       PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds));
     211         848 :       PetscCall(STMatMult(pep->st,1,v,t));
     212         848 :       PetscCall(VecAXPY(q,pep->sfactor,t));
     213         848 :       if (ctx->beta && ctx->alpha) {
     214           0 :         PetscCall(STMatMult(pep->st,2,v,t));
     215           0 :         PetscCall(VecAXPY(q,-pep->sfactor*pep->sfactor*ctx->beta/ctx->alpha,t));
     216             :       }
     217         848 :       PetscCall(STMatSolve(pep->st,q,t));
     218         848 :       PetscCall(VecScale(t,-1.0/(pep->sfactor*pep->sfactor)));
     219             :     } else {
     220         276 :       PetscCall(STMatMult(pep->st,1,v,q));
     221         276 :       PetscCall(STMatMult(pep->st,2,v,t));
     222         276 :       PetscCall(VecAXPY(q,sigma*pep->sfactor,t));
     223         276 :       PetscCall(VecScale(q,pep->sfactor));
     224         276 :       PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds));
     225         276 :       PetscCall(STMatMult(pep->st,2,v,t));
     226         276 :       PetscCall(VecAXPY(q,pep->sfactor*pep->sfactor,t));
     227         276 :       PetscCall(STMatSolve(pep->st,q,t));
     228         276 :       PetscCall(VecScale(t,-1.0));
     229             :     }
     230        1124 :     PetscCall(BVRestoreColumn(pep->V,nqt,&t));
     231             : 
     232             :     /* orthogonalize */
     233        1124 :     if (!sinvert) x = S+offq+(j+1)*lds;
     234         276 :     else x = S+(j+1)*lds;
     235        1124 :     PetscCall(BVOrthogonalizeColumn(pep->V,nqt,x,&norm,&lindep));
     236             : 
     237        1124 :     if (!lindep) {
     238        1124 :       if (!sinvert) *(S+offq+(j+1)*lds+nqt) = norm;
     239         276 :       else *(S+(j+1)*lds+nqt) = norm;
     240        1124 :       PetscCall(BVScaleColumn(pep->V,nqt,1.0/norm));
     241        1124 :       nqt++;
     242             :     }
     243        1124 :     if (!sinvert) {
     244       13963 :       for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+i) = *(S+offq+j*lds+i);
     245         848 :       if (ctx->beta && ctx->alpha) {
     246           0 :         for (i=0;i<=nqt-1;i++) *(S+(j+1)*lds+offq+i) -= *(S+(j+1)*lds+i)*ctx->beta/ctx->alpha;
     247             :       }
     248        4367 :     } else for (i=0;i<nqt;i++) *(S+(j+1)*lds+offq+i) = *(S+j*lds+i)+sigma*(*(S+(j+1)*lds+i));
     249        1124 :     PetscCall(BVSetActiveColumns(pep->V,0,nqt));
     250        1124 :     PetscCall(MatDenseRestoreArray(MS,&S));
     251        1124 :     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
     252             : 
     253             :     /* level-2 orthogonalization */
     254        1124 :     PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep));
     255        1124 :     a[j] = PetscRealPart(y[j]);
     256        1124 :     omega[j+1] = (norm > 0)?1.0:-1.0;
     257        1124 :     PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
     258        1124 :     b[j] = PetscAbsReal(norm);
     259             : 
     260             :     /* check symmetry */
     261        1124 :     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&f));
     262        1124 :     if (j==k) {
     263         904 :       for (i=l;i<j-1;i++) y[i] = PetscAbsScalar(y[i])-PetscAbsReal(f[2*ldds+i]);
     264         127 :       for (i=0;i<l;i++) y[i] = 0.0;
     265             :     }
     266        1124 :     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&f));
     267        1124 :     if (j>0) y[j-1] = PetscAbsScalar(y[j-1])-PetscAbsReal(b[j-1]);
     268        1124 :     PetscCall(PetscBLASIntCast(j,&j_));
     269        1124 :     sym = SlepcAbs(BLASnrm2_(&j_,y,&one),sym);
     270        1124 :     fro = SlepcAbs(fro,SlepcAbs(a[j],b[j]));
     271        1124 :     if (j>0) fro = SlepcAbs(fro,b[j-1]);
     272        1963 :     if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*pep->tol)) {
     273           0 :       *symmlost = PETSC_TRUE;
     274           0 :       *M=j;
     275           0 :       break;
     276             :     }
     277             :   }
     278         112 :   PetscCall(BVSetActiveColumns(pep->V,lock,nqt));
     279         112 :   PetscCall(BVSetActiveColumns(ctx->V,0,*M));
     280         112 :   PetscCall(PetscFree(y));
     281         112 :   PetscFunctionReturn(PETSC_SUCCESS);
     282             : }
     283             : 
     284             : #if 0
     285             : static PetscErrorCode PEPSTOARpreKConvergence(PEP pep,PetscInt nv,PetscReal *norm,Vec *w)
     286             : {
     287             :   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
     288             :   PetscBLASInt   n_,one=1;
     289             :   PetscInt       lds=2*ctx->ld;
     290             :   PetscReal      t1,t2;
     291             :   PetscScalar    *S=ctx->S;
     292             : 
     293             :   PetscFunctionBegin;
     294             :   PetscCall(PetscBLASIntCast(nv+2,&n_));
     295             :   t1 = BLASnrm2_(&n_,S+nv*2*ctx->ld,&one);
     296             :   t2 = BLASnrm2_(&n_,S+(nv*2+1)*ctx->ld,&one);
     297             :   *norm = SlepcAbs(t1,t2);
     298             :   PetscCall(BVSetActiveColumns(pep->V,0,nv+2));
     299             :   PetscCall(BVMultVec(pep->V,1.0,0.0,w[1],S+nv*lds));
     300             :   PetscCall(STMatMult(pep->st,0,w[1],w[2]));
     301             :   PetscCall(VecNorm(w[2],NORM_2,&t1));
     302             :   PetscCall(BVMultVec(pep->V,1.0,0.0,w[1],S+ctx->ld+nv*lds));
     303             :   PetscCall(STMatMult(pep->st,2,w[1],w[2]));
     304             :   PetscCall(VecNorm(w[2],NORM_2,&t2));
     305             :   t2 *= pep->sfactor*pep->sfactor;
     306             :   *norm = PetscMax(*norm,SlepcAbs(t1,t2));
     307             :   PetscFunctionReturn(PETSC_SUCCESS);
     308             : }
     309             : #endif
     310             : 
     311          13 : PetscErrorCode PEPSolve_STOAR(PEP pep)
     312             : {
     313          13 :   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
     314          13 :   PetscInt       j,k,l,nv=0,ld,ldds,t,nq=0;
     315          13 :   PetscInt       nconv=0,deg=pep->nmat-1;
     316          13 :   PetscScalar    sigma;
     317          13 :   PetscReal      beta,norm=1.0,*omega,*a,*b;
     318          13 :   PetscBool      breakdown,symmlost=PETSC_FALSE,sinv=PETSC_FALSE,falselock=PETSC_TRUE,flg;
     319          13 :   Mat            MQ,A,D;
     320          13 :   Vec            vomega;
     321             : 
     322          13 :   PetscFunctionBegin;
     323          13 :   PetscCall(PetscCitationsRegister(citation,&cited));
     324          13 :   PetscCall(PEPSTOARSetUpInnerMatrix(pep,&A));
     325          13 :   PetscCall(BVSetMatrix(ctx->V,A,PETSC_TRUE));
     326          13 :   PetscCall(MatDestroy(&A));
     327          13 :   if (ctx->lock) {
     328             :     /* undocumented option to use a cheaper locking instead of the true locking */
     329          12 :     PetscCall(PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL));
     330             :   }
     331          13 :   PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
     332          13 :   PetscCall(STGetShift(pep->st,&sigma));
     333          13 :   PetscCall(STGetTransform(pep->st,&flg));
     334          13 :   if (pep->sfactor!=1.0) {
     335           5 :     if (!flg) {
     336           4 :       pep->target /= pep->sfactor;
     337           4 :       PetscCall(RGPushScale(pep->rg,1.0/pep->sfactor));
     338           4 :       PetscCall(STScaleShift(pep->st,1.0/pep->sfactor));
     339           4 :       sigma /= pep->sfactor;
     340             :     } else {
     341           1 :       PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
     342           1 :       pep->target = sinv?pep->target*pep->sfactor:pep->target/pep->sfactor;
     343           1 :       PetscCall(RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor));
     344           1 :       PetscCall(STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor));
     345             :     }
     346             :   }
     347          13 :   if (flg) sigma = 0.0;
     348             : 
     349             :   /* Get the starting Arnoldi vector */
     350          13 :   PetscCall(BVTensorBuildFirstColumn(ctx->V,pep->nini));
     351          13 :   PetscCall(DSSetDimensions(pep->ds,1,PETSC_DEFAULT,PETSC_DEFAULT));
     352          13 :   PetscCall(BVSetActiveColumns(ctx->V,0,1));
     353          13 :   PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
     354          13 :   PetscCall(BVGetSignature(ctx->V,vomega));
     355          13 :   PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
     356             : 
     357             :   /* Restart loop */
     358          13 :   l = 0;
     359          13 :   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
     360         125 :   while (pep->reason == PEP_CONVERGED_ITERATING) {
     361         112 :     pep->its++;
     362         112 :     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&a));
     363         112 :     b = a+ldds;
     364         112 :     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_D,&omega));
     365             : 
     366             :     /* Compute an nv-step Lanczos factorization */
     367         112 :     nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
     368         112 :     PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
     369         112 :     PetscCall(PEPSTOARrun(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work));
     370         112 :     beta = b[nv-1];
     371         112 :     if (symmlost && nv==pep->nconv+l) {
     372           0 :       pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
     373           0 :       pep->nconv = nconv;
     374           0 :       if (falselock || !ctx->lock) {
     375           0 :        PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
     376           0 :        PetscCall(BVTensorCompress(ctx->V,0));
     377             :       }
     378             :       break;
     379             :     }
     380         112 :     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&a));
     381         112 :     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega));
     382         112 :     PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
     383         112 :     if (l==0) PetscCall(DSSetState(pep->ds,DS_STATE_INTERMEDIATE));
     384          99 :     else PetscCall(DSSetState(pep->ds,DS_STATE_RAW));
     385             : 
     386             :     /* Solve projected problem */
     387         112 :     PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
     388         112 :     PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
     389         112 :     PetscCall(DSUpdateExtraRow(pep->ds));
     390         112 :     PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
     391             : 
     392             :     /* Check convergence */
     393             :     /* PetscCall(PEPSTOARpreKConvergence(pep,nv,&norm,pep->work));*/
     394         112 :     norm = 1.0;
     395         112 :     PetscCall(DSGetDimensions(pep->ds,NULL,NULL,NULL,&t));
     396         112 :     PetscCall(PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k));
     397         112 :     PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx));
     398             : 
     399             :     /* Update l */
     400         112 :     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
     401             :     else {
     402          99 :       l = PetscMax(1,(PetscInt)((nv-k)/2));
     403          99 :       l = PetscMin(l,t);
     404          99 :       PetscCall(DSGetTruncateSize(pep->ds,k,t,&l));
     405          99 :       if (!breakdown) {
     406             :         /* Prepare the Rayleigh quotient for restart */
     407          99 :         PetscCall(DSTruncate(pep->ds,k+l,PETSC_FALSE));
     408             :       }
     409             :     }
     410         112 :     nconv = k;
     411         112 :     if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
     412         112 :     if (l) PetscCall(PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
     413             : 
     414             :     /* Update S */
     415         112 :     PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&MQ));
     416         112 :     PetscCall(BVMultInPlace(ctx->V,MQ,pep->nconv,k+l));
     417         112 :     PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&MQ));
     418             : 
     419             :     /* Copy last column of S */
     420         112 :     PetscCall(BVCopyColumn(ctx->V,nv,k+l));
     421         112 :     PetscCall(BVSetActiveColumns(ctx->V,0,k+l));
     422         112 :     PetscCall(DSSetDimensions(pep->ds,k+l,PETSC_DEFAULT,PETSC_DEFAULT));
     423         112 :     PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
     424         112 :     PetscCall(BVSetSignature(ctx->V,vomega));
     425         112 :     PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
     426             : 
     427         112 :     if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
     428             :       /* stop if breakdown */
     429           0 :       PetscCall(PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta));
     430           0 :       pep->reason = PEP_DIVERGED_BREAKDOWN;
     431             :     }
     432         112 :     if (pep->reason != PEP_CONVERGED_ITERATING) l--;
     433         112 :     PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
     434         112 :     if (k+l+deg<=nq) {
     435         112 :       PetscCall(BVSetActiveColumns(ctx->V,pep->nconv,k+l+1));
     436         112 :       if (!falselock && ctx->lock) PetscCall(BVTensorCompress(ctx->V,k-pep->nconv));
     437         112 :       else PetscCall(BVTensorCompress(ctx->V,0));
     438             :     }
     439         112 :     pep->nconv = k;
     440         125 :     PetscCall(PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv));
     441             :   }
     442             : 
     443          13 :   if (pep->nconv>0) {
     444          13 :     PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
     445          13 :     PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
     446          13 :     PetscCall(BVSetActiveColumns(pep->V,0,nq));
     447          13 :     if (nq>pep->nconv) {
     448          13 :       PetscCall(BVTensorCompress(ctx->V,pep->nconv));
     449          13 :       PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
     450             :     }
     451             :   }
     452          13 :   PetscCall(STGetTransform(pep->st,&flg));
     453          13 :   if (!flg) PetscTryTypeMethod(pep,backtransform);
     454          13 :   if (pep->sfactor!=1.0) {
     455          19 :     for (j=0;j<pep->nconv;j++) {
     456          14 :       pep->eigr[j] *= pep->sfactor;
     457          14 :       pep->eigi[j] *= pep->sfactor;
     458             :     }
     459             :   }
     460             :   /* restore original values */
     461          13 :   if (!flg) {
     462          11 :     pep->target *= pep->sfactor;
     463          11 :     PetscCall(STScaleShift(pep->st,pep->sfactor));
     464             :   } else {
     465           2 :     PetscCall(STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor));
     466           2 :     pep->target = (sinv)?pep->target/pep->sfactor:pep->target*pep->sfactor;
     467             :   }
     468          13 :   if (pep->sfactor!=1.0) PetscCall(RGPopScale(pep->rg));
     469             : 
     470          13 :   PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
     471          13 :   PetscFunctionReturn(PETSC_SUCCESS);
     472             : }
     473             : 
     474          20 : static PetscErrorCode PEPSetFromOptions_STOAR(PEP pep,PetscOptionItems *PetscOptionsObject)
     475             : {
     476          20 :   PetscBool      flg,lock,b,f1,f2,f3;
     477          20 :   PetscInt       i,j,k;
     478          20 :   PetscReal      array[2]={0,0};
     479          20 :   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
     480             : 
     481          20 :   PetscFunctionBegin;
     482          20 :   PetscOptionsHeadBegin(PetscOptionsObject,"PEP STOAR Options");
     483             : 
     484          20 :     PetscCall(PetscOptionsBool("-pep_stoar_locking","Choose between locking and non-locking variants","PEPSTOARSetLocking",PETSC_FALSE,&lock,&flg));
     485          20 :     if (flg) PetscCall(PEPSTOARSetLocking(pep,lock));
     486             : 
     487          20 :     b = ctx->detect;
     488          20 :     PetscCall(PetscOptionsBool("-pep_stoar_detect_zeros","Check zeros during factorizations at interval boundaries","PEPSTOARSetDetectZeros",ctx->detect,&b,&flg));
     489          20 :     if (flg) PetscCall(PEPSTOARSetDetectZeros(pep,b));
     490             : 
     491          20 :     i = 1;
     492          20 :     j = k = PETSC_DECIDE;
     493          20 :     PetscCall(PetscOptionsInt("-pep_stoar_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",20,&i,&f1));
     494          20 :     PetscCall(PetscOptionsInt("-pep_stoar_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&j,&f2));
     495          20 :     PetscCall(PetscOptionsInt("-pep_stoar_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","PEPSTOARSetDimensions",40,&k,&f3));
     496          20 :     if (f1 || f2 || f3) PetscCall(PEPSTOARSetDimensions(pep,i,j,k));
     497             : 
     498          20 :     k = 2;
     499          20 :     PetscCall(PetscOptionsRealArray("-pep_stoar_linearization","Parameters of the linearization","PEPSTOARSetLinearization",array,&k,&flg));
     500          20 :     if (flg) PetscCall(PEPSTOARSetLinearization(pep,array[0],array[1]));
     501             : 
     502          20 :     b = ctx->checket;
     503          20 :     PetscCall(PetscOptionsBool("-pep_stoar_check_eigenvalue_type","Check eigenvalue type during spectrum slicing","PEPSTOARSetCheckEigenvalueType",ctx->checket,&b,&flg));
     504          20 :     if (flg) PetscCall(PEPSTOARSetCheckEigenvalueType(pep,b));
     505             : 
     506          20 :   PetscOptionsHeadEnd();
     507          20 :   PetscFunctionReturn(PETSC_SUCCESS);
     508             : }
     509             : 
     510           2 : static PetscErrorCode PEPSTOARSetLocking_STOAR(PEP pep,PetscBool lock)
     511             : {
     512           2 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     513             : 
     514           2 :   PetscFunctionBegin;
     515           2 :   ctx->lock = lock;
     516           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     517             : }
     518             : 
     519             : /*@
     520             :    PEPSTOARSetLocking - Choose between locking and non-locking variants of
     521             :    the STOAR method.
     522             : 
     523             :    Logically Collective
     524             : 
     525             :    Input Parameters:
     526             : +  pep  - the eigenproblem solver context
     527             : -  lock - true if the locking variant must be selected
     528             : 
     529             :    Options Database Key:
     530             : .  -pep_stoar_locking - Sets the locking flag
     531             : 
     532             :    Notes:
     533             :    The default is to lock converged eigenpairs when the method restarts.
     534             :    This behaviour can be changed so that all directions are kept in the
     535             :    working subspace even if already converged to working accuracy (the
     536             :    non-locking variant).
     537             : 
     538             :    Level: advanced
     539             : 
     540             : .seealso: PEPSTOARGetLocking()
     541             : @*/
     542           2 : PetscErrorCode PEPSTOARSetLocking(PEP pep,PetscBool lock)
     543             : {
     544           2 :   PetscFunctionBegin;
     545           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     546           8 :   PetscValidLogicalCollectiveBool(pep,lock,2);
     547           2 :   PetscTryMethod(pep,"PEPSTOARSetLocking_C",(PEP,PetscBool),(pep,lock));
     548           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     549             : }
     550             : 
     551           2 : static PetscErrorCode PEPSTOARGetLocking_STOAR(PEP pep,PetscBool *lock)
     552             : {
     553           2 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     554             : 
     555           2 :   PetscFunctionBegin;
     556           2 :   *lock = ctx->lock;
     557           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     558             : }
     559             : 
     560             : /*@
     561             :    PEPSTOARGetLocking - Gets the locking flag used in the STOAR method.
     562             : 
     563             :    Not Collective
     564             : 
     565             :    Input Parameter:
     566             : .  pep - the eigenproblem solver context
     567             : 
     568             :    Output Parameter:
     569             : .  lock - the locking flag
     570             : 
     571             :    Level: advanced
     572             : 
     573             : .seealso: PEPSTOARSetLocking()
     574             : @*/
     575           2 : PetscErrorCode PEPSTOARGetLocking(PEP pep,PetscBool *lock)
     576             : {
     577           2 :   PetscFunctionBegin;
     578           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     579           2 :   PetscAssertPointer(lock,2);
     580           2 :   PetscUseMethod(pep,"PEPSTOARGetLocking_C",(PEP,PetscBool*),(pep,lock));
     581           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     582             : }
     583             : 
     584           2 : static PetscErrorCode PEPSTOARGetInertias_STOAR(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
     585             : {
     586           2 :   PetscInt       i,numsh;
     587           2 :   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
     588           2 :   PEP_SR         sr = ctx->sr;
     589             : 
     590           2 :   PetscFunctionBegin;
     591           2 :   PetscCheck(pep->state,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Must call PEPSetUp() first");
     592           2 :   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see PEPSetInterval()");
     593           2 :   switch (pep->state) {
     594             :   case PEP_STATE_INITIAL:
     595             :     break;
     596           2 :   case PEP_STATE_SETUP:
     597           2 :     if (n) *n = 2;
     598           2 :     if (shifts) {
     599           2 :       PetscCall(PetscMalloc1(2,shifts));
     600           2 :       (*shifts)[0] = pep->inta;
     601           2 :       (*shifts)[1] = pep->intb;
     602             :     }
     603           2 :     if (inertias) {
     604           2 :       PetscCall(PetscMalloc1(2,inertias));
     605           2 :       (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
     606           2 :       (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
     607             :     }
     608             :     break;
     609           0 :   case PEP_STATE_SOLVED:
     610             :   case PEP_STATE_EIGENVECTORS:
     611           0 :     numsh = ctx->nshifts;
     612           0 :     if (n) *n = numsh;
     613           0 :     if (shifts) {
     614           0 :       PetscCall(PetscMalloc1(numsh,shifts));
     615           0 :       for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
     616             :     }
     617           0 :     if (inertias) {
     618           0 :       PetscCall(PetscMalloc1(numsh,inertias));
     619           0 :       for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
     620             :     }
     621             :     break;
     622             :   }
     623           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     624             : }
     625             : 
     626             : /*@C
     627             :    PEPSTOARGetInertias - Gets the values of the shifts and their
     628             :    corresponding inertias in case of doing spectrum slicing for a
     629             :    computational interval.
     630             : 
     631             :    Not Collective
     632             : 
     633             :    Input Parameter:
     634             : .  pep - the eigenproblem solver context
     635             : 
     636             :    Output Parameters:
     637             : +  n        - number of shifts, including the endpoints of the interval
     638             : .  shifts   - the values of the shifts used internally in the solver
     639             : -  inertias - the values of the inertia in each shift
     640             : 
     641             :    Notes:
     642             :    If called after PEPSolve(), all shifts used internally by the solver are
     643             :    returned (including both endpoints and any intermediate ones). If called
     644             :    before PEPSolve() and after PEPSetUp() then only the information of the
     645             :    endpoints of subintervals is available.
     646             : 
     647             :    This function is only available for spectrum slicing runs.
     648             : 
     649             :    The returned arrays should be freed by the user. Can pass NULL in any of
     650             :    the two arrays if not required.
     651             : 
     652             :    Fortran Notes:
     653             :    The calling sequence from Fortran is
     654             : .vb
     655             :    PEPSTOARGetInertias(pep,n,shifts,inertias,ierr)
     656             :    integer n
     657             :    double precision shifts(*)
     658             :    integer inertias(*)
     659             : .ve
     660             :    The arrays should be at least of length n. The value of n can be determined
     661             :    by an initial call
     662             : .vb
     663             :    PEPSTOARGetInertias(pep,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
     664             : .ve
     665             : 
     666             :    Level: advanced
     667             : 
     668             : .seealso: PEPSetInterval()
     669             : @*/
     670           2 : PetscErrorCode PEPSTOARGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
     671             : {
     672           2 :   PetscFunctionBegin;
     673           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     674           2 :   PetscAssertPointer(n,2);
     675           2 :   PetscUseMethod(pep,"PEPSTOARGetInertias_C",(PEP,PetscInt*,PetscReal**,PetscInt**),(pep,n,shifts,inertias));
     676           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     677             : }
     678             : 
     679           1 : static PetscErrorCode PEPSTOARSetDetectZeros_STOAR(PEP pep,PetscBool detect)
     680             : {
     681           1 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     682             : 
     683           1 :   PetscFunctionBegin;
     684           1 :   ctx->detect = detect;
     685           1 :   pep->state  = PEP_STATE_INITIAL;
     686           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     687             : }
     688             : 
     689             : /*@
     690             :    PEPSTOARSetDetectZeros - Sets a flag to enforce detection of
     691             :    zeros during the factorizations throughout the spectrum slicing computation.
     692             : 
     693             :    Logically Collective
     694             : 
     695             :    Input Parameters:
     696             : +  pep    - the eigenproblem solver context
     697             : -  detect - check for zeros
     698             : 
     699             :    Options Database Key:
     700             : .  -pep_stoar_detect_zeros - Check for zeros; this takes an optional
     701             :    bool value (0/1/no/yes/true/false)
     702             : 
     703             :    Notes:
     704             :    A zero in the factorization indicates that a shift coincides with an eigenvalue.
     705             : 
     706             :    This flag is turned off by default, and may be necessary in some cases.
     707             :    This feature currently requires an external package for factorizations
     708             :    with support for zero detection, e.g. MUMPS.
     709             : 
     710             :    Level: advanced
     711             : 
     712             : .seealso: PEPSetInterval()
     713             : @*/
     714           1 : PetscErrorCode PEPSTOARSetDetectZeros(PEP pep,PetscBool detect)
     715             : {
     716           1 :   PetscFunctionBegin;
     717           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     718           4 :   PetscValidLogicalCollectiveBool(pep,detect,2);
     719           1 :   PetscTryMethod(pep,"PEPSTOARSetDetectZeros_C",(PEP,PetscBool),(pep,detect));
     720           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     721             : }
     722             : 
     723           2 : static PetscErrorCode PEPSTOARGetDetectZeros_STOAR(PEP pep,PetscBool *detect)
     724             : {
     725           2 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     726             : 
     727           2 :   PetscFunctionBegin;
     728           2 :   *detect = ctx->detect;
     729           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     730             : }
     731             : 
     732             : /*@
     733             :    PEPSTOARGetDetectZeros - Gets the flag that enforces zero detection
     734             :    in spectrum slicing.
     735             : 
     736             :    Not Collective
     737             : 
     738             :    Input Parameter:
     739             : .  pep - the eigenproblem solver context
     740             : 
     741             :    Output Parameter:
     742             : .  detect - whether zeros detection is enforced during factorizations
     743             : 
     744             :    Level: advanced
     745             : 
     746             : .seealso: PEPSTOARSetDetectZeros()
     747             : @*/
     748           2 : PetscErrorCode PEPSTOARGetDetectZeros(PEP pep,PetscBool *detect)
     749             : {
     750           2 :   PetscFunctionBegin;
     751           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     752           2 :   PetscAssertPointer(detect,2);
     753           2 :   PetscUseMethod(pep,"PEPSTOARGetDetectZeros_C",(PEP,PetscBool*),(pep,detect));
     754           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     755             : }
     756             : 
     757           4 : static PetscErrorCode PEPSTOARSetLinearization_STOAR(PEP pep,PetscReal alpha,PetscReal beta)
     758             : {
     759           4 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     760             : 
     761           4 :   PetscFunctionBegin;
     762           4 :   PetscCheck(beta!=0.0 || alpha!=0.0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Parameters alpha and beta cannot be zero simultaneously");
     763           4 :   ctx->alpha = alpha;
     764           4 :   ctx->beta  = beta;
     765           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     766             : }
     767             : 
     768             : /*@
     769             :    PEPSTOARSetLinearization - Set the coefficients that define
     770             :    the linearization of a quadratic eigenproblem.
     771             : 
     772             :    Logically Collective
     773             : 
     774             :    Input Parameters:
     775             : +  pep   - polynomial eigenvalue solver
     776             : .  alpha - first parameter of the linearization
     777             : -  beta  - second parameter of the linearization
     778             : 
     779             :    Options Database Key:
     780             : .  -pep_stoar_linearization <alpha,beta> - Sets the coefficients
     781             : 
     782             :    Notes:
     783             :    Cannot pass zero for both alpha and beta. The default values are
     784             :    alpha=1 and beta=0.
     785             : 
     786             :    Level: advanced
     787             : 
     788             : .seealso: PEPSTOARGetLinearization()
     789             : @*/
     790           4 : PetscErrorCode PEPSTOARSetLinearization(PEP pep,PetscReal alpha,PetscReal beta)
     791             : {
     792           4 :   PetscFunctionBegin;
     793           4 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     794          16 :   PetscValidLogicalCollectiveReal(pep,alpha,2);
     795          16 :   PetscValidLogicalCollectiveReal(pep,beta,3);
     796           4 :   PetscTryMethod(pep,"PEPSTOARSetLinearization_C",(PEP,PetscReal,PetscReal),(pep,alpha,beta));
     797           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     798             : }
     799             : 
     800           1 : static PetscErrorCode PEPSTOARGetLinearization_STOAR(PEP pep,PetscReal *alpha,PetscReal *beta)
     801             : {
     802           1 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     803             : 
     804           1 :   PetscFunctionBegin;
     805           1 :   if (alpha) *alpha = ctx->alpha;
     806           1 :   if (beta)  *beta  = ctx->beta;
     807           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     808             : }
     809             : 
     810             : /*@
     811             :    PEPSTOARGetLinearization - Returns the coefficients that define
     812             :    the linearization of a quadratic eigenproblem.
     813             : 
     814             :    Not Collective
     815             : 
     816             :    Input Parameter:
     817             : .  pep  - polynomial eigenvalue solver
     818             : 
     819             :    Output Parameters:
     820             : +  alpha - the first parameter of the linearization
     821             : -  beta  - the second parameter of the linearization
     822             : 
     823             :    Level: advanced
     824             : 
     825             : .seealso: PEPSTOARSetLinearization()
     826             : @*/
     827           1 : PetscErrorCode PEPSTOARGetLinearization(PEP pep,PetscReal *alpha,PetscReal *beta)
     828             : {
     829           1 :   PetscFunctionBegin;
     830           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     831           1 :   PetscUseMethod(pep,"PEPSTOARGetLinearization_C",(PEP,PetscReal*,PetscReal*),(pep,alpha,beta));
     832           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     833             : }
     834             : 
     835           3 : static PetscErrorCode PEPSTOARSetDimensions_STOAR(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
     836             : {
     837           3 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     838             : 
     839           3 :   PetscFunctionBegin;
     840           3 :   PetscCheck(nev>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
     841           3 :   ctx->nev = nev;
     842           3 :   if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
     843           2 :     ctx->ncv = PETSC_DEFAULT;
     844             :   } else {
     845           1 :     PetscCheck(ncv>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
     846           1 :     ctx->ncv = ncv;
     847             :   }
     848           3 :   if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
     849           2 :     ctx->mpd = PETSC_DEFAULT;
     850             :   } else {
     851           1 :     PetscCheck(mpd>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
     852           1 :     ctx->mpd = mpd;
     853             :   }
     854           3 :   pep->state = PEP_STATE_INITIAL;
     855           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     856             : }
     857             : 
     858             : /*@
     859             :    PEPSTOARSetDimensions - Sets the dimensions used for each subsolve
     860             :    step in case of doing spectrum slicing for a computational interval.
     861             :    The meaning of the parameters is the same as in PEPSetDimensions().
     862             : 
     863             :    Logically Collective
     864             : 
     865             :    Input Parameters:
     866             : +  pep - the eigenproblem solver context
     867             : .  nev - number of eigenvalues to compute
     868             : .  ncv - the maximum dimension of the subspace to be used by the subsolve
     869             : -  mpd - the maximum dimension allowed for the projected problem
     870             : 
     871             :    Options Database Key:
     872             : +  -eps_stoar_nev <nev> - Sets the number of eigenvalues
     873             : .  -eps_stoar_ncv <ncv> - Sets the dimension of the subspace
     874             : -  -eps_stoar_mpd <mpd> - Sets the maximum projected dimension
     875             : 
     876             :    Level: advanced
     877             : 
     878             : .seealso: PEPSTOARGetDimensions(), PEPSetDimensions(), PEPSetInterval()
     879             : @*/
     880           3 : PetscErrorCode PEPSTOARSetDimensions(PEP pep,PetscInt nev,PetscInt ncv,PetscInt mpd)
     881             : {
     882           3 :   PetscFunctionBegin;
     883           3 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     884          12 :   PetscValidLogicalCollectiveInt(pep,nev,2);
     885          12 :   PetscValidLogicalCollectiveInt(pep,ncv,3);
     886          12 :   PetscValidLogicalCollectiveInt(pep,mpd,4);
     887           3 :   PetscTryMethod(pep,"PEPSTOARSetDimensions_C",(PEP,PetscInt,PetscInt,PetscInt),(pep,nev,ncv,mpd));
     888           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     889             : }
     890             : 
     891           2 : static PetscErrorCode PEPSTOARGetDimensions_STOAR(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
     892             : {
     893           2 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     894             : 
     895           2 :   PetscFunctionBegin;
     896           2 :   if (nev) *nev = ctx->nev;
     897           2 :   if (ncv) *ncv = ctx->ncv;
     898           2 :   if (mpd) *mpd = ctx->mpd;
     899           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     900             : }
     901             : 
     902             : /*@
     903             :    PEPSTOARGetDimensions - Gets the dimensions used for each subsolve
     904             :    step in case of doing spectrum slicing for a computational interval.
     905             : 
     906             :    Not Collective
     907             : 
     908             :    Input Parameter:
     909             : .  pep - the eigenproblem solver context
     910             : 
     911             :    Output Parameters:
     912             : +  nev - number of eigenvalues to compute
     913             : .  ncv - the maximum dimension of the subspace to be used by the subsolve
     914             : -  mpd - the maximum dimension allowed for the projected problem
     915             : 
     916             :    Level: advanced
     917             : 
     918             : .seealso: PEPSTOARSetDimensions()
     919             : @*/
     920           2 : PetscErrorCode PEPSTOARGetDimensions(PEP pep,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
     921             : {
     922           2 :   PetscFunctionBegin;
     923           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     924           2 :   PetscUseMethod(pep,"PEPSTOARGetDimensions_C",(PEP,PetscInt*,PetscInt*,PetscInt*),(pep,nev,ncv,mpd));
     925           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     926             : }
     927             : 
     928           1 : static PetscErrorCode PEPSTOARSetCheckEigenvalueType_STOAR(PEP pep,PetscBool checket)
     929             : {
     930           1 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     931             : 
     932           1 :   PetscFunctionBegin;
     933           1 :   ctx->checket = checket;
     934           1 :   pep->state   = PEP_STATE_INITIAL;
     935           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     936             : }
     937             : 
     938             : /*@
     939             :    PEPSTOARSetCheckEigenvalueType - Sets a flag to check that all the eigenvalues
     940             :    obtained throughout the spectrum slicing computation have the same definite type.
     941             : 
     942             :    Logically Collective
     943             : 
     944             :    Input Parameters:
     945             : +  pep     - the eigenproblem solver context
     946             : -  checket - check eigenvalue type
     947             : 
     948             :    Options Database Key:
     949             : .  -pep_stoar_check_eigenvalue_type - Check eigenvalue type; this takes an optional
     950             :    bool value (0/1/no/yes/true/false)
     951             : 
     952             :    Notes:
     953             :    This option is relevant only for spectrum slicing computations, but it is
     954             :    ignored if the problem type is PEP_HYPERBOLIC.
     955             : 
     956             :    This flag is turned on by default, to guarantee that the computed eigenvalues
     957             :    have the same type (otherwise the computed solution might be wrong). But since
     958             :    the check is computationally quite expensive, the check may be turned off if
     959             :    the user knows for sure that all eigenvalues in the requested interval have
     960             :    the same type.
     961             : 
     962             :    Level: advanced
     963             : 
     964             : .seealso: PEPSetProblemType(), PEPSetInterval()
     965             : @*/
     966           1 : PetscErrorCode PEPSTOARSetCheckEigenvalueType(PEP pep,PetscBool checket)
     967             : {
     968           1 :   PetscFunctionBegin;
     969           1 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
     970           4 :   PetscValidLogicalCollectiveBool(pep,checket,2);
     971           1 :   PetscTryMethod(pep,"PEPSTOARSetCheckEigenvalueType_C",(PEP,PetscBool),(pep,checket));
     972           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     973             : }
     974             : 
     975           2 : static PetscErrorCode PEPSTOARGetCheckEigenvalueType_STOAR(PEP pep,PetscBool *checket)
     976             : {
     977           2 :   PEP_STOAR *ctx = (PEP_STOAR*)pep->data;
     978             : 
     979           2 :   PetscFunctionBegin;
     980           2 :   *checket = ctx->checket;
     981           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     982             : }
     983             : 
     984             : /*@
     985             :    PEPSTOARGetCheckEigenvalueType - Gets the flag for the eigenvalue type
     986             :    check in spectrum slicing.
     987             : 
     988             :    Not Collective
     989             : 
     990             :    Input Parameter:
     991             : .  pep - the eigenproblem solver context
     992             : 
     993             :    Output Parameter:
     994             : .  checket - whether eigenvalue type must be checked during spectrum slcing
     995             : 
     996             :    Level: advanced
     997             : 
     998             : .seealso: PEPSTOARSetCheckEigenvalueType()
     999             : @*/
    1000           2 : PetscErrorCode PEPSTOARGetCheckEigenvalueType(PEP pep,PetscBool *checket)
    1001             : {
    1002           2 :   PetscFunctionBegin;
    1003           2 :   PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
    1004           2 :   PetscAssertPointer(checket,2);
    1005           2 :   PetscUseMethod(pep,"PEPSTOARGetCheckEigenvalueType_C",(PEP,PetscBool*),(pep,checket));
    1006           2 :   PetscFunctionReturn(PETSC_SUCCESS);
    1007             : }
    1008             : 
    1009           0 : static PetscErrorCode PEPView_STOAR(PEP pep,PetscViewer viewer)
    1010             : {
    1011           0 :   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
    1012           0 :   PetscBool      isascii;
    1013             : 
    1014           0 :   PetscFunctionBegin;
    1015           0 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
    1016           0 :   if (isascii) {
    1017           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-"));
    1018           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"  linearization parameters: alpha=%g beta=%g\n",(double)ctx->alpha,(double)ctx->beta));
    1019           0 :     if (pep->which==PEP_ALL && !ctx->hyperbolic) PetscCall(PetscViewerASCIIPrintf(viewer,"  checking eigenvalue type: %s\n",ctx->checket?"enabled":"disabled"));
    1020             :   }
    1021           0 :   PetscFunctionReturn(PETSC_SUCCESS);
    1022             : }
    1023             : 
    1024          21 : static PetscErrorCode PEPReset_STOAR(PEP pep)
    1025             : {
    1026          21 :   PetscFunctionBegin;
    1027          21 :   if (pep->which==PEP_ALL) PetscCall(PEPReset_STOAR_QSlice(pep));
    1028          21 :   PetscFunctionReturn(PETSC_SUCCESS);
    1029             : }
    1030             : 
    1031          20 : static PetscErrorCode PEPDestroy_STOAR(PEP pep)
    1032             : {
    1033          20 :   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
    1034             : 
    1035          20 :   PetscFunctionBegin;
    1036          20 :   PetscCall(BVDestroy(&ctx->V));
    1037          20 :   PetscCall(PetscFree(pep->data));
    1038          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",NULL));
    1039          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",NULL));
    1040          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",NULL));
    1041          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",NULL));
    1042          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",NULL));
    1043          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",NULL));
    1044          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",NULL));
    1045          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",NULL));
    1046          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",NULL));
    1047          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",NULL));
    1048          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",NULL));
    1049          20 :   PetscFunctionReturn(PETSC_SUCCESS);
    1050             : }
    1051             : 
    1052          20 : SLEPC_EXTERN PetscErrorCode PEPCreate_STOAR(PEP pep)
    1053             : {
    1054          20 :   PEP_STOAR      *ctx;
    1055             : 
    1056          20 :   PetscFunctionBegin;
    1057          20 :   PetscCall(PetscNew(&ctx));
    1058          20 :   pep->data = (void*)ctx;
    1059             : 
    1060          20 :   pep->lineariz = PETSC_TRUE;
    1061          20 :   ctx->lock     = PETSC_TRUE;
    1062          20 :   ctx->nev      = 1;
    1063          20 :   ctx->ncv      = PETSC_DEFAULT;
    1064          20 :   ctx->mpd      = PETSC_DEFAULT;
    1065          20 :   ctx->alpha    = 1.0;
    1066          20 :   ctx->beta     = 0.0;
    1067          20 :   ctx->checket  = PETSC_TRUE;
    1068             : 
    1069          20 :   pep->ops->setup          = PEPSetUp_STOAR;
    1070          20 :   pep->ops->setfromoptions = PEPSetFromOptions_STOAR;
    1071          20 :   pep->ops->destroy        = PEPDestroy_STOAR;
    1072          20 :   pep->ops->view           = PEPView_STOAR;
    1073          20 :   pep->ops->backtransform  = PEPBackTransform_Default;
    1074          20 :   pep->ops->computevectors = PEPComputeVectors_Default;
    1075          20 :   pep->ops->extractvectors = PEPExtractVectors_TOAR;
    1076          20 :   pep->ops->reset          = PEPReset_STOAR;
    1077             : 
    1078          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLocking_C",PEPSTOARSetLocking_STOAR));
    1079          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLocking_C",PEPSTOARGetLocking_STOAR));
    1080          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDetectZeros_C",PEPSTOARSetDetectZeros_STOAR));
    1081          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDetectZeros_C",PEPSTOARGetDetectZeros_STOAR));
    1082          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetInertias_C",PEPSTOARGetInertias_STOAR));
    1083          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetDimensions_C",PEPSTOARGetDimensions_STOAR));
    1084          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetDimensions_C",PEPSTOARSetDimensions_STOAR));
    1085          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetLinearization_C",PEPSTOARSetLinearization_STOAR));
    1086          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetLinearization_C",PEPSTOARGetLinearization_STOAR));
    1087          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARSetCheckEigenvalueType_C",PEPSTOARSetCheckEigenvalueType_STOAR));
    1088          20 :   PetscCall(PetscObjectComposeFunction((PetscObject)pep,"PEPSTOARGetCheckEigenvalueType_C",PEPSTOARGetCheckEigenvalueType_STOAR));
    1089          20 :   PetscFunctionReturn(PETSC_SUCCESS);
    1090             : }

Generated by: LCOV version 1.14