Subversion Repositories slepc-dev

/*                      

   SLEPc eigensolver: "krylovschur"

   Method: Krylov-Schur with spectrum slicing for symmetric eigenproblems

   References:

       [1] R.G. Grimes et al., "A shifted block Lanczos algorithm for solving

           sparse symmetric generalized eigenproblems", SIAM J. Matrix Analysis

           and App., 15(1), pp. 228–272, 1994.

       [2] G.W. Stewart, "A Krylov-Schur Algorithm for Large Eigenproblems",

           SIAM J. Matrix Analysis and App., 23(3), pp. 601-614, 2001.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

   SLEPc - Scalable Library for Eigenvalue Problem Computations

   Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain

   This file is part of SLEPc.

   SLEPc is free software: you can redistribute it and/or modify it under  the

   terms of version 3 of the GNU Lesser General Public License as published by

   the Free Software Foundation.

   SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY

   WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS

   FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for

   more details.

   You  should have received a copy of the GNU Lesser General  Public  License

   along with SLEPc. If not, see <http://www.gnu.org/licenses/>.

   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

*/

#include <private/epsimpl.h>                /*I "slepceps.h" I*/

#include <slepcblaslapack.h>

extern PetscErrorCode EPSProjectedKSSym(EPS,PetscInt,PetscInt,PetscReal*,PetscReal*,PetscScalar*,PetscScalar*,PetscReal*,PetscInt*);

/**/

PetscInt   allKs,def,deg,db;

/* Type of data characterizing a shift (place from where an eps is applied) */

typedef struct _n_shift *shift;

struct _n_shift{

  PetscReal     value;

  PetscInt      inertia;

  PetscBool     comp[2]; //

  shift         neighb[2];//

  PetscInt      index;    //

  PetscInt      neigs;    //

  PetscReal     ext[2];   //

  PetscInt      nsch[2];  //

  PetscReal     pert;     //

  PetscBool     deg;  //not used

};

/* Type of data  for storing the state of spectrum slicing*/

struct _n_SR{

  PetscReal       int0,int1; // extremes of the interval

  PetscInt        dir; // determines the order of values in eig (+1 incr, -1 decr)

  PetscBool       hasEnd; // tells whether the interval has an end

  PetscInt        inertia0,inertia1;

  Vec             *V;

  PetscScalar     *eig;

  PetscInt        *perm;// permutation for keeping the eigenvalues in order

  PetscInt        numEigs; // number of eigenvalues in the interval

  PetscInt        indexEig;

  shift           sPres; // present shift 

  shift           *pending;//pending shifts array

  PetscInt        nPend;// number of pending shifts

  PetscInt        maxPend;// size of "pending" array

  Vec             *VDef; // vector for deflation

  PetscInt        *idxDef;//for deflation

  PetscInt        nMAXCompl;

  PetscInt        iterCompl;

  PetscInt        itsKs;

  PetscInt        nShifts;//number of computed shifts

  shift           s0;// initial shift

  PetscReal       tolDeg;

  PetscInt        nDeg;//number of values coinciding with a shift

  PetscInt        defMin; //minimum amount of values for deflation

};

typedef struct _n_SR  *SR;

/*

   Fills the fields of a shift structure

*/

#undef __FUNCT__

#define __FUNCT__ "EPSCreateShift"

static PetscErrorCode EPSCreateShift(EPS eps,PetscReal val, shift neighb0,shift neighb1)

