LCOV - code coverage report
Current view: top level - pep/impls/krylov/stoar - stoar.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 510 551 92.6 %
Date: 2020-10-22 06:47:52 Functions: 32 33 97.0 %

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

Generated by: LCOV version 1.13