LCOV - code coverage report
Current view: top level - pep/impls/krylov/stoar - qslice.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 902 1007 89.6 %
Date: 2019-11-15 09:58:42 Functions: 21 22 95.5 %

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

Generated by: LCOV version 1.13