{

  PetscErrorCode   ierr;

  shift            s,*pending2;

  PetscInt         i;

  SR               sr;

  PetscFunctionBegin;

  sr = (SR)(eps->data);

  ierr = PetscMalloc(sizeof(struct _n_shift),&s);CHKERRQ(ierr);

  s->value = val;

  s->neighb[0] = neighb0;

  if(neighb0) neighb0->neighb[1] = s;

  s->neighb[1] = neighb1;

  if(neighb1) neighb1->neighb[0] = s;

  s->comp[0] = PETSC_FALSE;

  s->comp[1] = PETSC_FALSE;

  s->index = -1;

  s->neigs = 0;

  s->deg = PETSC_FALSE;

  s->nsch[0] = s->nsch[1]=0;

  // inserts in the stack of pending shifts

  // if needed, the array is resized

  if(sr->nPend >= sr->maxPend){

    if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD,"resizing pending shifts array\n");CHKERRQ(ierr);}

    sr->maxPend *= 2;

    ierr = PetscMalloc((sr->maxPend)*sizeof(shift),&pending2);CHKERRQ(ierr);

    for(i=0;i < sr->nPend; i++)pending2[i] = sr->pending[i];

    ierr = PetscFree(sr->pending);CHKERRQ(ierr);

    sr->pending = pending2;

  }

  sr->pending[sr->nPend++]=s;

  sr->nShifts++;

  PetscFunctionReturn(0);

}

/* Provides next shift to be computed */

#undef __FUNCT__

#define __FUNCT__ "EPSExtractShift"

static PetscErrorCode EPSExtractShift(EPS eps){

  PetscErrorCode   ierr;

  PetscInt         iner;

  Mat              F;

  PC               pc;

  KSP              ksp;

  SR               sr;

  PetscFunctionBegin;  

  sr = (SR)(eps->data);

  if(sr->nPend > 0){

    sr->sPres = sr->pending[--sr->nPend];

    ierr = STSetShift(eps->OP, sr->sPres->value);CHKERRQ(ierr);

    ierr = STGetKSP(eps->OP, &ksp);CHKERRQ(ierr);

    ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);

    ierr = PCFactorGetMatrix(pc,&F);CHKERRQ(ierr);

    ierr = MatGetInertia(F,&iner,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);

    sr->sPres->inertia = iner;

    eps->target = sr->sPres->value;

    eps->nconv = 0;

    eps->reason = EPS_CONVERGED_ITERATING;

    eps->its =0;

   }else sr->sPres = PETSC_NULL;

   PetscFunctionReturn(0);

}

/*

   Symmetric KrylovSchur adapted to spectrum slicing:

   allows searching an specific amount of eigenvalues in the subintervals left and right.

   returns whether the search has succeed

*/

#undef __FUNCT__

#define __FUNCT__ "EPSKrylovSchur_Slice"

static PetscErrorCode EPSKrylovSchur_Slice(EPS eps)

