LCOV - code coverage report
Current view: top level - pep/impls/krylov/stoar - qslice.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 957 1039 92.1 %
Date: 2024-03-28 00:28:38 Functions: 21 22 95.5 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       6             :    This file is part of SLEPc.
       7             :    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
       8             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       9             : */
      10             : /*
      11             :    SLEPc polynomial eigensolver: "stoar"
      12             : 
      13             :    Method: S-TOAR with spectrum slicing for symmetric quadratic eigenproblems
      14             : 
      15             :    Algorithm:
      16             : 
      17             :        Symmetric Two-Level Orthogonal Arnoldi.
      18             : 
      19             :    References:
      20             : 
      21             :        [1] C. Campos and J.E. Roman, "Inertia-based spectrum slicing
      22             :            for symmetric quadratic eigenvalue problems", Numer. Linear
      23             :            Algebra Appl. 27(4):e2293, 2020.
      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-slice-qep,\n"
      33             :   "   author = \"C. Campos and J. E. Roman\",\n"
      34             :   "   title = \"Inertia-based spectrum slicing for symmetric quadratic eigenvalue problems\",\n"
      35             :   "   journal = \"Numer. Linear Algebra Appl.\",\n"
      36             :   "   volume = \"27\",\n"
      37             :   "   number = \"4\",\n"
      38             :   "   pages = \"e2293\",\n"
      39             :   "   year = \"2020,\"\n"
      40             :   "   doi = \"https://doi.org/10.1002/nla.2293\"\n"
      41             :   "}\n";
      42             : 
      43             : #define SLICE_PTOL PETSC_SQRT_MACHINE_EPSILON
      44             : 
      45          16 : static PetscErrorCode PEPQSliceResetSR(PEP pep)
      46             : {
      47          16 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
      48          16 :   PEP_SR         sr=ctx->sr;
      49          16 :   PEP_shift      s;
      50          16 :   PetscInt       i;
      51             : 
      52          16 :   PetscFunctionBegin;
      53          16 :   if (sr) {
      54             :     /* Reviewing list of shifts to free memory */
      55           8 :     s = sr->s0;
      56           8 :     if (s) {
      57          21 :       while (s->neighb[1]) {
      58          14 :         s = s->neighb[1];
      59          21 :         PetscCall(PetscFree(s->neighb[0]));
      60             :       }
      61           7 :       PetscCall(PetscFree(s));
      62             :     }
      63           8 :     PetscCall(PetscFree(sr->S));
      64         177 :     for (i=0;i<pep->nconv;i++) PetscCall(PetscFree(sr->qinfo[i].q));
      65           8 :     PetscCall(PetscFree(sr->qinfo));
      66          32 :     for (i=0;i<3;i++) PetscCall(VecDestroy(&sr->v[i]));
      67           8 :     PetscCall(EPSDestroy(&sr->eps));
      68           8 :     PetscCall(PetscFree(sr));
      69             :   }
      70          16 :   ctx->sr = NULL;
      71          16 :   PetscFunctionReturn(PETSC_SUCCESS);
      72             : }
      73             : 
      74           8 : PetscErrorCode PEPReset_STOAR_QSlice(PEP pep)
      75             : {
      76           8 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
      77             : 
      78           8 :   PetscFunctionBegin;
      79           8 :   PetscCall(PEPQSliceResetSR(pep));
      80           8 :   PetscCall(PetscFree(ctx->inertias));
      81           8 :   PetscCall(PetscFree(ctx->shifts));
      82           8 :   PetscFunctionReturn(PETSC_SUCCESS);
      83             : }
      84             : 
      85             : /*
      86             :   PEPQSliceAllocateSolution - Allocate memory storage for common variables such
      87             :   as eigenvalues and eigenvectors.
      88             : */
      89           7 : static PetscErrorCode PEPQSliceAllocateSolution(PEP pep)
      90             : {
      91           7 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
      92           7 :   PetscInt       k;
      93           7 :   BVType         type;
      94           7 :   Vec            t;
      95           7 :   PEP_SR         sr = ctx->sr;
      96             : 
      97           7 :   PetscFunctionBegin;
      98             :   /* allocate space for eigenvalues and friends */
      99           7 :   k = PetscMax(1,sr->numEigs);
     100           7 :   PetscCall(PetscFree4(sr->eigr,sr->eigi,sr->errest,sr->perm));
     101           7 :   PetscCall(PetscCalloc4(k,&sr->eigr,k,&sr->eigi,k,&sr->errest,k,&sr->perm));
     102           7 :   PetscCall(PetscFree(sr->qinfo));
     103           7 :   PetscCall(PetscCalloc1(k,&sr->qinfo));
     104             : 
     105             :   /* allocate sr->V and transfer options from pep->V */
     106           7 :   PetscCall(BVDestroy(&sr->V));
     107           7 :   PetscCall(BVCreate(PetscObjectComm((PetscObject)pep),&sr->V));
     108           7 :   if (!pep->V) PetscCall(PEPGetBV(pep,&pep->V));
     109           7 :   if (!((PetscObject)pep->V)->type_name) PetscCall(BVSetType(sr->V,BVMAT));
     110             :   else {
     111           7 :     PetscCall(BVGetType(pep->V,&type));
     112           7 :     PetscCall(BVSetType(sr->V,type));
     113             :   }
     114           7 :   PetscCall(STMatCreateVecsEmpty(pep->st,&t,NULL));
     115           7 :   PetscCall(BVSetSizesFromVec(sr->V,t,k+1));
     116           7 :   PetscCall(VecDestroy(&t));
     117           7 :   sr->ld = k;
     118           7 :   PetscCall(PetscFree(sr->S));
     119           7 :   PetscCall(PetscMalloc1((k+1)*sr->ld*(pep->nmat-1),&sr->S));
     120           7 :   PetscFunctionReturn(PETSC_SUCCESS);
     121             : }
     122             : 
     123             : /* Convergence test to compute positive Ritz values */
     124          31 : static PetscErrorCode ConvergedPositive(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
     125             : {
     126          31 :   PetscFunctionBegin;
     127          31 :   *errest = (PetscRealPart(eigr)>0.0)?0.0:res;
     128          31 :   PetscFunctionReturn(PETSC_SUCCESS);
     129             : }
     130             : 
     131          33 : static PetscErrorCode PEPQSliceMatGetInertia(PEP pep,PetscReal shift,PetscInt *inertia,PetscInt *zeros)
     132             : {
     133          33 :   KSP            ksp,kspr;
     134          33 :   PC             pc;
     135          33 :   Mat            F;
     136          33 :   PetscBool      flg;
     137             : 
     138          33 :   PetscFunctionBegin;
     139          33 :   if (!pep->solvematcoeffs) PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
     140          33 :   if (shift==PETSC_MAX_REAL) { /* Inertia of matrix A[2] */
     141           0 :     pep->solvematcoeffs[0] = 0.0; pep->solvematcoeffs[1] = 0.0; pep->solvematcoeffs[2] = 1.0;
     142          33 :   } else PetscCall(PEPEvaluateBasis(pep,shift,0,pep->solvematcoeffs,NULL));
     143          33 :   PetscCall(STMatSetUp(pep->st,pep->sfactor,pep->solvematcoeffs));
     144          33 :   PetscCall(STGetKSP(pep->st,&ksp));
     145          33 :   PetscCall(KSPGetPC(ksp,&pc));
     146          33 :   PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCREDUNDANT,&flg));
     147          33 :   if (flg) {
     148           0 :     PetscCall(PCRedundantGetKSP(pc,&kspr));
     149           0 :     PetscCall(KSPGetPC(kspr,&pc));
     150             :   }
     151          33 :   PetscCall(PCFactorGetMatrix(pc,&F));
     152          33 :   PetscCall(MatGetInertia(F,inertia,zeros,NULL));
     153          33 :   PetscFunctionReturn(PETSC_SUCCESS);
     154             : }
     155             : 
     156          31 : static PetscErrorCode PEPQSliceGetInertia(PEP pep,PetscReal shift,PetscInt *inertia,PetscInt *zeros,PetscInt correction)
     157             : {
     158          31 :   KSP            ksp;
     159          31 :   Mat            P;
     160          31 :   PetscReal      nzshift=0.0,dot;
     161          31 :   PetscRandom    rand;
     162          31 :   PetscInt       nconv;
     163          31 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
     164          31 :   PEP_SR         sr=ctx->sr;
     165             : 
     166          31 :   PetscFunctionBegin;
     167          31 :   if (shift >= PETSC_MAX_REAL) { /* Right-open interval */
     168           0 :     *inertia = 0;
     169          31 :   } else if (shift <= PETSC_MIN_REAL) {
     170           0 :     *inertia = 0;
     171           0 :     if (zeros) *zeros = 0;
     172             :   } else {
     173             :     /* If the shift is zero, perturb it to a very small positive value.
     174             :        The goal is that the nonzero pattern is the same in all cases and reuse
     175             :        the symbolic factorizations */
     176          31 :     nzshift = (shift==0.0)? 10.0/PETSC_MAX_REAL: shift;
     177          31 :     PetscCall(PEPQSliceMatGetInertia(pep,nzshift,inertia,zeros));
     178          31 :     PetscCall(STSetShift(pep->st,nzshift));
     179             :   }
     180          31 :   if (!correction) {
     181          17 :     if (shift >= PETSC_MAX_REAL) *inertia = 2*pep->n;
     182          17 :     else if (shift>PETSC_MIN_REAL) {
     183          17 :       PetscCall(STGetKSP(pep->st,&ksp));
     184          17 :       PetscCall(KSPGetOperators(ksp,&P,NULL));
     185          17 :       if (*inertia!=pep->n && !sr->v[0]) {
     186           6 :         PetscCall(MatCreateVecs(P,&sr->v[0],NULL));
     187           6 :         PetscCall(VecDuplicate(sr->v[0],&sr->v[1]));
     188           6 :         PetscCall(VecDuplicate(sr->v[0],&sr->v[2]));
     189           6 :         PetscCall(BVGetRandomContext(pep->V,&rand));
     190           6 :         PetscCall(VecSetRandom(sr->v[0],rand));
     191             :       }
     192          17 :       if (*inertia<pep->n && *inertia>0) {
     193           9 :         if (!sr->eps) {
     194           5 :           PetscCall(EPSCreate(PetscObjectComm((PetscObject)pep),&sr->eps));
     195           5 :           PetscCall(EPSSetProblemType(sr->eps,EPS_HEP));
     196           5 :           PetscCall(EPSSetWhichEigenpairs(sr->eps,EPS_LARGEST_REAL));
     197             :         }
     198           9 :         PetscCall(EPSSetConvergenceTestFunction(sr->eps,ConvergedPositive,NULL,NULL));
     199           9 :         PetscCall(EPSSetOperators(sr->eps,P,NULL));
     200           9 :         PetscCall(EPSSolve(sr->eps));
     201           9 :         PetscCall(EPSGetConverged(sr->eps,&nconv));
     202           9 :         PetscCheck(nconv,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)nzshift);
     203           9 :         PetscCall(EPSGetEigenpair(sr->eps,0,NULL,NULL,sr->v[0],sr->v[1]));
     204             :       }
     205          17 :       if (*inertia!=pep->n) {
     206          12 :         PetscCall(MatMult(pep->A[1],sr->v[0],sr->v[1]));
     207          12 :         PetscCall(MatMult(pep->A[2],sr->v[0],sr->v[2]));
     208          12 :         PetscCall(VecAXPY(sr->v[1],2*nzshift,sr->v[2]));
     209          12 :         PetscCall(VecDotRealPart(sr->v[1],sr->v[0],&dot));
     210          12 :         if (dot>0.0) *inertia = 2*pep->n-*inertia;
     211             :       }
     212             :     }
     213          14 :   } else if (correction<0) *inertia = 2*pep->n-*inertia;
     214          31 :   PetscFunctionReturn(PETSC_SUCCESS);
     215             : }
     216             : 
     217             : /*
     218             :    Check eigenvalue type - used only in non-hyperbolic problems.
     219             :    All computed eigenvalues must have the same definite type (positive or negative).
     220             :    If ini=TRUE the type is available in omega, otherwise we compute an eigenvalue
     221             :    closest to shift and determine its type.
     222             :  */
     223         279 : static PetscErrorCode PEPQSliceCheckEigenvalueType(PEP pep,PetscReal shift,PetscReal omega,PetscBool ini)
     224             : {
     225         279 :   PEP            pep2;
     226         279 :   ST             st;
     227         279 :   PetscInt       nconv;
     228         279 :   PetscScalar    lambda;
     229         279 :   PetscReal      dot;
     230         279 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
     231         279 :   PEP_SR         sr=ctx->sr;
     232             : 
     233         279 :   PetscFunctionBegin;
     234         279 :   if (!ini) {
     235         275 :     PetscCheck(-(omega/(shift*ctx->alpha+ctx->beta))*sr->type>=0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Different positive/negative type detected in eigenvalue %g",(double)shift);
     236             :   } else {
     237           4 :     PetscCall(PEPCreate(PetscObjectComm((PetscObject)pep),&pep2));
     238           4 :     PetscCall(PEPSetOptionsPrefix(pep2,((PetscObject)pep)->prefix));
     239           4 :     PetscCall(PEPAppendOptionsPrefix(pep2,"pep_eigenvalue_type_"));
     240           4 :     PetscCall(PEPSetTolerances(pep2,PETSC_DEFAULT,pep->max_it/4));
     241           4 :     PetscCall(PEPSetType(pep2,PEPTOAR));
     242           4 :     PetscCall(PEPSetOperators(pep2,pep->nmat,pep->A));
     243           4 :     PetscCall(PEPSetWhichEigenpairs(pep2,PEP_TARGET_MAGNITUDE));
     244           4 :     PetscCall(PEPGetRG(pep2,&pep2->rg));
     245           4 :     PetscCall(RGSetType(pep2->rg,RGINTERVAL));
     246             : #if defined(PETSC_USE_COMPLEX)
     247           4 :     PetscCall(RGIntervalSetEndpoints(pep2->rg,pep->inta,pep->intb,-PETSC_SQRT_MACHINE_EPSILON,PETSC_SQRT_MACHINE_EPSILON));
     248             : #else
     249             :     PetscCall(RGIntervalSetEndpoints(pep2->rg,pep->inta,pep->intb,0.0,0.0));
     250             : #endif
     251           4 :     pep2->target = shift;
     252           4 :     st = pep2->st;
     253           4 :     pep2->st = pep->st;
     254           4 :     PetscCall(PEPSolve(pep2));
     255           4 :     PetscCall(PEPGetConverged(pep2,&nconv));
     256           4 :     if (nconv) {
     257           4 :       PetscCall(PEPGetEigenpair(pep2,0,&lambda,NULL,pep2->work[0],NULL));
     258           4 :       PetscCall(MatMult(pep->A[1],pep2->work[0],pep2->work[1]));
     259           4 :       PetscCall(MatMult(pep->A[2],pep2->work[0],pep2->work[2]));
     260           4 :       PetscCall(VecAXPY(pep2->work[1],2.0*lambda*pep->sfactor,pep2->work[2]));
     261           4 :       PetscCall(VecDotRealPart(pep2->work[1],pep2->work[0],&dot));
     262           8 :       PetscCall(PetscInfo(pep,"lambda=%g, %s type\n",(double)PetscRealPart(lambda),(dot>0.0)?"positive":"negative"));
     263           6 :       if (!sr->type) sr->type = (dot>0.0)?1:-1;
     264           2 :       else PetscCheck(sr->type*dot>=0.0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Different positive/negative type detected in eigenvalue %g",(double)PetscRealPart(lambda));
     265             :     }
     266           4 :     pep2->st = st;
     267           4 :     PetscCall(PEPDestroy(&pep2));
     268             :   }
     269         279 :   PetscFunctionReturn(PETSC_SUCCESS);
     270             : }
     271             : 
     272           3 : static inline PetscErrorCode PEPQSliceDiscriminant(PEP pep,Vec u,Vec w,PetscReal *d,PetscReal *smas,PetscReal *smenos)
     273             : {
     274           3 :   PetscReal      ap,bp,cp,dis;
     275             : 
     276           3 :   PetscFunctionBegin;
     277           3 :   PetscCall(MatMult(pep->A[0],u,w));
     278           3 :   PetscCall(VecDotRealPart(w,u,&cp));
     279           3 :   PetscCall(MatMult(pep->A[1],u,w));
     280           3 :   PetscCall(VecDotRealPart(w,u,&bp));
     281           3 :   PetscCall(MatMult(pep->A[2],u,w));
     282           3 :   PetscCall(VecDotRealPart(w,u,&ap));
     283           3 :   dis = bp*bp-4*ap*cp;
     284           3 :   if (dis>=0.0 && smas) {
     285           3 :     if (ap>0) *smas = (-bp+PetscSqrtReal(dis))/(2*ap);
     286           3 :     else if (ap<0) *smas = (-bp-PetscSqrtReal(dis))/(2*ap);
     287             :     else {
     288           0 :       if (bp >0) *smas = -cp/bp;
     289           0 :       else *smas = PETSC_MAX_REAL;
     290             :     }
     291             :   }
     292           3 :   if (dis>=0.0 && smenos) {
     293           2 :     if (ap>0) *smenos = (-bp-PetscSqrtReal(dis))/(2*ap);
     294           2 :     else if (ap<0) *smenos = (-bp+PetscSqrtReal(dis))/(2*ap);
     295             :     else {
     296           0 :       if (bp<0) *smenos = -cp/bp;
     297           0 :       else *smenos = PETSC_MAX_REAL;
     298             :     }
     299             :   }
     300           3 :   if (d) *d = dis;
     301           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     302             : }
     303             : 
     304           2 : static inline PetscErrorCode PEPQSliceEvaluateQEP(PEP pep,PetscScalar x,Mat M,MatStructure str)
     305             : {
     306           2 :   PetscFunctionBegin;
     307           2 :   PetscCall(MatCopy(pep->A[0],M,SAME_NONZERO_PATTERN));
     308           2 :   PetscCall(MatAXPY(M,x,pep->A[1],str));
     309           2 :   PetscCall(MatAXPY(M,x*x,pep->A[2],str));
     310           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     311             : }
     312             : 
     313             : /*@
     314             :    PEPCheckDefiniteQEP - Determines if a symmetric/Hermitian quadratic eigenvalue problem
     315             :    is definite or not.
     316             : 
     317             :    Collective
     318             : 
     319             :    Input Parameter:
     320             : .  pep  - eigensolver context
     321             : 
     322             :    Output Parameters:
     323             : +  xi - first computed parameter
     324             : .  mu - second computed parameter
     325             : .  definite - flag indicating that the problem is definite
     326             : -  hyperbolic - flag indicating that the problem is hyperbolic
     327             : 
     328             :    Notes:
     329             :    This function is intended for quadratic eigenvalue problems, Q(lambda)=A*lambda^2+B*lambda+C,
     330             :    with symmetric (or Hermitian) coefficient matrices A,B,C.
     331             : 
     332             :    On output, the flag 'definite' may have the values -1 (meaning that the QEP is not
     333             :    definite), 1 (if the problem is definite), or 0 if the algorithm was not able to
     334             :    determine whether the problem is definite or not.
     335             : 
     336             :    If definite=1, the output flag 'hyperbolic' informs in a similar way about whether the
     337             :    problem is hyperbolic or not.
     338             : 
     339             :    If definite=1, the computed values xi and mu satisfy Q(xi)<0 and Q(mu)>0, as
     340             :    obtained via the method proposed in [Niendorf and Voss, LAA 2010]. Furthermore, if
     341             :    hyperbolic=1 then only xi is computed.
     342             : 
     343             :    Level: advanced
     344             : 
     345             : .seealso: PEPSetProblemType()
     346             : @*/
     347           1 : PetscErrorCode PEPCheckDefiniteQEP(PEP pep,PetscReal *xi,PetscReal *mu,PetscInt *definite,PetscInt *hyperbolic)
     348             : {
     349           1 :   PetscRandom    rand;
     350           1 :   Vec            u,w;
     351           1 :   PetscReal      d=0.0,s=0.0,sp,mut=0.0,omg=0.0,omgp;
     352           1 :   PetscInt       k,its=10,hyp=0,check=0,nconv,inertia,n;
     353           1 :   Mat            M=NULL;
     354           1 :   MatStructure   str;
     355           1 :   EPS            eps;
     356           1 :   PetscBool      transform,ptypehyp;
     357             : 
     358           1 :   PetscFunctionBegin;
     359           1 :   PetscCheck(pep->problem_type==PEP_HERMITIAN || pep->problem_type==PEP_HYPERBOLIC,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Only available for Hermitian (or hyperbolic) problems");
     360           1 :   ptypehyp = (pep->problem_type==PEP_HYPERBOLIC)? PETSC_TRUE: PETSC_FALSE;
     361           1 :   if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
     362           1 :   PetscCall(PEPSetDefaultST(pep));
     363           1 :   PetscCall(STSetMatrices(pep->st,pep->nmat,pep->A));
     364           1 :   PetscCall(MatGetSize(pep->A[0],&n,NULL));
     365           1 :   PetscCall(STGetTransform(pep->st,&transform));
     366           1 :   PetscCall(STSetTransform(pep->st,PETSC_FALSE));
     367           1 :   PetscCall(STSetUp(pep->st));
     368           1 :   PetscCall(MatCreateVecs(pep->A[0],&u,&w));
     369           1 :   PetscCall(PEPGetBV(pep,&pep->V));
     370           1 :   PetscCall(BVGetRandomContext(pep->V,&rand));
     371           1 :   PetscCall(VecSetRandom(u,rand));
     372           1 :   PetscCall(VecNormalize(u,NULL));
     373           1 :   PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,&s,NULL));
     374           1 :   if (d<0.0) check = -1;
     375           1 :   if (!check) {
     376           1 :     PetscCall(EPSCreate(PetscObjectComm((PetscObject)pep),&eps));
     377           1 :     PetscCall(EPSSetProblemType(eps,EPS_HEP));
     378           1 :     PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
     379           1 :     PetscCall(EPSSetTolerances(eps,PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON),PETSC_DECIDE));
     380           1 :     PetscCall(MatDuplicate(pep->A[0],MAT_DO_NOT_COPY_VALUES,&M));
     381           1 :     PetscCall(STGetMatStructure(pep->st,&str));
     382             :   }
     383           1 :   for (k=0;k<its&&!check;k++) {
     384           1 :     PetscCall(PEPQSliceEvaluateQEP(pep,s,M,str));
     385           1 :     PetscCall(EPSSetOperators(eps,M,NULL));
     386           1 :     PetscCall(EPSSolve(eps));
     387           1 :     PetscCall(EPSGetConverged(eps,&nconv));
     388           1 :     if (!nconv) break;
     389           1 :     PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
     390           1 :     sp = s;
     391           1 :     PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,&s,&omg));
     392           2 :     if (d<0.0) {check = -1; break;}
     393           1 :     if (PetscAbsReal((s-sp)/s)<100*PETSC_MACHINE_EPSILON) break;
     394           1 :     if (s>sp) {hyp = -1;}
     395           1 :     mut = 2*s-sp;
     396           1 :     PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
     397           1 :     if (inertia == n) {check = 1; break;}
     398             :   }
     399           1 :   for (;k<its&&!check;k++) {
     400           0 :     mut = (s-omg)/2;
     401           0 :     PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
     402           0 :     if (inertia == n) {check = 1; break;}
     403           0 :     if (PetscAbsReal((s-omg)/omg)<100*PETSC_MACHINE_EPSILON) break;
     404           0 :     PetscCall(PEPQSliceEvaluateQEP(pep,omg,M,str));
     405           0 :     PetscCall(EPSSetOperators(eps,M,NULL));
     406           0 :     PetscCall(EPSSolve(eps));
     407           0 :     PetscCall(EPSGetConverged(eps,&nconv));
     408           0 :     if (!nconv) break;
     409           0 :     PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
     410           0 :     omgp = omg;
     411           0 :     PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,NULL,&omg));
     412           0 :     if (d<0.0) {check = -1; break;}
     413           0 :     if (omg<omgp) hyp = -1;
     414             :   }
     415           1 :   if (check==1) *xi = mut;
     416           1 :   PetscCheck(hyp!=-1 || !ptypehyp,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Problem does not satisfy hyperbolic test; consider removing the hyperbolicity flag");
     417           1 :   if (check==1 && hyp==0) {
     418           0 :     PetscCall(PEPQSliceMatGetInertia(pep,PETSC_MAX_REAL,&inertia,NULL));
     419           0 :     if (inertia == 0) hyp = 1;
     420           0 :     else hyp = -1;
     421             :   }
     422           1 :   if (check==1 && hyp!=1) {
     423           1 :     check = 0;
     424           1 :     PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
     425           1 :     for (;k<its&&!check;k++) {
     426           1 :       PetscCall(PEPQSliceEvaluateQEP(pep,s,M,str));
     427           1 :       PetscCall(EPSSetOperators(eps,M,NULL));
     428           1 :       PetscCall(EPSSolve(eps));
     429           1 :       PetscCall(EPSGetConverged(eps,&nconv));
     430           1 :       if (!nconv) break;
     431           1 :       PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
     432           1 :       sp = s;
     433           1 :       PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,&s,&omg));
     434           2 :       if (d<0.0) {check = -1; break;}
     435           1 :       if (PetscAbsReal((s-sp)/s)<100*PETSC_MACHINE_EPSILON) break;
     436           1 :       mut = 2*s-sp;
     437           1 :       PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
     438           1 :       if (inertia == 0) {check = 1; break;}
     439             :     }
     440           1 :     for (;k<its&&!check;k++) {
     441           0 :       mut = (s-omg)/2;
     442           0 :       PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
     443           0 :       if (inertia == 0) {check = 1; break;}
     444           0 :       if (PetscAbsReal((s-omg)/omg)<100*PETSC_MACHINE_EPSILON) break;
     445           0 :       PetscCall(PEPQSliceEvaluateQEP(pep,omg,M,str));
     446           0 :       PetscCall(EPSSetOperators(eps,M,NULL));
     447           0 :       PetscCall(EPSSolve(eps));
     448           0 :       PetscCall(EPSGetConverged(eps,&nconv));
     449           0 :       if (!nconv) break;
     450           0 :       PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
     451           0 :       PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,NULL,&omg));
     452           0 :       if (d<0.0) {check = -1; break;}
     453             :     }
     454             :   }
     455           1 :   if (check==1) *mu = mut;
     456           1 :   *definite = check;
     457           1 :   *hyperbolic = hyp;
     458           1 :   if (M) PetscCall(MatDestroy(&M));
     459           1 :   PetscCall(VecDestroy(&u));
     460           1 :   PetscCall(VecDestroy(&w));
     461           1 :   PetscCall(EPSDestroy(&eps));
     462           1 :   PetscCall(STSetTransform(pep->st,transform));
     463           1 :   PetscFunctionReturn(PETSC_SUCCESS);
     464             : }
     465             : 
     466             : /*
     467             :    Dummy backtransform operation
     468             :  */
     469           0 : static PetscErrorCode PEPBackTransform_Skip(PEP pep)
     470             : {
     471           0 :   PetscFunctionBegin;
     472           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     473             : }
     474             : 
     475           8 : PetscErrorCode PEPSetUp_STOAR_QSlice(PEP pep)
     476             : {
     477           8 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
     478           8 :   PEP_SR         sr;
     479           8 :   PetscInt       ld,i,zeros=0;
     480           8 :   SlepcSC        sc;
     481           8 :   PetscReal      r;
     482             : 
     483           8 :   PetscFunctionBegin;
     484           8 :   PEPCheckSinvertCayley(pep);
     485           8 :   PetscCheck(pep->inta<pep->intb,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues unless you provide a computational interval with PEPSetInterval()");
     486           8 :   PetscCheck(pep->intb<PETSC_MAX_REAL || pep->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The defined computational interval should have at least one of their sides bounded");
     487           8 :   PEPCheckUnsupportedCondition(pep,PEP_FEATURE_STOPPING,PETSC_TRUE," (with spectrum slicing)");
     488           8 :   if (pep->tol==(PetscReal)PETSC_DEFAULT) {
     489             : #if defined(PETSC_USE_REAL_SINGLE)
     490             :     pep->tol = SLEPC_DEFAULT_TOL;
     491             : #else
     492             :     /* use tighter tolerance */
     493           8 :     pep->tol = SLEPC_DEFAULT_TOL*1e-2;
     494             : #endif
     495             :   }
     496           8 :   if (ctx->nev==1) ctx->nev = PetscMin(20,pep->n);  /* nev not set, use default value */
     497           8 :   PetscCheck(pep->n<=10 || ctx->nev>=10,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"nev cannot be less than 10 in spectrum slicing runs");
     498           8 :   pep->ops->backtransform = PEPBackTransform_Skip;
     499           8 :   if (pep->max_it==PETSC_DEFAULT) pep->max_it = 100;
     500             : 
     501             :   /* create spectrum slicing context and initialize it */
     502           8 :   PetscCall(PEPQSliceResetSR(pep));
     503           8 :   PetscCall(PetscNew(&sr));
     504           8 :   ctx->sr   = sr;
     505           8 :   sr->itsKs = 0;
     506           8 :   sr->nleap = 0;
     507           8 :   sr->sPres = NULL;
     508             : 
     509           8 :   if (pep->solvematcoeffs) PetscCall(PetscFree(pep->solvematcoeffs));
     510           8 :   PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
     511           8 :   if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
     512           8 :   PetscCall(STSetTransform(pep->st,PETSC_FALSE));
     513           8 :   PetscCall(STSetUp(pep->st));
     514             : 
     515           8 :   ctx->hyperbolic = (pep->problem_type==PEP_HYPERBOLIC)? PETSC_TRUE: PETSC_FALSE;
     516             : 
     517             :   /* check presence of ends and finding direction */
     518           8 :   if (pep->inta > PETSC_MIN_REAL || pep->intb >= PETSC_MAX_REAL) {
     519           8 :     sr->int0 = pep->inta;
     520           8 :     sr->int1 = pep->intb;
     521           8 :     sr->dir = 1;
     522           8 :     if (pep->intb >= PETSC_MAX_REAL) { /* Right-open interval */
     523           0 :       sr->hasEnd = PETSC_FALSE;
     524           8 :     } else sr->hasEnd = PETSC_TRUE;
     525             :   } else {
     526           0 :     sr->int0 = pep->intb;
     527           0 :     sr->int1 = pep->inta;
     528           0 :     sr->dir = -1;
     529           0 :     sr->hasEnd = PetscNot(pep->inta <= PETSC_MIN_REAL);
     530             :   }
     531             : 
     532             :   /* compute inertia0 */
     533          15 :   PetscCall(PEPQSliceGetInertia(pep,sr->int0,&sr->inertia0,ctx->detect?&zeros:NULL,ctx->hyperbolic?0:1));
     534           8 :   PetscCheck(!zeros || (sr->int0!=pep->inta && sr->int0!=pep->intb),((PetscObject)pep)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in the interval endpoint");
     535           8 :   if (!ctx->hyperbolic && ctx->checket) PetscCall(PEPQSliceCheckEigenvalueType(pep,sr->int0,0.0,PETSC_TRUE));
     536             : 
     537             :   /* compute inertia1 */
     538          15 :   PetscCall(PEPQSliceGetInertia(pep,sr->int1,&sr->inertia1,ctx->detect?&zeros:NULL,ctx->hyperbolic?0:1));
     539           8 :   PetscCheck(!zeros,((PetscObject)pep)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in an interval endpoint defined by user");
     540           8 :   if (!ctx->hyperbolic && ctx->checket && sr->hasEnd) {
     541           2 :     PetscCall(PEPQSliceCheckEigenvalueType(pep,sr->int1,0.0,PETSC_TRUE));
     542           2 :     PetscCheck(sr->type || sr->inertia1==sr->inertia0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"No information of eigenvalue type in interval");
     543           2 :     PetscCheck(!sr->type || sr->inertia1!=sr->inertia0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Different positive/negative type detected");
     544           2 :     if (sr->dir*(sr->inertia1-sr->inertia0)<0) {
     545           0 :       sr->intcorr = -1;
     546           0 :       sr->inertia0 = 2*pep->n-sr->inertia0;
     547           0 :       sr->inertia1 = 2*pep->n-sr->inertia1;
     548           2 :     } else sr->intcorr = 1;
     549             :   } else {
     550           6 :     if (sr->inertia0<=pep->n && sr->inertia1<=pep->n) sr->intcorr = 1;
     551           4 :     else if (sr->inertia0>=pep->n && sr->inertia1>=pep->n) sr->intcorr = -1;
     552             :   }
     553             : 
     554           8 :   if (sr->hasEnd) {
     555           8 :     sr->dir = -sr->dir; r = sr->int0; sr->int0 = sr->int1; sr->int1 = r;
     556           8 :     i = sr->inertia0; sr->inertia0 = sr->inertia1; sr->inertia1 = i;
     557             :   }
     558             : 
     559             :   /* number of eigenvalues in interval */
     560           8 :   sr->numEigs = (sr->dir)*(sr->inertia1 - sr->inertia0);
     561           8 :   PetscCall(PetscInfo(pep,"QSlice setup: allocating for %" PetscInt_FMT " eigenvalues in [%g,%g]\n",sr->numEigs,(double)pep->inta,(double)pep->intb));
     562           8 :   if (sr->numEigs) {
     563           7 :     PetscCall(PEPQSliceAllocateSolution(pep));
     564           7 :     PetscCall(PEPSetDimensions_Default(pep,ctx->nev,&ctx->ncv,&ctx->mpd));
     565           7 :     pep->nev = ctx->nev; pep->ncv = ctx->ncv; pep->mpd = ctx->mpd;
     566           7 :     ld   = ctx->ncv+2;
     567           7 :     PetscCall(DSSetType(pep->ds,DSGHIEP));
     568           7 :     PetscCall(DSSetCompact(pep->ds,PETSC_TRUE));
     569           7 :     PetscCall(DSSetExtraRow(pep->ds,PETSC_TRUE));
     570           7 :     PetscCall(DSAllocate(pep->ds,ld));
     571           7 :     PetscCall(DSGetSlepcSC(pep->ds,&sc));
     572           7 :     sc->rg            = NULL;
     573           7 :     sc->comparison    = SlepcCompareLargestMagnitude;
     574           7 :     sc->comparisonctx = NULL;
     575           7 :     sc->map           = NULL;
     576           7 :     sc->mapobj        = NULL;
     577           1 :   } else {pep->ncv = 0; pep->nev = 0; pep->mpd = 0;}
     578           8 :   PetscFunctionReturn(PETSC_SUCCESS);
     579             : }
     580             : 
     581             : /*
     582             :    Fills the fields of a shift structure
     583             : */
     584          23 : static PetscErrorCode PEPCreateShift(PEP pep,PetscReal val,PEP_shift neighb0,PEP_shift neighb1)
     585             : {
     586          23 :   PEP_shift      s,*pending2;
     587          23 :   PetscInt       i;
     588          23 :   PEP_SR         sr;
     589          23 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
     590             : 
     591          23 :   PetscFunctionBegin;
     592          23 :   sr = ctx->sr;
     593          23 :   PetscCall(PetscNew(&s));
     594          23 :   s->value = val;
     595          23 :   s->neighb[0] = neighb0;
     596          23 :   if (neighb0) neighb0->neighb[1] = s;
     597          23 :   s->neighb[1] = neighb1;
     598          23 :   if (neighb1) neighb1->neighb[0] = s;
     599          23 :   s->comp[0] = PETSC_FALSE;
     600          23 :   s->comp[1] = PETSC_FALSE;
     601          23 :   s->index = -1;
     602          23 :   s->neigs = 0;
     603          23 :   s->nconv[0] = s->nconv[1] = 0;
     604          23 :   s->nsch[0] = s->nsch[1]=0;
     605             :   /* Inserts in the stack of pending shifts */
     606             :   /* If needed, the array is resized */
     607          23 :   if (sr->nPend >= sr->maxPend) {
     608           0 :     sr->maxPend *= 2;
     609           0 :     PetscCall(PetscMalloc1(sr->maxPend,&pending2));
     610           0 :     for (i=0;i<sr->nPend;i++) pending2[i] = sr->pending[i];
     611           0 :     PetscCall(PetscFree(sr->pending));
     612           0 :     sr->pending = pending2;
     613             :   }
     614          23 :   sr->pending[sr->nPend++]=s;
     615          23 :   PetscFunctionReturn(PETSC_SUCCESS);
     616             : }
     617             : 
     618             : /* Provides next shift to be computed */
     619          22 : static PetscErrorCode PEPExtractShift(PEP pep)
     620             : {
     621          22 :   PetscInt       iner,zeros=0;
     622          22 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
     623          22 :   PEP_SR         sr;
     624          22 :   PetscReal      newShift,aux;
     625          22 :   PEP_shift      sPres;
     626             : 
     627          22 :   PetscFunctionBegin;
     628          22 :   sr = ctx->sr;
     629          22 :   if (sr->nPend > 0) {
     630          15 :     if (sr->dirch) {
     631           1 :       aux = sr->int1; sr->int1 = sr->int0; sr->int0 = aux;
     632           1 :       iner = sr->inertia1; sr->inertia1 = sr->inertia0; sr->inertia0 = iner;
     633           1 :       sr->dir *= -1;
     634           1 :       PetscCall(PetscFree(sr->s0->neighb[1]));
     635           1 :       PetscCall(PetscFree(sr->s0));
     636           1 :       sr->nPend--;
     637           1 :       PetscCall(PEPCreateShift(pep,sr->int0,NULL,NULL));
     638           1 :       sr->sPrev = NULL;
     639           1 :       sr->sPres = sr->pending[--sr->nPend];
     640           1 :       pep->target = sr->sPres->value;
     641           1 :       sr->s0 = sr->sPres;
     642           1 :       pep->reason = PEP_CONVERGED_ITERATING;
     643             :     } else {
     644          14 :       sr->sPrev = sr->sPres;
     645          14 :       sr->sPres = sr->pending[--sr->nPend];
     646             :     }
     647          15 :     sPres = sr->sPres;
     648          30 :     PetscCall(PEPQSliceGetInertia(pep,sPres->value,&iner,ctx->detect?&zeros:NULL,sr->intcorr));
     649          15 :     if (zeros) {
     650           0 :       newShift = sPres->value*(1.0+SLICE_PTOL);
     651           0 :       if (sr->dir*(sPres->neighb[0] && newShift-sPres->neighb[0]->value) < 0) newShift = (sPres->value+sPres->neighb[0]->value)/2;
     652           0 :       else if (sPres->neighb[1] && sr->dir*(sPres->neighb[1]->value-newShift) < 0) newShift = (sPres->value+sPres->neighb[1]->value)/2;
     653           0 :       PetscCall(PEPQSliceGetInertia(pep,newShift,&iner,&zeros,sr->intcorr));
     654           0 :       PetscCheck(!zeros,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)newShift);
     655           0 :       sPres->value = newShift;
     656             :     }
     657          15 :     sr->sPres->inertia = iner;
     658          15 :     pep->target = sr->sPres->value;
     659          15 :     pep->reason = PEP_CONVERGED_ITERATING;
     660          15 :     pep->its = 0;
     661           7 :   } else sr->sPres = NULL;
     662          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     663             : }
     664             : 
     665             : /*
     666             :   Obtains value of subsequent shift
     667             : */
     668          15 : static PetscErrorCode PEPGetNewShiftValue(PEP pep,PetscInt side,PetscReal *newS)
     669             : {
     670          15 :   PetscReal lambda,d_prev;
     671          15 :   PetscInt  i,idxP;
     672          15 :   PEP_SR    sr;
     673          15 :   PEP_shift sPres,s;
     674          15 :   PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
     675             : 
     676          15 :   PetscFunctionBegin;
     677          15 :   sr = ctx->sr;
     678          15 :   sPres = sr->sPres;
     679          15 :   if (sPres->neighb[side]) {
     680             :   /* Completing a previous interval */
     681           3 :     if (!sPres->neighb[side]->neighb[side] && sPres->neighb[side]->nconv[side]==0) { /* One of the ends might be too far from eigenvalues */
     682           0 :       if (side) *newS = (sPres->value + PetscRealPart(sr->eigr[sr->perm[sr->indexEig-1]]))/2;
     683           0 :       else *newS = (sPres->value + PetscRealPart(sr->eigr[sr->perm[0]]))/2;
     684           3 :     } else *newS=(sPres->value + sPres->neighb[side]->value)/2;
     685             :   } else { /* (Only for side=1). Creating a new interval. */
     686          12 :     if (sPres->neigs==0) {/* No value has been accepted*/
     687           5 :       if (sPres->neighb[0]) {
     688             :         /* Multiplying by 10 the previous distance */
     689           4 :         *newS = sPres->value + 10*sr->dir*PetscAbsReal(sPres->value - sPres->neighb[0]->value);
     690           4 :         sr->nleap++;
     691             :         /* Stops when the interval is open and no values are found in the last 5 shifts (there might be infinite eigenvalues) */
     692           4 :         PetscCheck(sr->hasEnd || sr->nleap<=5,PetscObjectComm((PetscObject)pep),PETSC_ERR_CONV_FAILED,"Unable to compute the wanted eigenvalues with open interval");
     693             :       } else { /* First shift */
     694           1 :         if (pep->nconv != 0) {
     695             :           /* Unaccepted values give information for next shift */
     696             :           idxP=0;/* Number of values left from shift */
     697           0 :           for (i=0;i<pep->nconv;i++) {
     698           0 :             lambda = PetscRealPart(pep->eigr[i]);
     699           0 :             if (sr->dir*(lambda - sPres->value) <0) idxP++;
     700             :             else break;
     701             :           }
     702             :           /* Avoiding subtraction of eigenvalues (might be the same).*/
     703           0 :           if (idxP>0) {
     704           0 :             d_prev = PetscAbsReal(sPres->value - PetscRealPart(pep->eigr[0]))/(idxP+0.3);
     705             :           } else {
     706           0 :             d_prev = PetscAbsReal(sPres->value - PetscRealPart(pep->eigr[pep->nconv-1]))/(pep->nconv+0.3);
     707             :           }
     708           0 :           *newS = sPres->value + (sr->dir*d_prev*pep->nev)/2;
     709           0 :           sr->dirch = PETSC_FALSE;
     710             :         } else { /* No values found, no information for next shift */
     711           1 :           PetscCheck(!sr->dirch,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"First shift renders no information");
     712           1 :           sr->dirch = PETSC_TRUE;
     713           1 :           *newS = sr->int1;
     714             :         }
     715             :       }
     716             :     } else { /* Accepted values found */
     717           7 :       sr->dirch = PETSC_FALSE;
     718           7 :       sr->nleap = 0;
     719             :       /* Average distance of values in previous subinterval */
     720           7 :       s = sPres->neighb[0];
     721           7 :       while (s && PetscAbs(s->inertia - sPres->inertia)==0) {
     722           0 :         s = s->neighb[0];/* Looking for previous shifts with eigenvalues within */
     723             :       }
     724           7 :       if (s) {
     725           3 :         d_prev = PetscAbsReal((sPres->value - s->value)/(sPres->inertia - s->inertia));
     726             :       } else { /* First shift. Average distance obtained with values in this shift */
     727             :         /* first shift might be too far from first wanted eigenvalue (no values found outside the interval)*/
     728           4 :         if (sr->dir*(PetscRealPart(sr->eigr[0])-sPres->value)>0 && PetscAbsReal((PetscRealPart(sr->eigr[sr->indexEig-1]) - PetscRealPart(sr->eigr[0]))/PetscRealPart(sr->eigr[0])) > PetscSqrtReal(pep->tol)) {
     729           3 :           d_prev =  PetscAbsReal((PetscRealPart(sr->eigr[sr->indexEig-1]) - PetscRealPart(sr->eigr[0])))/(sPres->neigs+0.3);
     730             :         } else {
     731           1 :           d_prev = PetscAbsReal(PetscRealPart(sr->eigr[sr->indexEig-1]) - sPres->value)/(sPres->neigs+0.3);
     732             :         }
     733             :       }
     734             :       /* Average distance is used for next shift by adding it to value on the right or to shift */
     735           7 :       if (sr->dir*(PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1]) - sPres->value)>0) {
     736           7 :         *newS = PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1])+ (sr->dir*d_prev*pep->nev)/2;
     737             :       } else { /* Last accepted value is on the left of shift. Adding to shift */
     738           0 :         *newS = sPres->value + (sr->dir*d_prev*pep->nev)/2;
     739             :       }
     740             :     }
     741             :     /* End of interval can not be surpassed */
     742          12 :     if (sr->dir*(sr->int1 - *newS) < 0) *newS = sr->int1;
     743             :   }/* of neighb[side]==null */
     744          15 :   PetscFunctionReturn(PETSC_SUCCESS);
     745             : }
     746             : 
     747             : /*
     748             :   Function for sorting an array of real values
     749             : */
     750          43 : static PetscErrorCode sortRealEigenvalues(PetscScalar *r,PetscInt *perm,PetscInt nr,PetscBool prev,PetscInt dir)
     751             : {
     752          43 :   PetscReal re;
     753          43 :   PetscInt  i,j,tmp;
     754             : 
     755          43 :   PetscFunctionBegin;
     756        1006 :   if (!prev) for (i=0;i<nr;i++) perm[i] = i;
     757             :   /* Insertion sort */
     758         964 :   for (i=1;i<nr;i++) {
     759         921 :     re = PetscRealPart(r[perm[i]]);
     760         921 :     j = i-1;
     761        3498 :     while (j>=0 && dir*(re - PetscRealPart(r[perm[j]])) <= 0) {
     762        2577 :       tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
     763             :     }
     764             :   }
     765          43 :   PetscFunctionReturn(PETSC_SUCCESS);
     766             : }
     767             : 
     768             : /* Stores the pairs obtained since the last shift in the global arrays */
     769          22 : static PetscErrorCode PEPStoreEigenpairs(PEP pep)
     770             : {
     771          22 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
     772          22 :   PetscReal      lambda,err,*errest;
     773          22 :   PetscInt       i,*aux,count=0,ndef,ld,nconv=pep->nconv,d=pep->nmat-1,idx;
     774          22 :   PetscBool      iscayley,divide=PETSC_FALSE;
     775          22 :   PEP_SR         sr = ctx->sr;
     776          22 :   PEP_shift      sPres;
     777          22 :   Vec            w,vomega;
     778          22 :   Mat            MS;
     779          22 :   BV             tV;
     780          22 :   PetscScalar    *S,*eigr,*tS,*omega;
     781             : 
     782          22 :   PetscFunctionBegin;
     783          22 :   sPres = sr->sPres;
     784          22 :   sPres->index = sr->indexEig;
     785             : 
     786          22 :   if (nconv>sr->ndef0+sr->ndef1) {
     787             :     /* Back-transform */
     788          21 :     PetscCall(STBackTransform(pep->st,nconv,pep->eigr,pep->eigi));
     789         450 :     for (i=0;i<nconv;i++) {
     790             : #if defined(PETSC_USE_COMPLEX)
     791         429 :       if (PetscImaginaryPart(pep->eigr[i])) pep->eigr[i] = sr->int0-sr->dir;
     792             : #else
     793             :       if (pep->eigi[i]) pep->eigr[i] = sr->int0-sr->dir;
     794             : #endif
     795             :     }
     796          21 :     PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STCAYLEY,&iscayley));
     797             :     /* Sort eigenvalues */
     798          21 :     PetscCall(sortRealEigenvalues(pep->eigr,pep->perm,nconv,PETSC_FALSE,sr->dir));
     799          21 :     PetscCall(VecCreateSeq(PETSC_COMM_SELF,nconv,&vomega));
     800          21 :     PetscCall(BVGetSignature(ctx->V,vomega));
     801          21 :     PetscCall(VecGetArray(vomega,&omega));
     802          21 :     PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
     803          21 :     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
     804          21 :     PetscCall(MatDenseGetArray(MS,&S));
     805             :     /* Values stored in global array */
     806          21 :     PetscCall(PetscCalloc4(nconv,&eigr,nconv,&errest,nconv*nconv*d,&tS,nconv,&aux));
     807          21 :     ndef = sr->ndef0+sr->ndef1;
     808         450 :     for (i=0;i<nconv;i++) {
     809         429 :       lambda = PetscRealPart(pep->eigr[pep->perm[i]]);
     810         429 :       err = pep->errest[pep->perm[i]];
     811         429 :       if (sr->dir*(lambda - sPres->ext[0]) > 0 && (sr->dir)*(sPres->ext[1] - lambda) > 0) {/* Valid value */
     812         275 :         PetscCheck(sr->indexEig+count-ndef<sr->numEigs,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Unexpected error in Spectrum Slicing");
     813         275 :         PetscCall(PEPQSliceCheckEigenvalueType(pep,lambda,PetscRealPart(omega[pep->perm[i]]),PETSC_FALSE));
     814         275 :         eigr[count] = lambda;
     815         275 :         errest[count] = err;
     816         275 :         if ((sr->dir*(sPres->value - lambda) > 0) && (sr->dir*(lambda - sPres->ext[0]) > 0)) sPres->nconv[0]++;
     817         275 :         if ((sr->dir*(lambda - sPres->value) > 0) && (sr->dir*(sPres->ext[1] - lambda) > 0)) sPres->nconv[1]++;
     818         275 :         PetscCall(PetscArraycpy(tS+count*(d*nconv),S+pep->perm[i]*(d*ld),nconv));
     819         275 :         PetscCall(PetscArraycpy(tS+count*(d*nconv)+nconv,S+pep->perm[i]*(d*ld)+ld,nconv));
     820         275 :         count++;
     821             :       }
     822             :     }
     823          21 :     PetscCall(VecRestoreArray(vomega,&omega));
     824          21 :     PetscCall(VecDestroy(&vomega));
     825         296 :     for (i=0;i<count;i++) {
     826         275 :       PetscCall(PetscArraycpy(S+i*(d*ld),tS+i*nconv*d,nconv));
     827         275 :       PetscCall(PetscArraycpy(S+i*(d*ld)+ld,tS+i*nconv*d+nconv,nconv));
     828             :     }
     829          21 :     PetscCall(MatDenseRestoreArray(MS,&S));
     830          21 :     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
     831          21 :     PetscCall(BVSetActiveColumns(ctx->V,0,count));
     832          21 :     PetscCall(BVTensorCompress(ctx->V,count));
     833          21 :     if (sr->sPres->nconv[0] && sr->sPres->nconv[1]) {
     834           6 :       divide = PETSC_TRUE;
     835           6 :       PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
     836           6 :       PetscCall(MatDenseGetArray(MS,&S));
     837           6 :       PetscCall(PetscArrayzero(tS,nconv*nconv*d));
     838         137 :       for (i=0;i<count;i++) {
     839         131 :         PetscCall(PetscArraycpy(tS+i*nconv*d,S+i*(d*ld),count));
     840         131 :         PetscCall(PetscArraycpy(tS+i*nconv*d+nconv,S+i*(d*ld)+ld,count));
     841             :       }
     842           6 :       PetscCall(MatDenseRestoreArray(MS,&S));
     843           6 :       PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
     844           6 :       PetscCall(BVSetActiveColumns(pep->V,0,count));
     845           6 :       PetscCall(BVDuplicateResize(pep->V,count,&tV));
     846           6 :       PetscCall(BVCopy(pep->V,tV));
     847             :     }
     848          21 :     if (sr->sPres->nconv[0]) {
     849          11 :       if (divide) {
     850           6 :         PetscCall(BVSetActiveColumns(ctx->V,0,sr->sPres->nconv[0]));
     851           6 :         PetscCall(BVTensorCompress(ctx->V,sr->sPres->nconv[0]));
     852             :       }
     853          94 :       for (i=0;i<sr->ndef0;i++) aux[i] = sr->idxDef0[i];
     854          61 :       for (i=sr->ndef0;i<sr->sPres->nconv[0];i++) aux[i] = sr->indexEig+i-sr->ndef0;
     855          11 :       PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
     856          11 :       PetscCall(MatDenseGetArray(MS,&S));
     857         144 :       for (i=0;i<sr->sPres->nconv[0];i++) {
     858         133 :         sr->eigr[aux[i]] = eigr[i];
     859         133 :         sr->errest[aux[i]] = errest[i];
     860         133 :         PetscCall(BVGetColumn(pep->V,i,&w));
     861         133 :         PetscCall(BVInsertVec(sr->V,aux[i],w));
     862         133 :         PetscCall(BVRestoreColumn(pep->V,i,&w));
     863         133 :         idx = sr->ld*d*aux[i];
     864         133 :         PetscCall(PetscArrayzero(sr->S+idx,sr->ld*d));
     865         133 :         PetscCall(PetscArraycpy(sr->S+idx,S+i*(ld*d),sr->sPres->nconv[0]));
     866         133 :         PetscCall(PetscArraycpy(sr->S+idx+sr->ld,S+i*(ld*d)+ld,sr->sPres->nconv[0]));
     867         133 :         PetscCall(PetscFree(sr->qinfo[aux[i]].q));
     868         133 :         PetscCall(PetscMalloc1(sr->sPres->nconv[0],&sr->qinfo[aux[i]].q));
     869         133 :         PetscCall(PetscArraycpy(sr->qinfo[aux[i]].q,aux,sr->sPres->nconv[0]));
     870         133 :         sr->qinfo[aux[i]].nq = sr->sPres->nconv[0];
     871             :       }
     872          11 :       PetscCall(MatDenseRestoreArray(MS,&S));
     873          11 :       PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
     874             :     }
     875             : 
     876          21 :     if (sr->sPres->nconv[1]) {
     877          13 :       if (divide) {
     878           6 :         PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
     879           6 :         PetscCall(MatDenseGetArray(MS,&S));
     880          73 :         for (i=0;i<sr->sPres->nconv[1];i++) {
     881          67 :           PetscCall(PetscArraycpy(S+i*(d*ld),tS+(sr->sPres->nconv[0]+i)*nconv*d,count));
     882          67 :           PetscCall(PetscArraycpy(S+i*(d*ld)+ld,tS+(sr->sPres->nconv[0]+i)*nconv*d+nconv,count));
     883             :         }
     884           6 :         PetscCall(MatDenseRestoreArray(MS,&S));
     885           6 :         PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
     886           6 :         PetscCall(BVSetActiveColumns(pep->V,0,count));
     887           6 :         PetscCall(BVCopy(tV,pep->V));
     888           6 :         PetscCall(BVSetActiveColumns(ctx->V,0,sr->sPres->nconv[1]));
     889           6 :         PetscCall(BVTensorCompress(ctx->V,sr->sPres->nconv[1]));
     890             :       }
     891          36 :       for (i=0;i<sr->ndef1;i++) aux[i] = sr->idxDef1[i];
     892         132 :       for (i=sr->ndef1;i<sr->sPres->nconv[1];i++) aux[i] = sr->indexEig+sr->sPres->nconv[0]-sr->ndef0+i-sr->ndef1;
     893          13 :       PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
     894          13 :       PetscCall(MatDenseGetArray(MS,&S));
     895         155 :       for (i=0;i<sr->sPres->nconv[1];i++) {
     896         142 :         sr->eigr[aux[i]] = eigr[sr->sPres->nconv[0]+i];
     897         142 :         sr->errest[aux[i]] = errest[sr->sPres->nconv[0]+i];
     898         142 :         PetscCall(BVGetColumn(pep->V,i,&w));
     899         142 :         PetscCall(BVInsertVec(sr->V,aux[i],w));
     900         142 :         PetscCall(BVRestoreColumn(pep->V,i,&w));
     901         142 :         idx = sr->ld*d*aux[i];
     902         142 :         PetscCall(PetscArrayzero(sr->S+idx,sr->ld*d));
     903         142 :         PetscCall(PetscArraycpy(sr->S+idx,S+i*(ld*d),sr->sPres->nconv[1]));
     904         142 :         PetscCall(PetscArraycpy(sr->S+idx+sr->ld,S+i*(ld*d)+ld,sr->sPres->nconv[1]));
     905         142 :         PetscCall(PetscFree(sr->qinfo[aux[i]].q));
     906         142 :         PetscCall(PetscMalloc1(sr->sPres->nconv[1],&sr->qinfo[aux[i]].q));
     907         142 :         PetscCall(PetscArraycpy(sr->qinfo[aux[i]].q,aux,sr->sPres->nconv[1]));
     908         142 :         sr->qinfo[aux[i]].nq = sr->sPres->nconv[1];
     909             :       }
     910          13 :       PetscCall(MatDenseRestoreArray(MS,&S));
     911          13 :       PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
     912             :     }
     913          21 :     sPres->neigs = count-sr->ndef0-sr->ndef1;
     914          21 :     sr->indexEig += sPres->neigs;
     915          21 :     sPres->nconv[0]-= sr->ndef0;
     916          21 :     sPres->nconv[1]-= sr->ndef1;
     917          21 :     PetscCall(PetscFree4(eigr,errest,tS,aux));
     918             :   } else {
     919           1 :     sPres->neigs = 0;
     920           1 :     sPres->nconv[0]= 0;
     921           1 :     sPres->nconv[1]= 0;
     922             :   }
     923             :   /* Global ordering array updating */
     924          22 :   PetscCall(sortRealEigenvalues(sr->eigr,sr->perm,sr->indexEig,PETSC_FALSE,sr->dir));
     925             :   /* Check for completion */
     926          22 :   sPres->comp[0] = PetscNot(sPres->nconv[0] < sPres->nsch[0]);
     927          22 :   sPres->comp[1] = PetscNot(sPres->nconv[1] < sPres->nsch[1]);
     928          22 :   PetscCheck(sPres->nconv[0]<=sPres->nsch[0] && sPres->nconv[1]<=sPres->nsch[1],PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia");
     929          22 :   if (divide) PetscCall(BVDestroy(&tV));
     930          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     931             : }
     932             : 
     933          22 : static PetscErrorCode PEPLookForDeflation(PEP pep)
     934             : {
     935          22 :   PetscReal val;
     936          22 :   PetscInt  i,count0=0,count1=0;
     937          22 :   PEP_shift sPres;
     938          22 :   PetscInt  ini,fin;
     939          22 :   PEP_SR    sr;
     940          22 :   PEP_STOAR *ctx=(PEP_STOAR*)pep->data;
     941             : 
     942          22 :   PetscFunctionBegin;
     943          22 :   sr = ctx->sr;
     944          22 :   sPres = sr->sPres;
     945             : 
     946          22 :   if (sPres->neighb[0]) ini = (sr->dir)*(sPres->neighb[0]->inertia - sr->inertia0);
     947             :   else ini = 0;
     948          22 :   fin = sr->indexEig;
     949             :   /* Selection of ends for searching new values */
     950          22 :   if (!sPres->neighb[0]) sPres->ext[0] = sr->int0;/* First shift */
     951          14 :   else sPres->ext[0] = sPres->neighb[0]->value;
     952          22 :   if (!sPres->neighb[1]) {
     953          19 :     if (sr->hasEnd) sPres->ext[1] = sr->int1;
     954           0 :     else sPres->ext[1] = (sr->dir > 0)?PETSC_MAX_REAL:PETSC_MIN_REAL;
     955           3 :   } else sPres->ext[1] = sPres->neighb[1]->value;
     956             :   /* Selection of values between right and left ends */
     957         128 :   for (i=ini;i<fin;i++) {
     958         108 :     val=PetscRealPart(sr->eigr[sr->perm[i]]);
     959             :     /* Values to the right of left shift */
     960         108 :     if (sr->dir*(val - sPres->ext[1]) < 0) {
     961         106 :       if (sr->dir*(val - sPres->value) < 0) count0++;
     962          23 :       else count1++;
     963             :     } else break;
     964             :   }
     965             :   /* The number of values on each side are found */
     966          22 :   if (sPres->neighb[0]) {
     967          14 :     sPres->nsch[0] = (sr->dir)*(sPres->inertia - sPres->neighb[0]->inertia)-count0;
     968          14 :     PetscCheck(sPres->nsch[0]>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia");
     969           8 :   } else sPres->nsch[0] = 0;
     970             : 
     971          22 :   if (sPres->neighb[1]) {
     972           3 :     sPres->nsch[1] = (sr->dir)*(sPres->neighb[1]->inertia - sPres->inertia) - count1;
     973           3 :     PetscCheck(sPres->nsch[1]>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia");
     974          19 :   } else sPres->nsch[1] = (sr->dir)*(sr->inertia1 - sPres->inertia);
     975             : 
     976             :   /* Completing vector of indexes for deflation */
     977         105 :   for (i=0;i<count0;i++) sr->idxDef0[i] = sr->perm[ini+i];
     978          22 :   sr->ndef0 = count0;
     979          45 :   for (i=0;i<count1;i++) sr->idxDef1[i] = sr->perm[ini+count0+i];
     980          22 :   sr->ndef1 = count1;
     981          22 :   PetscFunctionReturn(PETSC_SUCCESS);
     982             : }
     983             : 
     984             : /*
     985             :   Compute a run of Lanczos iterations
     986             : */
     987         194 : static PetscErrorCode PEPSTOARrun_QSlice(PEP pep,PetscReal *a,PetscReal *b,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,Vec *t_)
     988             : {
     989         194 :   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
     990         194 :   PetscInt       i,j,m=*M,l,lock;
     991         194 :   PetscInt       lds,d,ld,offq,nqt,ldds;
     992         194 :   Vec            v=t_[0],t=t_[1],q=t_[2];
     993         194 :   PetscReal      norm,sym=0.0,fro=0.0,*f;
     994         194 :   PetscScalar    *y,*S,sigma;
     995         194 :   PetscBLASInt   j_,one=1;
     996         194 :   PetscBool      lindep;
     997         194 :   Mat            MS;
     998             : 
     999         194 :   PetscFunctionBegin;
    1000         194 :   PetscCall(PetscMalloc1(*M,&y));
    1001         194 :   PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
    1002         194 :   PetscCall(BVTensorGetDegree(ctx->V,&d));
    1003         194 :   PetscCall(BVGetActiveColumns(pep->V,&lock,&nqt));
    1004         194 :   lds = d*ld;
    1005         194 :   offq = ld;
    1006         194 :   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
    1007             : 
    1008         194 :   *breakdown = PETSC_FALSE; /* ----- */
    1009         194 :   PetscCall(STGetShift(pep->st,&sigma));
    1010         194 :   PetscCall(DSGetDimensions(pep->ds,NULL,&l,NULL,NULL));
    1011         194 :   PetscCall(BVSetActiveColumns(ctx->V,0,m));
    1012         194 :   PetscCall(BVSetActiveColumns(pep->V,0,nqt));
    1013        3370 :   for (j=k;j<m;j++) {
    1014             :     /* apply operator */
    1015        3178 :     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1016        3178 :     PetscCall(MatDenseGetArray(MS,&S));
    1017        3178 :     PetscCall(BVGetColumn(pep->V,nqt,&t));
    1018        3178 :     PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+j*lds));
    1019        3178 :     PetscCall(MatMult(pep->A[1],v,q));
    1020        3178 :     PetscCall(MatMult(pep->A[2],v,t));
    1021        3178 :     PetscCall(VecAXPY(q,sigma*pep->sfactor,t));
    1022        3178 :     PetscCall(VecScale(q,pep->sfactor));
    1023        3178 :     PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds));
    1024        3178 :     PetscCall(MatMult(pep->A[2],v,t));
    1025        3178 :     PetscCall(VecAXPY(q,pep->sfactor*pep->sfactor,t));
    1026        3178 :     PetscCall(STMatSolve(pep->st,q,t));
    1027        3178 :     PetscCall(VecScale(t,-1.0));
    1028        3178 :     PetscCall(BVRestoreColumn(pep->V,nqt,&t));
    1029             : 
    1030             :     /* orthogonalize */
    1031        3178 :     PetscCall(BVOrthogonalizeColumn(pep->V,nqt,S+(j+1)*lds,&norm,&lindep));
    1032        3178 :     if (!lindep) {
    1033        3178 :       *(S+(j+1)*lds+nqt) = norm;
    1034        3178 :       PetscCall(BVScaleColumn(pep->V,nqt,1.0/norm));
    1035        3178 :       nqt++;
    1036             :     }
    1037       88994 :     for (i=0;i<nqt;i++) *(S+(j+1)*lds+offq+i) = *(S+j*lds+i)+sigma*(*(S+(j+1)*lds+i));
    1038        3178 :     PetscCall(BVSetActiveColumns(pep->V,0,nqt));
    1039        3178 :     PetscCall(MatDenseRestoreArray(MS,&S));
    1040        3178 :     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    1041             : 
    1042             :     /* level-2 orthogonalization */
    1043        3178 :     PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep));
    1044        3178 :     a[j] = PetscRealPart(y[j]);
    1045        3178 :     omega[j+1] = (norm > 0)?1.0:-1.0;
    1046        3178 :     PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
    1047        3178 :     b[j] = PetscAbsReal(norm);
    1048             : 
    1049             :     /* check symmetry */
    1050        3178 :     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&f));
    1051        3178 :     if (j==k) {
    1052        3140 :       for (i=l;i<j-1;i++) y[i] = PetscAbsScalar(y[i])-PetscAbsReal(f[2*ldds+i]);
    1053         194 :       for (i=0;i<l;i++) y[i] = 0.0;
    1054             :     }
    1055        3178 :     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&f));
    1056        3178 :     if (j>0) y[j-1] = PetscAbsScalar(y[j-1])-PetscAbsReal(b[j-1]);
    1057        3178 :     PetscCall(PetscBLASIntCast(j,&j_));
    1058        3178 :     sym = SlepcAbs(BLASnrm2_(&j_,y,&one),sym);
    1059        3178 :     fro = SlepcAbs(fro,SlepcAbs(a[j],b[j]));
    1060        3178 :     if (j>0) fro = SlepcAbs(fro,b[j-1]);
    1061        3178 :     if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*pep->tol)) {
    1062           2 :       *symmlost = PETSC_TRUE;
    1063           2 :       *M=j;
    1064           2 :       break;
    1065             :     }
    1066             :   }
    1067         194 :   PetscCall(BVSetActiveColumns(pep->V,lock,nqt));
    1068         194 :   PetscCall(BVSetActiveColumns(ctx->V,0,*M));
    1069         194 :   PetscCall(PetscFree(y));
    1070         194 :   PetscFunctionReturn(PETSC_SUCCESS);
    1071             : }
    1072             : 
    1073          22 : static PetscErrorCode PEPSTOAR_QSlice(PEP pep,Mat B)
    1074             : {
    1075          22 :   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
    1076          22 :   PetscInt       j,k,l,nv=0,ld,ldds,t,nq=0,idx;
    1077          22 :   PetscInt       nconv=0,deg=pep->nmat-1,count0=0,count1=0;
    1078          22 :   PetscScalar    *om,sigma,*back,*S,*pQ;
    1079          22 :   PetscReal      beta,norm=1.0,*omega,*a,*b,eta,lambda;
    1080          22 :   PetscBool      breakdown,symmlost=PETSC_FALSE,sinv,falselock=PETSC_TRUE;
    1081          22 :   Mat            MS,MQ,D;
    1082          22 :   Vec            v,vomega;
    1083          22 :   PEP_SR         sr;
    1084          22 :   BVOrthogType   otype;
    1085          22 :   BVOrthogBlockType obtype;
    1086             : 
    1087          22 :   PetscFunctionBegin;
    1088             :   /* Resize if needed for deflating vectors  */
    1089          22 :   sr = ctx->sr;
    1090          22 :   sigma = sr->sPres->value;
    1091          22 :   k = sr->ndef0+sr->ndef1;
    1092          22 :   pep->ncv = ctx->ncv+k;
    1093          22 :   pep->nev = ctx->nev+k;
    1094          22 :   PetscCall(PEPAllocateSolution(pep,3));
    1095          22 :   PetscCall(BVDestroy(&ctx->V));
    1096          22 :   PetscCall(BVCreateTensor(pep->V,pep->nmat-1,&ctx->V));
    1097          22 :   PetscCall(BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype));
    1098          22 :   PetscCall(BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
    1099          22 :   PetscCall(DSAllocate(pep->ds,pep->ncv+2));
    1100          22 :   PetscCall(PetscMalloc1(pep->ncv,&back));
    1101          22 :   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
    1102          22 :   PetscCall(BVSetMatrix(ctx->V,B,PETSC_TRUE));
    1103          22 :   PetscCheck(ctx->lock,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"A locking variant is needed for spectrum slicing");
    1104             :   /* undocumented option to use a cheaper locking instead of the true locking */
    1105          22 :   PetscCall(PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL));
    1106          22 :   PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
    1107          22 :   PetscCall(RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor));
    1108          22 :   PetscCall(STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor));
    1109             : 
    1110             :   /* Get the starting Arnoldi vector */
    1111          22 :   PetscCall(BVSetActiveColumns(pep->V,0,1));
    1112          22 :   PetscCall(BVTensorBuildFirstColumn(ctx->V,pep->nini));
    1113          22 :   PetscCall(BVSetActiveColumns(ctx->V,0,1));
    1114          22 :   if (k) {
    1115             :     /* Insert deflated vectors */
    1116          10 :     PetscCall(BVSetActiveColumns(pep->V,0,0));
    1117          10 :     idx = sr->ndef0?sr->idxDef0[0]:sr->idxDef1[0];
    1118         116 :     for (j=0;j<k;j++) {
    1119         106 :       PetscCall(BVGetColumn(pep->V,j,&v));
    1120         106 :       PetscCall(BVCopyVec(sr->V,sr->qinfo[idx].q[j],v));
    1121         106 :       PetscCall(BVRestoreColumn(pep->V,j,&v));
    1122             :     }
    1123             :     /* Update innerproduct matrix */
    1124          10 :     PetscCall(BVSetActiveColumns(ctx->V,0,0));
    1125          10 :     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1126          10 :     PetscCall(BVSetActiveColumns(pep->V,0,k));
    1127          10 :     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    1128             : 
    1129          10 :     PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
    1130          10 :     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1131          10 :     PetscCall(MatDenseGetArray(MS,&S));
    1132          93 :     for (j=0;j<sr->ndef0;j++) {
    1133          83 :       PetscCall(PetscArrayzero(S+j*ld*deg,ld*deg));
    1134          83 :       PetscCall(PetscArraycpy(S+j*ld*deg,sr->S+sr->idxDef0[j]*sr->ld*deg,k));
    1135          83 :       PetscCall(PetscArraycpy(S+j*ld*deg+ld,sr->S+sr->idxDef0[j]*sr->ld*deg+sr->ld,k));
    1136          83 :       pep->eigr[j] = sr->eigr[sr->idxDef0[j]];
    1137          83 :       pep->errest[j] = sr->errest[sr->idxDef0[j]];
    1138             :     }
    1139          33 :     for (j=0;j<sr->ndef1;j++) {
    1140          23 :       PetscCall(PetscArrayzero(S+(j+sr->ndef0)*ld*deg,ld*deg));
    1141          23 :       PetscCall(PetscArraycpy(S+(j+sr->ndef0)*ld*deg,sr->S+sr->idxDef1[j]*sr->ld*deg,k));
    1142          23 :       PetscCall(PetscArraycpy(S+(j+sr->ndef0)*ld*deg+ld,sr->S+sr->idxDef1[j]*sr->ld*deg+sr->ld,k));
    1143          23 :       pep->eigr[j+sr->ndef0] = sr->eigr[sr->idxDef1[j]];
    1144          23 :       pep->errest[j+sr->ndef0] = sr->errest[sr->idxDef1[j]];
    1145             :     }
    1146          10 :     PetscCall(MatDenseRestoreArray(MS,&S));
    1147          10 :     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    1148          10 :     PetscCall(BVSetActiveColumns(ctx->V,0,k+1));
    1149          10 :     PetscCall(VecCreateSeq(PETSC_COMM_SELF,k+1,&vomega));
    1150          10 :     PetscCall(VecGetArray(vomega,&om));
    1151         116 :     for (j=0;j<k;j++) {
    1152         106 :       PetscCall(BVOrthogonalizeColumn(ctx->V,j,NULL,&norm,NULL));
    1153         106 :       PetscCall(BVScaleColumn(ctx->V,j,1/norm));
    1154         125 :       om[j] = (norm>=0.0)?1.0:-1.0;
    1155             :     }
    1156          10 :     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1157          10 :     PetscCall(MatDenseGetArray(MS,&S));
    1158          30 :     for (j=0;j<deg;j++) {
    1159          20 :       PetscCall(BVSetRandomColumn(pep->V,k+j));
    1160          20 :       PetscCall(BVOrthogonalizeColumn(pep->V,k+j,S+k*ld*deg+j*ld,&norm,NULL));
    1161          20 :       PetscCall(BVScaleColumn(pep->V,k+j,1.0/norm));
    1162          20 :       S[k*ld*deg+j*ld+k+j] = norm;
    1163             :     }
    1164          10 :     PetscCall(MatDenseRestoreArray(MS,&S));
    1165          10 :     PetscCall(BVSetActiveColumns(pep->V,0,k+deg));
    1166          10 :     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    1167          10 :     PetscCall(BVOrthogonalizeColumn(ctx->V,k,NULL,&norm,NULL));
    1168          10 :     PetscCall(BVScaleColumn(ctx->V,k,1.0/norm));
    1169          10 :     om[k] = (norm>=0.0)?1.0:-1.0;
    1170          10 :     PetscCall(VecRestoreArray(vomega,&om));
    1171          10 :     PetscCall(BVSetSignature(ctx->V,vomega));
    1172          10 :     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&a));
    1173          10 :     PetscCall(VecGetArray(vomega,&om));
    1174         116 :     for (j=0;j<k;j++) a[j] = PetscRealPart(om[j]/(pep->eigr[j]-sigma));
    1175          10 :     PetscCall(VecRestoreArray(vomega,&om));
    1176          10 :     PetscCall(VecDestroy(&vomega));
    1177          10 :     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&a));
    1178          10 :     PetscCall(DSGetArray(pep->ds,DS_MAT_Q,&pQ));
    1179          10 :     PetscCall(PetscArrayzero(pQ,ldds*k));
    1180         116 :     for (j=0;j<k;j++) pQ[j+j*ldds] = 1.0;
    1181          10 :     PetscCall(DSRestoreArray(pep->ds,DS_MAT_Q,&pQ));
    1182             :   }
    1183          22 :   PetscCall(BVSetActiveColumns(ctx->V,0,k+1));
    1184          22 :   PetscCall(DSSetDimensions(pep->ds,k+1,PETSC_DEFAULT,PETSC_DEFAULT));
    1185          22 :   PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
    1186          22 :   PetscCall(BVGetSignature(ctx->V,vomega));
    1187          22 :   PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
    1188             : 
    1189          22 :   PetscCall(PetscInfo(pep,"Start STOAR: sigma=%g in [%g,%g], for deflation: left=%" PetscInt_FMT " right=%" PetscInt_FMT ", searching: left=%" PetscInt_FMT " right=%" PetscInt_FMT "\n",(double)sr->sPres->value,(double)(sr->sPres->neighb[0]?sr->sPres->neighb[0]->value:sr->int0),(double)(sr->sPres->neighb[1]?sr->sPres->neighb[1]->value:sr->int1),sr->ndef0,sr->ndef1,sr->sPres->nsch[0],sr->sPres->nsch[1]));
    1190             : 
    1191             :   /* Restart loop */
    1192          22 :   l = 0;
    1193          22 :   pep->nconv = k;
    1194          22 :   while (pep->reason == PEP_CONVERGED_ITERATING) {
    1195         194 :     pep->its++;
    1196         194 :     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&a));
    1197         194 :     b = a+ldds;
    1198         194 :     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_D,&omega));
    1199             : 
    1200             :     /* Compute an nv-step Lanczos factorization */
    1201         194 :     nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
    1202         194 :     PetscCall(PEPSTOARrun_QSlice(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work));
    1203         194 :     beta = b[nv-1];
    1204         194 :     if (symmlost && nv==pep->nconv+l) {
    1205           1 :       pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
    1206           1 :       pep->nconv = nconv;
    1207           1 :       PetscCall(PetscInfo(pep,"Symmetry lost in STOAR sigma=%g nconv=%" PetscInt_FMT "\n",(double)sr->sPres->value,nconv));
    1208           1 :       if (falselock || !ctx->lock) {
    1209           1 :         PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
    1210           1 :         PetscCall(BVTensorCompress(ctx->V,0));
    1211             :       }
    1212             :       break;
    1213             :     }
    1214         193 :     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&a));
    1215         193 :     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega));
    1216         193 :     PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
    1217         193 :     if (l==0) PetscCall(DSSetState(pep->ds,DS_STATE_INTERMEDIATE));
    1218         171 :     else PetscCall(DSSetState(pep->ds,DS_STATE_RAW));
    1219             : 
    1220             :     /* Solve projected problem */
    1221         193 :     PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
    1222         193 :     PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
    1223         193 :     PetscCall(DSUpdateExtraRow(pep->ds));
    1224         193 :     PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));
    1225             : 
    1226             :     /* Check convergence */
    1227             :     /* PetscCall(PEPSTOARpreKConvergence(pep,nv,&norm,pep->work));*/
    1228         193 :     norm = 1.0;
    1229         193 :     PetscCall(DSGetDimensions(pep->ds,NULL,NULL,NULL,&t));
    1230         193 :     PetscCall(PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k));
    1231         193 :     PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx));
    1232        1238 :     for (j=0;j<k;j++) back[j] = pep->eigr[j];
    1233         193 :     PetscCall(STBackTransform(pep->st,k,back,pep->eigi));
    1234             :     count0=count1=0;
    1235        1238 :     for (j=0;j<k;j++) {
    1236        1045 :       lambda = PetscRealPart(back[j]);
    1237        1045 :       if ((sr->dir*(sr->sPres->value - lambda) > 0) && (sr->dir*(lambda - sr->sPres->ext[0]) > 0)) count0++;
    1238        1045 :       if ((sr->dir*(lambda - sr->sPres->value) > 0) && (sr->dir*(sr->sPres->ext[1] - lambda) > 0)) count1++;
    1239             :     }
    1240         193 :     if ((count0-sr->ndef0 >= sr->sPres->nsch[0]) && (count1-sr->ndef1 >= sr->sPres->nsch[1])) pep->reason = PEP_CONVERGED_TOL;
    1241             :     /* Update l */
    1242         193 :     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
    1243             :     else {
    1244         172 :       l = PetscMax(1,(PetscInt)((nv-k)/2));
    1245         172 :       l = PetscMin(l,t);
    1246         172 :       PetscCall(DSGetTruncateSize(pep->ds,k,t,&l));
    1247         172 :       if (!breakdown) {
    1248             :         /* Prepare the Rayleigh quotient for restart */
    1249         172 :         PetscCall(DSTruncate(pep->ds,k+l,PETSC_FALSE));
    1250             :       }
    1251             :     }
    1252         193 :     nconv = k;
    1253         193 :     if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
    1254         193 :     if (l) PetscCall(PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
    1255             : 
    1256             :     /* Update S */
    1257         193 :     PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&MQ));
    1258         193 :     PetscCall(BVMultInPlace(ctx->V,MQ,pep->nconv,k+l));
    1259         193 :     PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&MQ));
    1260             : 
    1261             :     /* Copy last column of S */
    1262         193 :     PetscCall(BVCopyColumn(ctx->V,nv,k+l));
    1263         193 :     PetscCall(BVSetActiveColumns(ctx->V,0,k+l));
    1264         193 :     if (k+l) {
    1265         192 :       PetscCall(DSSetDimensions(pep->ds,k+l,PETSC_DEFAULT,PETSC_DEFAULT));
    1266         192 :       PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
    1267         192 :       PetscCall(BVSetSignature(ctx->V,vomega));
    1268         192 :       PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
    1269             :     }
    1270             : 
    1271         193 :     if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
    1272             :       /* stop if breakdown */
    1273           0 :       PetscCall(PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta));
    1274           0 :       pep->reason = PEP_DIVERGED_BREAKDOWN;
    1275             :     }
    1276         193 :     if (pep->reason != PEP_CONVERGED_ITERATING) l--;
    1277         193 :     PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
    1278         193 :     if (k+l+deg<=nq) {
    1279         193 :       PetscCall(BVSetActiveColumns(ctx->V,pep->nconv,k+l+1));
    1280         193 :       if (!falselock && ctx->lock) PetscCall(BVTensorCompress(ctx->V,k-pep->nconv));
    1281         193 :       else PetscCall(BVTensorCompress(ctx->V,0));
    1282             :     }
    1283         193 :     pep->nconv = k;
    1284         215 :     PetscCall(PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv));
    1285             :   }
    1286          22 :   sr->itsKs += pep->its;
    1287          22 :   if (pep->nconv>0) {
    1288          21 :     PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
    1289          21 :     PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
    1290          21 :     PetscCall(BVSetActiveColumns(pep->V,0,nq));
    1291          21 :     if (nq>pep->nconv) {
    1292          21 :       PetscCall(BVTensorCompress(ctx->V,pep->nconv));
    1293          21 :       PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
    1294             :     }
    1295         450 :     for (j=0;j<pep->nconv;j++) {
    1296         429 :       pep->eigr[j] *= pep->sfactor;
    1297         429 :       pep->eigi[j] *= pep->sfactor;
    1298             :     }
    1299             :   }
    1300          22 :   PetscCall(PetscInfo(pep,"Finished STOAR: nconv=%" PetscInt_FMT " (deflated=%" PetscInt_FMT ", left=%" PetscInt_FMT ", right=%" PetscInt_FMT ")\n",pep->nconv,sr->ndef0+sr->ndef1,count0-sr->ndef0,count1-sr->ndef1));
    1301          22 :   PetscCall(STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor));
    1302          22 :   PetscCall(RGPopScale(pep->rg));
    1303             : 
    1304          22 :   PetscCheck(pep->reason!=PEP_DIVERGED_SYMMETRY_LOST || nconv>=sr->ndef0+sr->ndef1,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Symmetry lost at sigma=%g",(double)sr->sPres->value);
    1305          22 :   if (pep->reason == PEP_DIVERGED_SYMMETRY_LOST && nconv==sr->ndef0+sr->ndef1) {
    1306           0 :     PetscCheck(++sr->symmlost<=10,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Symmetry lost at sigma=%g",(double)sr->sPres->value);
    1307          22 :   } else sr->symmlost = 0;
    1308             : 
    1309          22 :   PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
    1310          22 :   PetscCall(PetscFree(back));
    1311          22 :   PetscFunctionReturn(PETSC_SUCCESS);
    1312             : }
    1313             : 
    1314             : #define SWAP(a,b,t) do {t=a;a=b;b=t;} while (0)
    1315             : 
    1316           7 : static PetscErrorCode PEPQSliceGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
    1317             : {
    1318           7 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
    1319           7 :   PEP_SR          sr=ctx->sr;
    1320           7 :   PetscInt        i=0,j,tmpi;
    1321           7 :   PetscReal       v,tmpr;
    1322           7 :   PEP_shift       s;
    1323             : 
    1324           7 :   PetscFunctionBegin;
    1325           7 :   PetscCheck(pep->state,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Must call PEPSetUp() first");
    1326           7 :   PetscCheck(sr,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see PEPSetInterval()");
    1327           7 :   if (!sr->s0) {  /* PEPSolve not called yet */
    1328           0 :     *n = 2;
    1329             :   } else {
    1330           7 :     *n = 1;
    1331           7 :     s = sr->s0;
    1332          28 :     while (s) {
    1333          21 :       (*n)++;
    1334          21 :       s = s->neighb[1];
    1335             :     }
    1336             :   }
    1337           7 :   PetscCall(PetscMalloc1(*n,shifts));
    1338           7 :   PetscCall(PetscMalloc1(*n,inertias));
    1339           7 :   if (!sr->s0) {  /* PEPSolve not called yet */
    1340           0 :     (*shifts)[0]   = sr->int0;
    1341           0 :     (*shifts)[1]   = sr->int1;
    1342           0 :     (*inertias)[0] = sr->inertia0;
    1343           0 :     (*inertias)[1] = sr->inertia1;
    1344             :   } else {
    1345             :     s = sr->s0;
    1346          28 :     while (s) {
    1347          21 :       (*shifts)[i]     = s->value;
    1348          21 :       (*inertias)[i++] = s->inertia;
    1349          21 :       s = s->neighb[1];
    1350             :     }
    1351           7 :     (*shifts)[i]   = sr->int1;
    1352           7 :     (*inertias)[i] = sr->inertia1;
    1353             :   }
    1354             :   /* remove possible duplicate in last position */
    1355           7 :   if ((*shifts)[(*n)-1]==(*shifts)[(*n)-2]) (*n)--;
    1356             :   /* sort result */
    1357          32 :   for (i=0;i<*n;i++) {
    1358          25 :     v = (*shifts)[i];
    1359          81 :     for (j=i+1;j<*n;j++) {
    1360          56 :       if (v > (*shifts)[j]) {
    1361          20 :         SWAP((*shifts)[i],(*shifts)[j],tmpr);
    1362          20 :         SWAP((*inertias)[i],(*inertias)[j],tmpi);
    1363          20 :         v = (*shifts)[i];
    1364             :       }
    1365             :     }
    1366             :   }
    1367           7 :   PetscFunctionReturn(PETSC_SUCCESS);
    1368             : }
    1369             : 
    1370           8 : PetscErrorCode PEPSolve_STOAR_QSlice(PEP pep)
    1371             : {
    1372           8 :   PetscInt       i,j,ti,deg=pep->nmat-1;
    1373           8 :   PetscReal      newS;
    1374           8 :   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
    1375           8 :   PEP_SR         sr=ctx->sr;
    1376           8 :   Mat            S,B;
    1377           8 :   PetscScalar    *pS;
    1378             : 
    1379           8 :   PetscFunctionBegin;
    1380           8 :   PetscCall(PetscCitationsRegister(citation,&cited));
    1381             : 
    1382             :   /* Only with eigenvalues present in the interval ...*/
    1383           8 :   if (sr->numEigs==0) {
    1384           1 :     pep->reason = PEP_CONVERGED_TOL;
    1385           1 :     PetscFunctionReturn(PETSC_SUCCESS);
    1386             :   }
    1387             : 
    1388             :   /* Inner product matrix */
    1389           7 :   PetscCall(PEPSTOARSetUpInnerMatrix(pep,&B));
    1390             : 
    1391             :   /* Array of pending shifts */
    1392           7 :   sr->maxPend = 100; /* Initial size */
    1393           7 :   sr->nPend = 0;
    1394           7 :   PetscCall(PetscMalloc1(sr->maxPend,&sr->pending));
    1395           7 :   PetscCall(PEPCreateShift(pep,sr->int0,NULL,NULL));
    1396             :   /* extract first shift */
    1397           7 :   sr->sPrev = NULL;
    1398           7 :   sr->sPres = sr->pending[--sr->nPend];
    1399           7 :   sr->sPres->inertia = sr->inertia0;
    1400           7 :   pep->target = sr->sPres->value;
    1401           7 :   sr->s0 = sr->sPres;
    1402           7 :   sr->indexEig = 0;
    1403             : 
    1404         176 :   for (i=0;i<sr->numEigs;i++) {
    1405         169 :     sr->eigr[i]   = 0.0;
    1406         169 :     sr->eigi[i]   = 0.0;
    1407         169 :     sr->errest[i] = 0.0;
    1408         169 :     sr->perm[i]   = i;
    1409             :   }
    1410             :   /* Vectors for deflation */
    1411           7 :   PetscCall(PetscMalloc2(sr->numEigs,&sr->idxDef0,sr->numEigs,&sr->idxDef1));
    1412           7 :   sr->indexEig = 0;
    1413           7 :   while (sr->sPres) {
    1414             :     /* Search for deflation */
    1415          22 :     PetscCall(PEPLookForDeflation(pep));
    1416             :     /* KrylovSchur */
    1417          22 :     PetscCall(PEPSTOAR_QSlice(pep,B));
    1418             : 
    1419          22 :     PetscCall(PEPStoreEigenpairs(pep));
    1420             :     /* Select new shift */
    1421          22 :     if (!sr->sPres->comp[1]) {
    1422          12 :       PetscCall(PEPGetNewShiftValue(pep,1,&newS));
    1423          12 :       PetscCall(PEPCreateShift(pep,newS,sr->sPres,sr->sPres->neighb[1]));
    1424             :     }
    1425          22 :     if (!sr->sPres->comp[0]) {
    1426             :       /* Completing earlier interval */
    1427           3 :       PetscCall(PEPGetNewShiftValue(pep,0,&newS));
    1428           3 :       PetscCall(PEPCreateShift(pep,newS,sr->sPres->neighb[0],sr->sPres));
    1429             :     }
    1430             :     /* Preparing for a new search of values */
    1431          29 :     PetscCall(PEPExtractShift(pep));
    1432             :   }
    1433             : 
    1434             :   /* Updating pep values prior to exit */
    1435           7 :   PetscCall(PetscFree2(sr->idxDef0,sr->idxDef1));
    1436           7 :   PetscCall(PetscFree(sr->pending));
    1437           7 :   pep->nconv  = sr->indexEig;
    1438           7 :   pep->reason = PEP_CONVERGED_TOL;
    1439           7 :   pep->its    = sr->itsKs;
    1440           7 :   pep->nev    = sr->indexEig;
    1441           7 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,pep->nconv,pep->nconv,NULL,&S));
    1442           7 :   PetscCall(MatDenseGetArray(S,&pS));
    1443         176 :   for (i=0;i<pep->nconv;i++) {
    1444        2466 :     for (j=0;j<sr->qinfo[i].nq;j++) pS[i*pep->nconv+sr->qinfo[i].q[j]] = *(sr->S+i*sr->ld*deg+j);
    1445             :   }
    1446           7 :   PetscCall(MatDenseRestoreArray(S,&pS));
    1447           7 :   PetscCall(BVSetActiveColumns(sr->V,0,pep->nconv));
    1448           7 :   PetscCall(BVMultInPlace(sr->V,S,0,pep->nconv));
    1449           7 :   PetscCall(MatDestroy(&S));
    1450           7 :   PetscCall(BVDestroy(&pep->V));
    1451           7 :   pep->V = sr->V;
    1452           7 :   PetscCall(PetscFree4(pep->eigr,pep->eigi,pep->errest,pep->perm));
    1453           7 :   pep->eigr   = sr->eigr;
    1454           7 :   pep->eigi   = sr->eigi;
    1455           7 :   pep->perm   = sr->perm;
    1456           7 :   pep->errest = sr->errest;
    1457           7 :   if (sr->dir<0) {
    1458          48 :     for (i=0;i<pep->nconv/2;i++) {
    1459          42 :       ti = sr->perm[i]; sr->perm[i] = sr->perm[pep->nconv-1-i]; sr->perm[pep->nconv-1-i] = ti;
    1460             :     }
    1461             :   }
    1462           7 :   PetscCall(PetscFree(ctx->inertias));
    1463           7 :   PetscCall(PetscFree(ctx->shifts));
    1464           7 :   PetscCall(MatDestroy(&B));
    1465           7 :   PetscCall(PEPQSliceGetInertias(pep,&ctx->nshifts,&ctx->shifts,&ctx->inertias));
    1466           7 :   PetscFunctionReturn(PETSC_SUCCESS);
    1467             : }

Generated by: LCOV version 1.14