LCOV - code coverage report
Current view: top level - pep/impls/krylov/stoar - stoar.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 513 554 92.6 %
Date: 2019-11-15 09:58:42 Functions: 32 33 97.0 %

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

Generated by: LCOV version 1.13