{

  PetscErrorCode ierr;

  PetscInt       i,conv,k,l,lds,lt,nv,m,*iwork,p,j;

  Vec            u=eps->work[0];

  PetscScalar    *Q;

  PetscReal      *a,*b,*work,beta,*Qreal,rtmp;//

  PetscBool      breakdown;

  PetscInt       count0,count1;//nconv_prev;

  PetscReal      theta,lambda;

  shift          sPres;

  PetscBool      complIterating;/* shows whether iterations are made for completion */

  PetscBool       sch0,sch1;//shows whether values are looked after on each side

  PetscInt        iterCompl=0,n0,n1;

  //PetscReal       res;

  SR             sr;

  PetscFunctionBegin;

  /* spectrum slicing data */

  sr = (SR)eps->data;

  sPres = sr->sPres;

  complIterating =PETSC_FALSE;

  sch1 = sch0 = PETSC_TRUE;

  lds = PetscMin(eps->mpd,eps->ncv);

  ierr = PetscMalloc(lds*lds*sizeof(PetscReal),&work);CHKERRQ(ierr);

  ierr = PetscMalloc(lds*lds*sizeof(PetscScalar),&Q);CHKERRQ(ierr);

  ierr = PetscMalloc(2*lds*sizeof(PetscInt),&iwork);CHKERRQ(ierr);

  lt = PetscMin(eps->nev+eps->mpd,eps->ncv);

  ierr = PetscMalloc(lt*sizeof(PetscReal),&a);CHKERRQ(ierr);

  ierr = PetscMalloc(lt*sizeof(PetscReal),&b);CHKERRQ(ierr);

  count0=0;count1=0; // found on both sides

  /* Get the starting Lanczos vector */

  ierr = EPSGetStartVector(eps,0,eps->V[0],PETSC_NULL);CHKERRQ(ierr);

  l = 0;

  /* Restart loop */

  while (eps->reason == EPS_CONVERGED_ITERATING) {

    eps->its++;

    /* Compute an nv-step Lanczos factorization */

    m = PetscMin(eps->nconv+eps->mpd,eps->ncv);

    ierr = EPSFullLanczos(eps,a+l,b+l,eps->V,eps->nconv+l,&m,u,&breakdown);CHKERRQ(ierr);

    nv = m - eps->nconv;

    beta = b[nv-1];

    /* Solve projected problem and compute residual norm estimates */

    ierr = EPSProjectedKSSym(eps,nv,l,a,b,eps->eigr+eps->nconv,Q,work,iwork);CHKERRQ(ierr);

    /* Check convergence */

    ierr = EPSKrylovConvergence(eps,PETSC_TRUE,eps->nconv,nv,PETSC_NULL,nv,Q,eps->V+eps->nconv,nv,beta,1.0,&k,PETSC_NULL);CHKERRQ(ierr);

    if(allKs ==1){//option for accepting all converging values

      Qreal = (PetscReal*)Q;//

      conv=k=j=0;

      for(i=0;i<nv;i++)if(eps->errest[eps->nconv+i] < eps->tol)conv++;

      for(i=0;i<nv;i++){

        if(eps->errest[eps->nconv+i] < eps->tol){

          iwork[j++]=i;

        }else iwork[conv+k++]=i;

      }

      for(i=0;i<nv;i++) a[i]=PetscRealPart(eps->eigr[eps->nconv+i]);

      for(i=0;i<nv;i++){

        eps->eigr[eps->nconv+i] = a[iwork[i]];

      }

      for( i=0;i<nv;i++){

        p=iwork[i];

        if(p!=i){

          j=i+1;

          while(iwork[j]!=i)j++;

          iwork[j]=p;iwork[i]=i;

          for(k=0;k<nv;k++){

            rtmp=Qreal[k+p*nv];Qreal[k+p*nv]=Qreal[k+i*nv];Qreal[k+i*nv]=rtmp;

            //rtmp=Q[k+p*nv];Q[k+p*nv]=Q[k+i*nv];Q[k+i*nv]=rtmp; 

          }

        }

      }

      k=eps->nconv+conv;

    }

  /*checking proximity to an eigenvalue*/

  if(deg==1){

    for(i=0;i < k; i++){

      theta = PetscRealPart(eps->eigr[i]);

      if(PetscAbsReal(theta*sPres->value*eps->tol*10)>1){

        if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD,"DEGENERATED SHIFT\n");CHKERRQ(ierr);}

        sr->nDeg++;

        sPres->deg = PETSC_TRUE;

      }else break;

    }

  }

  /*

  if(deg == 1 && sr->nDeg > 0){

    eps->reason = EPS_CONVERGED_TOL;

  }else{

  */

    /* Checking values obtained for completing */

    count0=count1=0;

    for(i=0;i<k;i++){      

      theta = PetscRealPart(eps->eigr[i]);

      lambda = sPres->value + 1/theta;

      if( ((sr->dir)*theta < 0) && ((sr->dir)*(lambda - sPres->ext[0]) > 0))count0++;

      if( ((sr->dir)*theta > 0) && ((sr->dir)*(sPres->ext[1] - lambda) > 0))count1++;

    }

    /* checks completion */

    if( (!sch0||count0 >= sPres->nsch[0]) && (!sch1 ||count1 >= sPres->nsch[1]) ) {

      eps->reason = EPS_CONVERGED_TOL;

    }else {

      if(!complIterating && eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;

      if(complIterating){

        if(--iterCompl <= 0) eps->reason = EPS_DIVERGED_ITS;

      }else if (k >= eps->nev) {//eps->reason = EPS_CONVERGED_TOL;

        n0 = sPres->nsch[0]-count0;

        n1 = sPres->nsch[1]-count1;

        if( sr->iterCompl>0 && ( (n0>0&& n0<= sr->nMAXCompl)||(n1>0&&n1<=sr->nMAXCompl) )){

          complIterating = PETSC_TRUE;

          if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD,"iterating for completion nMAXComp=%d\n",sr->nMAXCompl);CHKERRQ(ierr);}

          if(n0 >sr->nMAXCompl)sch0 = PETSC_FALSE;

          if(n1 >sr->nMAXCompl)sch1 = PETSC_FALSE;

          iterCompl = sr->iterCompl;

        }else eps->reason = EPS_CONVERGED_TOL;

      }      

    }

    /* Update l */

    if (eps->reason != EPS_CONVERGED_ITERATING || breakdown) l = 0;

    else l = (eps->nconv+nv-k)/2;

    if (eps->reason == EPS_CONVERGED_ITERATING) {

      if (breakdown) {

        /* Start a new Lanczos factorization */

        ierr = PetscInfo2(eps,"Breakdown in Krylov-Schur method (it=%D norm=%G)\n",eps->its,beta);CHKERRQ(ierr);

        ierr = EPSGetStartVector(eps,k,eps->V[k],&breakdown);CHKERRQ(ierr);

        if (breakdown) {

          eps->reason = EPS_DIVERGED_BREAKDOWN;

          ierr = PetscInfo(eps,"Unable to generate more start vectors\n");CHKERRQ(ierr);

        }

      } else {

        /* Prepare the Rayleigh quotient for restart */

        for (i=0;i<l;i++) {

          a[i] = PetscRealPart(eps->eigr[i+k]);

          b[i] = PetscRealPart(Q[nv-1+(i+k-eps->nconv)*nv]*beta);

        }

      }

    }

    /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */

    ierr = SlepcUpdateVectors(nv,eps->V+eps->nconv,0,k+l-eps->nconv,Q,nv,PETSC_FALSE);CHKERRQ(ierr);

    /* Normalize u and append it to V */

    if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {

      ierr = VecAXPBY(eps->V[k+l],1.0/beta,0.0,u);CHKERRQ(ierr);

    }

    ierr = EPSMonitor(eps,eps->its,k,eps->eigr,eps->eigi,eps->errest,nv+eps->nconv);CHKERRQ(ierr);

    //nconv_prev = eps->nconv;//

    eps->nconv = k;

  }

  /* check for completion */

  sPres->comp[0] = (count0 >= sPres->nsch[0])?PETSC_TRUE:PETSC_FALSE;

  sPres->comp[1] = (count1 >= sPres->nsch[1])?PETSC_TRUE:PETSC_FALSE;

  if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD," found count0=%d(of %d) and count1=%d(of %d)\n",count0,sPres->nsch[0],count1,sPres->nsch[1]);CHKERRQ(ierr);}

  sr->itsKs += eps->its;  

  ierr = PetscFree(Q);CHKERRQ(ierr);

  ierr = PetscFree(a);CHKERRQ(ierr);

  ierr = PetscFree(b);CHKERRQ(ierr);

  ierr = PetscFree(work);CHKERRQ(ierr);

  ierr = PetscFree(iwork);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

/*

  obtains value of subsequent shift

*/

#undef __FUNCT__

#define __FUNCT__ "EPSGetNewShiftValue"

static PetscErrorCode EPSGetNewShiftValue(EPS eps,PetscInt side,PetscReal *newS){

  PetscReal   lambda,d_prev;//pert,d,d_avg;

  PetscInt    i,idxP;

  SR          sr;

  shift       sPres;

  PetscFunctionBegin;  

  sr = (SR)eps->data;

  sPres = sr->sPres;

/*  

  pert = 0;

  if(sPres->neigs >0){

      idxP=0;//number of computed eigenvalues previous to sPres->value

      for(i=sPres->index;i< sPres->index + sPres->neigs;i++){

        lambda = PetscRealPart(sr->eig[i]);

        if((sr->dir)*(lambda - sPres->value) < 0)idxP++;

        else break;

      }

      //middle point between shift and previous/posterior value

      pert = PetscAbs(sr->eig[sPres->index+idxP]- sr->sPres->value)/2;

  }

*/

  if( sPres->neighb[side]){

  /* completing a previous interval */

    *newS=(sPres->value + sPres->neighb[side]->value)/2;

  }else{ //(only for side=1). creating a new interval.

    if(sPres->neigs==0){// no value has been accepted

      if(sPres->neighb[0]){

        // multiplying by 10 the previous distance

        *newS = sPres->value + 10*(sr->dir)*PetscAbsReal(sPres->value - sPres->neighb[0]->value);

      }else {//first shift

        if(eps->nconv != 0){

           //unaccepted values give information for next shift

           idxP=0;//number of values left from shift 

           for(i=0;i<eps->nconv;i++){

             lambda = PetscRealPart(eps->eigr[i]);

             if( (sr->dir)*(lambda - sPres->value) <0)idxP++;

             else break;

           }

           //avoiding substraction of eigenvalues (might be the same).

           if(idxP>0){

             d_prev = PetscAbsReal(sPres->value - PetscRealPart(eps->eigr[0]))/(idxP+0.3);

           }else {

             d_prev = PetscAbsReal(sPres->value - PetscRealPart(eps->eigr[eps->nconv-1]))/(eps->nconv+0.3);

           }

           *newS = sPres->value + ((sr->dir)*d_prev*eps->nev)/2;

        }else{//no values found, no information for next shift

          // changing the end of the interval

        }

      }

    }else{// accepted values found

      //average distance of values in previous subinterval

      shift s = sPres->neighb[0];

      while(s && PetscAbs(s->inertia - sPres->inertia)==0){

        s = s->neighb[0];//looking for previous shifts with eigenvalues within

      }

      if(s){

        d_prev = PetscAbsReal( (sPres->value - s->value)/(sPres->inertia - s->inertia));

      }else{//first shift. average distance obtained with values in this shift        

        d_prev = PetscAbsReal( PetscRealPart(sr->eig[sPres->index+sPres->neigs-1]) - sPres->value)/(sPres->neigs+0.3);

      }

      // average distance is used for next shift by adding it to value on the right or to shift

      if( (sr->dir)*(PetscRealPart(sr->eig[sPres->index + sPres->neigs -1]) - sPres->value) >0){

        *newS = PetscRealPart(sr->eig[sPres->index + sPres->neigs -1])+ ((sr->dir)*d_prev*(eps->nev))/2;  

      }else{//last accepted value is on the left of shift. adding to shift.

        *newS = sPres->value + ((sr->dir)*d_prev*(eps->nev))/2;

      }

    //}

    }

    //end of interval can not be surpassed

    if(sr->hasEnd && ((sr->dir)*(*newS - sr->int1) > 0))*newS=sr->int1;

  }//of neighb[side]==null

  PetscFunctionReturn(0);

}

/*

  function for sorting an array of real values

*/

#undef __FUNCT__

#define __FUNCT__ "sortRealEigenvalues"

static PetscErrorCode sortRealEigenvalues(PetscScalar *r,PetscInt *perm,PetscInt nr,PetscBool prev,PetscInt dir)

{

  PetscReal      re;

  PetscInt       i,j,tmp;

  PetscFunctionBegin;

  if(!prev) for (i=0; i < nr; i++) { perm[i] = i; }

  /* insertion sort */

  for (i=1; i < nr; i++) {

    re = PetscRealPart(r[perm[i]]);

    j = i-1;

    while ( j>=0 && dir*(re - PetscRealPart(r[perm[j]])) <= 0 ) {

      tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;

    }

  }

  PetscFunctionReturn(0);

}

/* Stores the pairs obtained since the last shift in the global arrays */

#undef __FUNCT__

#define __FUNCT__ "EPSStoreEigenpairs"

PetscErrorCode EPSStoreEigenpairs(EPS eps)

{

  PetscErrorCode ierr;

  PetscReal      lambda,error;

  PetscInt       i,count;

  PetscBool      cond;

  SR             sr;

  shift          sPres;

  PetscFunctionBegin;

  sr = (SR)(eps->data);

  sPres = sr->sPres;

  sPres->index = sr->indexEig;

  count = sr->indexEig;

  error=0;

  /* backtransform */

  for(i=0;i < eps->nconv; i++) eps->eigr[i] =  sPres->value + 1.0/(eps->eigr[i]);

  /* sort eigenvalues */

  ierr = sortRealEigenvalues(eps->eigr,eps->perm,eps->nconv,PETSC_FALSE,sr->dir);

  /* values stored in global array */

  // condition for avoiding comparing whith a non-existing end.

  cond = (!sPres->neighb[1] && !sr->hasEnd)?PETSC_TRUE:PETSC_FALSE;

  for( i=0; i < eps->nconv ;i++ ){

    lambda = PetscRealPart(eps->eigr[eps->perm[i]]);

    if(db>1){ierr = EPSComputeRelativeError(eps,eps->perm[i],&error);CHKERRQ(ierr);}

    if( ( ((sr->dir)*(lambda - sPres->ext[0]) > 0) && ( cond || ((sr->dir)*(lambda - sPres->ext[1]) < 0)) ) ){

      if(count>=sr->numEigs){//Error found

        ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of reserved values exceeded  lambda=%.14g\n",lambda);

        break;

      }  

      sr->eig[count] = lambda;

      ierr = VecCopy(eps->V[eps->perm[i]], sr->V[count]);CHKERRQ(ierr);

      if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"i=%d perm[i]=%d lambda=%.14g error=%g indexEig=%d\n",i,eps->perm[i],lambda,error,count);CHKERRQ(ierr);}

      count++;

    }else if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"i=%d perm[i]=%d lambda=%.14g NOT VALID\n",i,eps->perm[i],lambda);CHKERRQ(ierr);}

  }

  sPres->neigs = count - sr->indexEig;

  if(db>=1){PetscPrintf(PETSC_COMM_WORLD," stored between %d and %d\n",sr->indexEig,count);CHKERRQ(ierr);}

  sr->indexEig = count;

  ierr = sortRealEigenvalues(sr->eig,sr->perm,count,PETSC_TRUE,sr->dir);CHKERRQ(ierr);

  PetscFunctionReturn(0);

}

#undef __FUNCT__

#define __FUNCT__ "EPSLookForDeflation"

PetscErrorCode EPSLookForDeflation(EPS eps)

{

  PetscErrorCode  ierr;

  PetscReal       val,lambda,lambda2;

  PetscInt        i,count0=0,count1=0;

  shift           sPres;

  PetscInt        ini,fin,defMin,k,idx0,idx1;

  PetscBool       consec;

  SR              sr;

  PetscFunctionBegin;

  sr = (SR)(eps->data);

  sPres = sr->sPres;

  defMin = sr->defMin;

  if(sPres->neighb[0]) ini = (sr->dir)*(sPres->neighb[0]->inertia - sr->inertia0);

  else ini = 0;

  fin = sr->indexEig;

  // selection of ends for searching new values

  // later modified with version def=0

  if(!sPres->neighb[0]) sPres->ext[0] = sr->int0;//first shift

  else sPres->ext[0] = sPres->neighb[0]->value;

  if(!sPres->neighb[1]) {

    if(sr->hasEnd) sPres->ext[1] = sr->int1;

    else sPres->ext[1] = (sr->dir > 0)?PETSC_MAX_REAL:PETSC_MIN_REAL;

  }else sPres->ext[1] = sPres->neighb[1]->value;

  if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"ext0=%g ext1=%g\n",sPres->ext[0],sPres->ext[1]);CHKERRQ(ierr);}

  //selection of values between right and left ends

  for(i=ini;i<fin;i++){

    val=PetscRealPart(sr->eig[sr->perm[i]]);

    //values to the right of left shift

    if( (sr->dir)*(val - sPres->ext[1]) < 0 ){

      if((sr->dir)*(val - sPres->value) < 0)count0++;

      else count1++;

    }else break;

  }

  // the number of values on each side are found

  if(sPres->neighb[0])

     sPres->nsch[0] = (sr->dir)*(sPres->inertia - sPres->neighb[0]->inertia)-count0;

  else sPres->nsch[0] = 0;

  if(sPres->neighb[1])

    sPres->nsch[1] = (sr->dir)*(sPres->neighb[1]->inertia - sPres->inertia) - count1;

  else sPres->nsch[1] = (sr->dir)*(sr->inertia1 - sPres->inertia);

  //completing vector of indexes for deflation

  if(def==0 && !sPres->neighb[1]){//new interval && no deflation

    if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"def=0 y neig1=null\n");CHKERRQ(ierr);}

    k=0;

    for(i=fin-1;i>ini;i--){

      k++;

      lambda = PetscRealPart(sr->eig[sr->perm[i]]);

      lambda2 = PetscRealPart(sr->eig[sr->perm[i-1]]);

      if( PetscAbsReal((lambda - lambda2)/lambda) > sr->tolDeg){//relative tolerance 

        break;

      }

    }

    // if i!=ini values for i and i-1 more than toldeg apart

    if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"lookDef i=%d ini=%d\n",i,ini);CHKERRQ(ierr);}

    if(i<=ini){

      sPres->ext[0] = sPres->value;

    }else{//middle point

       sPres->ext[0] = (PetscRealPart(sr->eig[sr->perm[i]])+PetscRealPart(sr->eig[sr->perm[i-1]]))/2;        

    }

    idx0=ini+count0-k;

    idx1=ini+count0;

    if(db>1){ierr = PetscPrintf(PETSC_COMM_WORLD,"ext0=%g ext1=%g idx0=%d idx1=%d count0=%d k=%d\n",sPres->ext[0],sPres->ext[1],idx0,idx1,count0,k);CHKERRQ(ierr);}

  }else{  //completing a subinterval or without deflation

    k = PetscMax(0,defMin-count0);

    idx0 = PetscMax(0,ini-k);

    if(def==0 && sPres->nsch[0]==0){//no deflation towards 0

      idx0 = ini + count0;

      sPres->ext[0] = sPres->value;

    }

    k = PetscMax(0,defMin-count1);

    idx1 = PetscMin(sr->indexEig,ini+count0+count1+k);

    if(def==0 && sPres->nsch[1]==0){//no deflation towards 1

      idx1 = ini + count0;

      sPres->ext[1] = sPres->value;

    }

  }

  k=0;

  for(i=idx0;i<idx1;i++)sr->idxDef[k++]=sr->perm[i];

   ///// info

  if(db>=1){ierr = PetscPrintf(PETSC_COMM_WORLD," deflated %d in [0] and %d in [1]",count0,count1);CHKERRQ(ierr);}

    /////  

  // check for consecutives

  consec=PETSC_TRUE;

  for(i=1;i<k;i++)if(sr->idxDef[i]!=sr->idxDef[i-1]+1){consec = PETSC_FALSE; break;}

  // if not consecutives, copied in array

//if(consec){

//  V o which 

//else{

  for(i=0;i<k;i++)sr->VDef[i]=sr->V[sr->idxDef[i]];

  eps->DS = sr->VDef;

//}

  eps->nds = k;

  //////info

  if(db>=1){

    if(consec){ ierr = PetscPrintf(PETSC_COMM_WORLD," (%d consecutive values)\n",k);CHKERRQ(ierr);}

    else{ ierr = PetscPrintf(PETSC_COMM_WORLD," (%d non consecutive values)\n",k);CHKERRQ(ierr);}

  }

  PetscFunctionReturn(0);

}

#undef __FUNCT__

#define __FUNCT__ "EPSSolve_KrylovSchur_Slice"