Actual source code: qslice.c

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      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"

     13:    Method: S-TOAR with spectrum slicing for symmetric quadratic eigenproblems

     15:    Algorithm:

     17:        Symmetric Two-Level Orthogonal Arnoldi.

     19:    References:

     21:        [1] C. Campos and J.E. Roman, "Inertia-based spectrum slicing
     22:            for symmetric quadratic eigenvalue problems", Numer. Linear
     23:            Algebra Appl. 27(4):e2293, 2020.
     24: */

     26: #include <slepc/private/pepimpl.h>
     27: #include "../src/pep/impls/krylov/pepkrylov.h"
     28: #include <slepcblaslapack.h>

     30: static PetscBool  cited = PETSC_FALSE;
     31: static const char citation[] =
     32:   "@Article{slepc-slice-qep,\n"
     33:   "   author = \"C. Campos and J. E. Roman\",\n"
     34:   "   title = \"Inertia-based spectrum slicing for symmetric quadratic eigenvalue problems\",\n"
     35:   "   journal = \"Numer. Linear Algebra Appl.\",\n"
     36:   "   volume = \"27\",\n"
     37:   "   number = \"4\",\n"
     38:   "   pages = \"e2293\",\n"
     39:   "   year = \"2020,\"\n"
     40:   "   doi = \"https://doi.org/10.1002/nla.2293\"\n"
     41:   "}\n";

     43: #define SLICE_PTOL PETSC_SQRT_MACHINE_EPSILON

     45: static PetscErrorCode PEPQSliceResetSR(PEP pep)
     46: {
     47:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
     48:   PEP_SR         sr=ctx->sr;
     49:   PEP_shift      s;
     50:   PetscInt       i;

     52:   PetscFunctionBegin;
     53:   if (sr) {
     54:     /* Reviewing list of shifts to free memory */
     55:     s = sr->s0;
     56:     if (s) {
     57:       while (s->neighb[1]) {
     58:         s = s->neighb[1];
     59:         PetscCall(PetscFree(s->neighb[0]));
     60:       }
     61:       PetscCall(PetscFree(s));
     62:     }
     63:     PetscCall(PetscFree(sr->S));
     64:     for (i=0;i<pep->nconv;i++) PetscCall(PetscFree(sr->qinfo[i].q));
     65:     PetscCall(PetscFree(sr->qinfo));
     66:     for (i=0;i<3;i++) PetscCall(VecDestroy(&sr->v[i]));
     67:     PetscCall(EPSDestroy(&sr->eps));
     68:     PetscCall(PetscFree(sr));
     69:   }
     70:   ctx->sr = NULL;
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: PetscErrorCode PEPReset_STOAR_QSlice(PEP pep)
     75: {
     76:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;

     78:   PetscFunctionBegin;
     79:   PetscCall(PEPQSliceResetSR(pep));
     80:   PetscCall(PetscFree(ctx->inertias));
     81:   PetscCall(PetscFree(ctx->shifts));
     82:   PetscFunctionReturn(PETSC_SUCCESS);
     83: }

     85: /*
     86:   PEPQSliceAllocateSolution - Allocate memory storage for common variables such
     87:   as eigenvalues and eigenvectors.
     88: */
     89: static PetscErrorCode PEPQSliceAllocateSolution(PEP pep)
     90: {
     91:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
     92:   PetscInt       k;
     93:   BVType         type;
     94:   Vec            t;
     95:   PEP_SR         sr = ctx->sr;

     97:   PetscFunctionBegin;
     98:   /* allocate space for eigenvalues and friends */
     99:   k = PetscMax(1,sr->numEigs);
    100:   PetscCall(PetscFree4(sr->eigr,sr->eigi,sr->errest,sr->perm));
    101:   PetscCall(PetscCalloc4(k,&sr->eigr,k,&sr->eigi,k,&sr->errest,k,&sr->perm));
    102:   PetscCall(PetscFree(sr->qinfo));
    103:   PetscCall(PetscCalloc1(k,&sr->qinfo));

    105:   /* allocate sr->V and transfer options from pep->V */
    106:   PetscCall(BVDestroy(&sr->V));
    107:   PetscCall(BVCreate(PetscObjectComm((PetscObject)pep),&sr->V));
    108:   if (!pep->V) PetscCall(PEPGetBV(pep,&pep->V));
    109:   if (!((PetscObject)pep->V)->type_name) PetscCall(BVSetType(sr->V,BVMAT));
    110:   else {
    111:     PetscCall(BVGetType(pep->V,&type));
    112:     PetscCall(BVSetType(sr->V,type));
    113:   }
    114:   PetscCall(STMatCreateVecsEmpty(pep->st,&t,NULL));
    115:   PetscCall(BVSetSizesFromVec(sr->V,t,k+1));
    116:   PetscCall(VecDestroy(&t));
    117:   sr->ld = k;
    118:   PetscCall(PetscFree(sr->S));
    119:   PetscCall(PetscMalloc1((k+1)*sr->ld*(pep->nmat-1),&sr->S));
    120:   PetscFunctionReturn(PETSC_SUCCESS);
    121: }

    123: /* Convergence test to compute positive Ritz values */
    124: static PetscErrorCode ConvergedPositive(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
    125: {
    126:   PetscFunctionBegin;
    127:   *errest = (PetscRealPart(eigr)>0.0)?0.0:res;
    128:   PetscFunctionReturn(PETSC_SUCCESS);
    129: }

    131: static PetscErrorCode PEPQSliceMatGetInertia(PEP pep,PetscReal shift,PetscInt *inertia,PetscInt *zeros)
    132: {
    133:   KSP            ksp,kspr;
    134:   PC             pc;
    135:   Mat            F;
    136:   PetscBool      flg;

    138:   PetscFunctionBegin;
    139:   if (!pep->solvematcoeffs) PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
    140:   if (shift==PETSC_MAX_REAL) { /* Inertia of matrix A[2] */
    141:     pep->solvematcoeffs[0] = 0.0; pep->solvematcoeffs[1] = 0.0; pep->solvematcoeffs[2] = 1.0;
    142:   } else PetscCall(PEPEvaluateBasis(pep,shift,0,pep->solvematcoeffs,NULL));
    143:   PetscCall(STMatSetUp(pep->st,pep->sfactor,pep->solvematcoeffs));
    144:   PetscCall(STGetKSP(pep->st,&ksp));
    145:   PetscCall(KSPGetPC(ksp,&pc));
    146:   PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCREDUNDANT,&flg));
    147:   if (flg) {
    148:     PetscCall(PCRedundantGetKSP(pc,&kspr));
    149:     PetscCall(KSPGetPC(kspr,&pc));
    150:   }
    151:   PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCTELESCOPE,&flg));
    152:   if (flg) {
    153:     PetscCall(PCTelescopeGetKSP(pc,&kspr));
    154:     PetscCall(KSPGetPC(kspr,&pc));
    155:   }
    156:   PetscCall(PCFactorGetMatrix(pc,&F));
    157:   PetscCall(PCSetUp(pc));
    158:   PetscCall(MatGetInertia(F,inertia,zeros,NULL));
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }

    162: static PetscErrorCode PEPQSliceGetInertia(PEP pep,PetscReal shift,PetscInt *inertia,PetscInt *zeros,PetscInt correction)
    163: {
    164:   KSP            ksp;
    165:   Mat            P;
    166:   PetscReal      nzshift=0.0,dot;
    167:   PetscRandom    rand;
    168:   PetscInt       nconv;
    169:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
    170:   PEP_SR         sr=ctx->sr;

    172:   PetscFunctionBegin;
    173:   if (shift >= PETSC_MAX_REAL) { /* Right-open interval */
    174:     *inertia = 0;
    175:   } else if (shift <= PETSC_MIN_REAL) {
    176:     *inertia = 0;
    177:     if (zeros) *zeros = 0;
    178:   } else {
    179:     /* If the shift is zero, perturb it to a very small positive value.
    180:        The goal is that the nonzero pattern is the same in all cases and reuse
    181:        the symbolic factorizations */
    182:     nzshift = (shift==0.0)? 10.0/PETSC_MAX_REAL: shift;
    183:     PetscCall(PEPQSliceMatGetInertia(pep,nzshift,inertia,zeros));
    184:     PetscCall(STSetShift(pep->st,nzshift));
    185:   }
    186:   if (!correction) {
    187:     if (shift >= PETSC_MAX_REAL) *inertia = 2*pep->n;
    188:     else if (shift>PETSC_MIN_REAL) {
    189:       PetscCall(STGetKSP(pep->st,&ksp));
    190:       PetscCall(KSPGetOperators(ksp,&P,NULL));
    191:       if (*inertia!=pep->n && !sr->v[0]) {
    192:         PetscCall(MatCreateVecs(P,&sr->v[0],NULL));
    193:         PetscCall(VecDuplicate(sr->v[0],&sr->v[1]));
    194:         PetscCall(VecDuplicate(sr->v[0],&sr->v[2]));
    195:         PetscCall(BVGetRandomContext(pep->V,&rand));
    196:         PetscCall(VecSetRandom(sr->v[0],rand));
    197:       }
    198:       if (*inertia<pep->n && *inertia>0) {
    199:         if (!sr->eps) {
    200:           PetscCall(EPSCreate(PetscObjectComm((PetscObject)pep),&sr->eps));
    201:           PetscCall(EPSSetProblemType(sr->eps,EPS_HEP));
    202:           PetscCall(EPSSetWhichEigenpairs(sr->eps,EPS_LARGEST_REAL));
    203:         }
    204:         PetscCall(EPSSetConvergenceTestFunction(sr->eps,ConvergedPositive,NULL,NULL));
    205:         PetscCall(EPSSetOperators(sr->eps,P,NULL));
    206:         PetscCall(EPSSolve(sr->eps));
    207:         PetscCall(EPSGetConverged(sr->eps,&nconv));
    208:         PetscCheck(nconv,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)nzshift);
    209:         PetscCall(EPSGetEigenpair(sr->eps,0,NULL,NULL,sr->v[0],sr->v[1]));
    210:       }
    211:       if (*inertia!=pep->n) {
    212:         PetscCall(MatMult(pep->A[1],sr->v[0],sr->v[1]));
    213:         PetscCall(MatMult(pep->A[2],sr->v[0],sr->v[2]));
    214:         PetscCall(VecAXPY(sr->v[1],2*nzshift,sr->v[2]));
    215:         PetscCall(VecDotRealPart(sr->v[1],sr->v[0],&dot));
    216:         if (dot>0.0) *inertia = 2*pep->n-*inertia;
    217:       }
    218:     }
    219:   } else if (correction<0) *inertia = 2*pep->n-*inertia;
    220:   PetscFunctionReturn(PETSC_SUCCESS);
    221: }

    223: /*
    224:    Check eigenvalue type - used only in non-hyperbolic problems.
    225:    All computed eigenvalues must have the same definite type (positive or negative).
    226:    If ini=TRUE the type is available in omega, otherwise we compute an eigenvalue
    227:    closest to shift and determine its type.
    228:  */
    229: static PetscErrorCode PEPQSliceCheckEigenvalueType(PEP pep,PetscReal shift,PetscReal omega,PetscBool ini)
    230: {
    231:   PEP            pep2;
    232:   ST             st;
    233:   PetscInt       nconv;
    234:   PetscScalar    lambda;
    235:   PetscReal      dot;
    236:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
    237:   PEP_SR         sr=ctx->sr;

    239:   PetscFunctionBegin;
    240:   if (!ini) {
    241:     PetscCheck(-(omega/(shift*ctx->alpha+ctx->beta))*sr->type>=0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Different positive/negative type detected in eigenvalue %g",(double)shift);
    242:   } else {
    243:     PetscCall(PEPCreate(PetscObjectComm((PetscObject)pep),&pep2));
    244:     PetscCall(PEPSetOptionsPrefix(pep2,((PetscObject)pep)->prefix));
    245:     PetscCall(PEPAppendOptionsPrefix(pep2,"pep_eigenvalue_type_"));
    246:     PetscCall(PEPSetTolerances(pep2,PETSC_CURRENT,pep->max_it/4));
    247:     PetscCall(PEPSetType(pep2,PEPTOAR));
    248:     PetscCall(PEPSetOperators(pep2,pep->nmat,pep->A));
    249:     PetscCall(PEPSetWhichEigenpairs(pep2,PEP_TARGET_MAGNITUDE));
    250:     PetscCall(PEPGetRG(pep2,&pep2->rg));
    251:     PetscCall(RGSetType(pep2->rg,RGINTERVAL));
    252: #if defined(PETSC_USE_COMPLEX)
    253:     PetscCall(RGIntervalSetEndpoints(pep2->rg,pep->inta,pep->intb,-PETSC_SQRT_MACHINE_EPSILON,PETSC_SQRT_MACHINE_EPSILON));
    254: #else
    255:     PetscCall(RGIntervalSetEndpoints(pep2->rg,pep->inta,pep->intb,0.0,0.0));
    256: #endif
    257:     pep2->target = shift;
    258:     st = pep2->st;
    259:     pep2->st = pep->st;
    260:     PetscCall(PEPSolve(pep2));
    261:     PetscCall(PEPGetConverged(pep2,&nconv));
    262:     if (nconv) {
    263:       PetscCall(PEPGetEigenpair(pep2,0,&lambda,NULL,pep2->work[0],NULL));
    264:       PetscCall(MatMult(pep->A[1],pep2->work[0],pep2->work[1]));
    265:       PetscCall(MatMult(pep->A[2],pep2->work[0],pep2->work[2]));
    266:       PetscCall(VecAXPY(pep2->work[1],2.0*lambda*pep->sfactor,pep2->work[2]));
    267:       PetscCall(VecDotRealPart(pep2->work[1],pep2->work[0],&dot));
    268:       PetscCall(PetscInfo(pep,"lambda=%g, %s type\n",(double)PetscRealPart(lambda),(dot>0.0)?"positive":"negative"));
    269:       if (!sr->type) sr->type = (dot>0.0)?1:-1;
    270:       else PetscCheck(sr->type*dot>=0.0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Different positive/negative type detected in eigenvalue %g",(double)PetscRealPart(lambda));
    271:     }
    272:     pep2->st = st;
    273:     PetscCall(PEPDestroy(&pep2));
    274:   }
    275:   PetscFunctionReturn(PETSC_SUCCESS);
    276: }

    278: static inline PetscErrorCode PEPQSliceDiscriminant(PEP pep,Vec u,Vec w,PetscReal *d,PetscReal *smas,PetscReal *smenos)
    279: {
    280:   PetscReal      ap,bp,cp,dis;

    282:   PetscFunctionBegin;
    283:   PetscCall(MatMult(pep->A[0],u,w));
    284:   PetscCall(VecDotRealPart(w,u,&cp));
    285:   PetscCall(MatMult(pep->A[1],u,w));
    286:   PetscCall(VecDotRealPart(w,u,&bp));
    287:   PetscCall(MatMult(pep->A[2],u,w));
    288:   PetscCall(VecDotRealPart(w,u,&ap));
    289:   dis = bp*bp-4*ap*cp;
    290:   if (dis>=0.0 && smas) {
    291:     if (ap>0) *smas = (-bp+PetscSqrtReal(dis))/(2*ap);
    292:     else if (ap<0) *smas = (-bp-PetscSqrtReal(dis))/(2*ap);
    293:     else {
    294:       if (bp >0) *smas = -cp/bp;
    295:       else *smas = PETSC_MAX_REAL;
    296:     }
    297:   }
    298:   if (dis>=0.0 && smenos) {
    299:     if (ap>0) *smenos = (-bp-PetscSqrtReal(dis))/(2*ap);
    300:     else if (ap<0) *smenos = (-bp+PetscSqrtReal(dis))/(2*ap);
    301:     else {
    302:       if (bp<0) *smenos = -cp/bp;
    303:       else *smenos = PETSC_MAX_REAL;
    304:     }
    305:   }
    306:   if (d) *d = dis;
    307:   PetscFunctionReturn(PETSC_SUCCESS);
    308: }

    310: static inline PetscErrorCode PEPQSliceEvaluateQEP(PEP pep,PetscScalar x,Mat M,MatStructure str)
    311: {
    312:   PetscFunctionBegin;
    313:   PetscCall(MatCopy(pep->A[0],M,SAME_NONZERO_PATTERN));
    314:   PetscCall(MatAXPY(M,x,pep->A[1],str));
    315:   PetscCall(MatAXPY(M,x*x,pep->A[2],str));
    316:   PetscFunctionReturn(PETSC_SUCCESS);
    317: }

    319: /*@
    320:    PEPCheckDefiniteQEP - Determines if a symmetric/Hermitian quadratic eigenvalue problem
    321:    is definite or not.

    323:    Collective

    325:    Input Parameter:
    326: .  pep - the polynomial eigensolver context

    328:    Output Parameters:
    329: +  xi - first computed parameter
    330: .  mu - second computed parameter
    331: .  definite - flag indicating that the problem is definite
    332: -  hyperbolic - flag indicating that the problem is hyperbolic

    334:    Notes:
    335:    This function is intended for quadratic eigenvalue problems, $Q(\lambda)=K+\lambda C+\lambda^2M$,
    336:    with symmetric (or Hermitian) coefficient matrices $K$, $C$, $M$.

    338:    On output, the flag `definite` may have the values -1 (meaning that the QEP is not
    339:    definite), 1 (if the problem is definite), or 0 if the algorithm was not able to
    340:    determine whether the problem is definite or not.

    342:    If `definite`=1, the output flag `hyperbolic` informs in a similar way about whether the
    343:    problem is hyperbolic or not.

    345:    If `definite`=1, the computed values `xi` and `mu` satisfy $Q(\xi)<0$ and $Q(\mu)>0$,
    346:    as obtained via the method proposed by {cite:t}`Nie10`. Furthermore, if
    347:    `hyperbolic`=1 then only `xi` is computed.

    349:    Level: advanced

    351: .seealso: [](ch:pep), `PEPSetProblemType()`
    352: @*/
    353: PetscErrorCode PEPCheckDefiniteQEP(PEP pep,PetscReal *xi,PetscReal *mu,PetscInt *definite,PetscInt *hyperbolic)
    354: {
    355:   PetscRandom    rand;
    356:   Vec            u,w;
    357:   PetscReal      d=0.0,s=0.0,sp,mut=0.0,omg=0.0,omgp;
    358:   PetscInt       k,its=10,hyp=0,check=0,nconv,inertia,n;
    359:   Mat            M=NULL;
    360:   MatStructure   str;
    361:   EPS            eps;
    362:   PetscBool      transform,ptypehyp;

    364:   PetscFunctionBegin;
    365:   PetscCheck(pep->problem_type==PEP_HERMITIAN || pep->problem_type==PEP_HYPERBOLIC,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Only available for Hermitian (or hyperbolic) problems");
    366:   ptypehyp = (pep->problem_type==PEP_HYPERBOLIC)? PETSC_TRUE: PETSC_FALSE;
    367:   if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
    368:   PetscCall(PEPSetDefaultST(pep));
    369:   PetscCall(STSetMatrices(pep->st,pep->nmat,pep->A));
    370:   PetscCall(MatGetSize(pep->A[0],&n,NULL));
    371:   PetscCall(STGetTransform(pep->st,&transform));
    372:   PetscCall(STSetTransform(pep->st,PETSC_FALSE));
    373:   PetscCall(STSetUp(pep->st));
    374:   PetscCall(MatCreateVecs(pep->A[0],&u,&w));
    375:   PetscCall(PEPGetBV(pep,&pep->V));
    376:   PetscCall(BVGetRandomContext(pep->V,&rand));
    377:   PetscCall(VecSetRandom(u,rand));
    378:   PetscCall(VecNormalize(u,NULL));
    379:   PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,&s,NULL));
    380:   if (d<0.0) check = -1;
    381:   if (!check) {
    382:     PetscCall(EPSCreate(PetscObjectComm((PetscObject)pep),&eps));
    383:     PetscCall(EPSSetProblemType(eps,EPS_HEP));
    384:     PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
    385:     PetscCall(EPSSetTolerances(eps,PetscSqrtReal(PETSC_SQRT_MACHINE_EPSILON),PETSC_DECIDE));
    386:     PetscCall(MatDuplicate(pep->A[0],MAT_DO_NOT_COPY_VALUES,&M));
    387:     PetscCall(STGetMatStructure(pep->st,&str));
    388:   }
    389:   for (k=0;k<its&&!check;k++) {
    390:     PetscCall(PEPQSliceEvaluateQEP(pep,s,M,str));
    391:     PetscCall(EPSSetOperators(eps,M,NULL));
    392:     PetscCall(EPSSolve(eps));
    393:     PetscCall(EPSGetConverged(eps,&nconv));
    394:     if (!nconv) break;
    395:     PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
    396:     sp = s;
    397:     PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,&s,&omg));
    398:     if (d<0.0) {check = -1; break;}
    399:     if (PetscAbsReal((s-sp)/s)<100*PETSC_MACHINE_EPSILON) break;
    400:     if (s>sp) hyp = -1;
    401:     mut = 2*s-sp;
    402:     PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
    403:     if (inertia == n) {check = 1; break;}
    404:   }
    405:   for (;k<its&&!check;k++) {
    406:     mut = (s-omg)/2;
    407:     PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
    408:     if (inertia == n) {check = 1; break;}
    409:     if (PetscAbsReal((s-omg)/omg)<100*PETSC_MACHINE_EPSILON) break;
    410:     PetscCall(PEPQSliceEvaluateQEP(pep,omg,M,str));
    411:     PetscCall(EPSSetOperators(eps,M,NULL));
    412:     PetscCall(EPSSolve(eps));
    413:     PetscCall(EPSGetConverged(eps,&nconv));
    414:     if (!nconv) break;
    415:     PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
    416:     omgp = omg;
    417:     PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,NULL,&omg));
    418:     if (d<0.0) {check = -1; break;}
    419:     if (omg<omgp) hyp = -1;
    420:   }
    421:   if (check==1) *xi = mut;
    422:   PetscCheck(hyp!=-1 || !ptypehyp,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Problem does not satisfy hyperbolic test; consider removing the hyperbolicity flag");
    423:   if (check==1 && hyp==0) {
    424:     PetscCall(PEPQSliceMatGetInertia(pep,PETSC_MAX_REAL,&inertia,NULL));
    425:     if (inertia == 0) hyp = 1;
    426:     else hyp = -1;
    427:   }
    428:   if (check==1 && hyp!=1) {
    429:     check = 0;
    430:     PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
    431:     for (;k<its&&!check;k++) {
    432:       PetscCall(PEPQSliceEvaluateQEP(pep,s,M,str));
    433:       PetscCall(EPSSetOperators(eps,M,NULL));
    434:       PetscCall(EPSSolve(eps));
    435:       PetscCall(EPSGetConverged(eps,&nconv));
    436:       if (!nconv) break;
    437:       PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
    438:       sp = s;
    439:       PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,&s,&omg));
    440:       if (d<0.0) {check = -1; break;}
    441:       if (PetscAbsReal((s-sp)/s)<100*PETSC_MACHINE_EPSILON) break;
    442:       mut = 2*s-sp;
    443:       PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
    444:       if (inertia == 0) {check = 1; break;}
    445:     }
    446:     for (;k<its&&!check;k++) {
    447:       mut = (s-omg)/2;
    448:       PetscCall(PEPQSliceMatGetInertia(pep,mut,&inertia,NULL));
    449:       if (inertia == 0) {check = 1; break;}
    450:       if (PetscAbsReal((s-omg)/omg)<100*PETSC_MACHINE_EPSILON) break;
    451:       PetscCall(PEPQSliceEvaluateQEP(pep,omg,M,str));
    452:       PetscCall(EPSSetOperators(eps,M,NULL));
    453:       PetscCall(EPSSolve(eps));
    454:       PetscCall(EPSGetConverged(eps,&nconv));
    455:       if (!nconv) break;
    456:       PetscCall(EPSGetEigenpair(eps,0,NULL,NULL,u,w));
    457:       PetscCall(PEPQSliceDiscriminant(pep,u,w,&d,NULL,&omg));
    458:       if (d<0.0) {check = -1; break;}
    459:     }
    460:   }
    461:   if (check==1) *mu = mut;
    462:   *definite = check;
    463:   *hyperbolic = hyp;
    464:   PetscCall(MatDestroy(&M));
    465:   PetscCall(VecDestroy(&u));
    466:   PetscCall(VecDestroy(&w));
    467:   PetscCall(EPSDestroy(&eps));
    468:   PetscCall(STSetTransform(pep->st,transform));
    469:   PetscFunctionReturn(PETSC_SUCCESS);
    470: }

    472: /*
    473:    Dummy backtransform operation
    474:  */
    475: static PetscErrorCode PEPBackTransform_Skip(PEP pep)
    476: {
    477:   PetscFunctionBegin;
    478:   PetscFunctionReturn(PETSC_SUCCESS);
    479: }

    481: PetscErrorCode PEPSetUp_STOAR_QSlice(PEP pep)
    482: {
    483:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
    484:   PEP_SR         sr;
    485:   PetscInt       ld,i,zeros=0;
    486:   SlepcSC        sc;
    487:   PetscReal      r;

    489:   PetscFunctionBegin;
    490:   PEPCheckSinvertCayley(pep);
    491:   PetscCheck(pep->inta<pep->intb,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues unless you provide a computational interval with PEPSetInterval()");
    492:   PetscCheck(pep->intb<PETSC_MAX_REAL || pep->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"The defined computational interval should have at least one of their sides bounded");
    493:   PEPCheckUnsupportedCondition(pep,PEP_FEATURE_STOPPING,PETSC_TRUE," (with spectrum slicing)");
    494:   if (pep->tol==(PetscReal)PETSC_DETERMINE) {
    495: #if defined(PETSC_USE_REAL_SINGLE)
    496:     pep->tol = SLEPC_DEFAULT_TOL;
    497: #else
    498:     /* use tighter tolerance */
    499:     pep->tol = SLEPC_DEFAULT_TOL*1e-2;
    500: #endif
    501:   }
    502:   if (ctx->nev==1) ctx->nev = PetscMin(20,pep->n);  /* nev not set, use default value */
    503:   PetscCheck(pep->n<=10 || ctx->nev>=10,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"nev cannot be less than 10 in spectrum slicing runs");
    504:   pep->ops->backtransform = PEPBackTransform_Skip;
    505:   if (pep->max_it==PETSC_DETERMINE) pep->max_it = 100;

    507:   /* create spectrum slicing context and initialize it */
    508:   PetscCall(PEPQSliceResetSR(pep));
    509:   PetscCall(PetscNew(&sr));
    510:   ctx->sr   = sr;
    511:   sr->itsKs = 0;
    512:   sr->nleap = 0;
    513:   sr->sPres = NULL;

    515:   PetscCall(PetscFree(pep->solvematcoeffs));
    516:   PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
    517:   if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
    518:   PetscCall(STSetTransform(pep->st,PETSC_FALSE));
    519:   PetscCall(STSetUp(pep->st));

    521:   ctx->hyperbolic = (pep->problem_type==PEP_HYPERBOLIC)? PETSC_TRUE: PETSC_FALSE;

    523:   /* check presence of ends and finding direction */
    524:   if (pep->inta > PETSC_MIN_REAL || pep->intb >= PETSC_MAX_REAL) {
    525:     sr->int0 = pep->inta;
    526:     sr->int1 = pep->intb;
    527:     sr->dir = 1;
    528:     if (pep->intb >= PETSC_MAX_REAL) { /* Right-open interval */
    529:       sr->hasEnd = PETSC_FALSE;
    530:     } else sr->hasEnd = PETSC_TRUE;
    531:   } else {
    532:     sr->int0 = pep->intb;
    533:     sr->int1 = pep->inta;
    534:     sr->dir = -1;
    535:     sr->hasEnd = PetscNot(pep->inta <= PETSC_MIN_REAL);
    536:   }

    538:   /* compute inertia0 */
    539:   PetscCall(PEPQSliceGetInertia(pep,sr->int0,&sr->inertia0,ctx->detect?&zeros:NULL,ctx->hyperbolic?0:1));
    540:   PetscCheck(!zeros || (sr->int0!=pep->inta && sr->int0!=pep->intb),((PetscObject)pep)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in the interval endpoint");
    541:   if (!ctx->hyperbolic && ctx->checket) PetscCall(PEPQSliceCheckEigenvalueType(pep,sr->int0,0.0,PETSC_TRUE));

    543:   /* compute inertia1 */
    544:   PetscCall(PEPQSliceGetInertia(pep,sr->int1,&sr->inertia1,ctx->detect?&zeros:NULL,ctx->hyperbolic?0:1));
    545:   PetscCheck(!zeros,((PetscObject)pep)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in an interval endpoint defined by user");
    546:   if (!ctx->hyperbolic && ctx->checket && sr->hasEnd) {
    547:     PetscCall(PEPQSliceCheckEigenvalueType(pep,sr->int1,0.0,PETSC_TRUE));
    548:     PetscCheck(sr->type || sr->inertia1==sr->inertia0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"No information of eigenvalue type in interval");
    549:     PetscCheck(!sr->type || sr->inertia1!=sr->inertia0,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Different positive/negative type detected");
    550:     if (sr->dir*(sr->inertia1-sr->inertia0)<0) {
    551:       sr->intcorr = -1;
    552:       sr->inertia0 = 2*pep->n-sr->inertia0;
    553:       sr->inertia1 = 2*pep->n-sr->inertia1;
    554:     } else sr->intcorr = 1;
    555:   } else {
    556:     if (sr->inertia0<=pep->n && sr->inertia1<=pep->n) sr->intcorr = 1;
    557:     else if (sr->inertia0>=pep->n && sr->inertia1>=pep->n) sr->intcorr = -1;
    558:   }

    560:   if (sr->hasEnd) {
    561:     sr->dir = -sr->dir; r = sr->int0; sr->int0 = sr->int1; sr->int1 = r;
    562:     i = sr->inertia0; sr->inertia0 = sr->inertia1; sr->inertia1 = i;
    563:   }

    565:   /* number of eigenvalues in interval */
    566:   sr->numEigs = (sr->dir)*(sr->inertia1 - sr->inertia0);
    567:   PetscCall(PetscInfo(pep,"QSlice setup: allocating for %" PetscInt_FMT " eigenvalues in [%g,%g]\n",sr->numEigs,(double)pep->inta,(double)pep->intb));
    568:   if (sr->numEigs) {
    569:     PetscCall(PEPQSliceAllocateSolution(pep));
    570:     PetscCall(PEPSetDimensions_Default(pep,ctx->nev,&ctx->ncv,&ctx->mpd));
    571:     pep->nev = ctx->nev; pep->ncv = ctx->ncv; pep->mpd = ctx->mpd;
    572:     ld   = ctx->ncv+2;
    573:     PetscCall(DSSetType(pep->ds,DSGHIEP));
    574:     PetscCall(DSSetCompact(pep->ds,PETSC_TRUE));
    575:     PetscCall(DSSetExtraRow(pep->ds,PETSC_TRUE));
    576:     PetscCall(DSAllocate(pep->ds,ld));
    577:     PetscCall(DSGetSlepcSC(pep->ds,&sc));
    578:     sc->rg            = NULL;
    579:     sc->comparison    = SlepcCompareLargestMagnitude;
    580:     sc->comparisonctx = NULL;
    581:     sc->map           = NULL;
    582:     sc->mapobj        = NULL;
    583:   } else {pep->ncv = 0; pep->nev = 0; pep->mpd = 0;}
    584:   PetscFunctionReturn(PETSC_SUCCESS);
    585: }

    587: /*
    588:    Fills the fields of a shift structure
    589: */
    590: static PetscErrorCode PEPCreateShift(PEP pep,PetscReal val,PEP_shift neighb0,PEP_shift neighb1)
    591: {
    592:   PEP_shift      s,*pending2;
    593:   PetscInt       i;
    594:   PEP_SR         sr;
    595:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;

    597:   PetscFunctionBegin;
    598:   sr = ctx->sr;
    599:   PetscCall(PetscNew(&s));
    600:   s->value = val;
    601:   s->neighb[0] = neighb0;
    602:   if (neighb0) neighb0->neighb[1] = s;
    603:   s->neighb[1] = neighb1;
    604:   if (neighb1) neighb1->neighb[0] = s;
    605:   s->comp[0] = PETSC_FALSE;
    606:   s->comp[1] = PETSC_FALSE;
    607:   s->index = -1;
    608:   s->neigs = 0;
    609:   s->nconv[0] = s->nconv[1] = 0;
    610:   s->nsch[0] = s->nsch[1]=0;
    611:   /* Inserts in the stack of pending shifts */
    612:   /* If needed, the array is resized */
    613:   if (sr->nPend >= sr->maxPend) {
    614:     sr->maxPend *= 2;
    615:     PetscCall(PetscMalloc1(sr->maxPend,&pending2));
    616:     for (i=0;i<sr->nPend;i++) pending2[i] = sr->pending[i];
    617:     PetscCall(PetscFree(sr->pending));
    618:     sr->pending = pending2;
    619:   }
    620:   sr->pending[sr->nPend++]=s;
    621:   PetscFunctionReturn(PETSC_SUCCESS);
    622: }

    624: /* Provides next shift to be computed */
    625: static PetscErrorCode PEPExtractShift(PEP pep)
    626: {
    627:   PetscInt       iner,zeros=0;
    628:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
    629:   PEP_SR         sr;
    630:   PetscReal      newShift,aux;
    631:   PEP_shift      sPres;

    633:   PetscFunctionBegin;
    634:   sr = ctx->sr;
    635:   if (sr->nPend > 0) {
    636:     if (sr->dirch) {
    637:       aux = sr->int1; sr->int1 = sr->int0; sr->int0 = aux;
    638:       iner = sr->inertia1; sr->inertia1 = sr->inertia0; sr->inertia0 = iner;
    639:       sr->dir *= -1;
    640:       PetscCall(PetscFree(sr->s0->neighb[1]));
    641:       PetscCall(PetscFree(sr->s0));
    642:       sr->nPend--;
    643:       PetscCall(PEPCreateShift(pep,sr->int0,NULL,NULL));
    644:       sr->sPrev = NULL;
    645:       sr->sPres = sr->pending[--sr->nPend];
    646:       pep->target = sr->sPres->value;
    647:       sr->s0 = sr->sPres;
    648:       pep->reason = PEP_CONVERGED_ITERATING;
    649:     } else {
    650:       sr->sPrev = sr->sPres;
    651:       sr->sPres = sr->pending[--sr->nPend];
    652:     }
    653:     sPres = sr->sPres;
    654:     PetscCall(PEPQSliceGetInertia(pep,sPres->value,&iner,ctx->detect?&zeros:NULL,sr->intcorr));
    655:     if (zeros) {
    656:       newShift = sPres->value*(1.0+SLICE_PTOL);
    657:       if (sr->dir*(sPres->neighb[0] && newShift-sPres->neighb[0]->value) < 0) newShift = (sPres->value+sPres->neighb[0]->value)/2;
    658:       else if (sPres->neighb[1] && sr->dir*(sPres->neighb[1]->value-newShift) < 0) newShift = (sPres->value+sPres->neighb[1]->value)/2;
    659:       PetscCall(PEPQSliceGetInertia(pep,newShift,&iner,&zeros,sr->intcorr));
    660:       PetscCheck(!zeros,((PetscObject)pep)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)newShift);
    661:       sPres->value = newShift;
    662:     }
    663:     sr->sPres->inertia = iner;
    664:     pep->target = sr->sPres->value;
    665:     pep->reason = PEP_CONVERGED_ITERATING;
    666:     pep->its = 0;
    667:   } else sr->sPres = NULL;
    668:   PetscFunctionReturn(PETSC_SUCCESS);
    669: }

    671: /*
    672:   Obtains value of subsequent shift
    673: */
    674: static PetscErrorCode PEPGetNewShiftValue(PEP pep,PetscInt side,PetscReal *newS)
    675: {
    676:   PetscReal lambda,d_prev;
    677:   PetscInt  i,idxP;
    678:   PEP_SR    sr;
    679:   PEP_shift sPres,s;
    680:   PEP_STOAR *ctx=(PEP_STOAR*)pep->data;

    682:   PetscFunctionBegin;
    683:   sr = ctx->sr;
    684:   sPres = sr->sPres;
    685:   if (sPres->neighb[side]) {
    686:   /* Completing a previous interval */
    687:     if (!sPres->neighb[side]->neighb[side] && sPres->neighb[side]->nconv[side]==0) { /* One of the ends might be too far from eigenvalues */
    688:       if (side) *newS = (sPres->value + PetscRealPart(sr->eigr[sr->perm[sr->indexEig-1]]))/2;
    689:       else *newS = (sPres->value + PetscRealPart(sr->eigr[sr->perm[0]]))/2;
    690:     } else *newS=(sPres->value + sPres->neighb[side]->value)/2;
    691:   } else { /* (Only for side=1). Creating a new interval. */
    692:     if (sPres->neigs==0) {/* No value has been accepted*/
    693:       if (sPres->neighb[0]) {
    694:         /* Multiplying by 10 the previous distance */
    695:         *newS = sPres->value + 10*sr->dir*PetscAbsReal(sPres->value - sPres->neighb[0]->value);
    696:         sr->nleap++;
    697:         /* Stops when the interval is open and no values are found in the last 5 shifts (there might be infinite eigenvalues) */
    698:         PetscCheck(sr->hasEnd || sr->nleap<=5,PetscObjectComm((PetscObject)pep),PETSC_ERR_CONV_FAILED,"Unable to compute the wanted eigenvalues with open interval");
    699:       } else { /* First shift */
    700:         if (pep->nconv != 0) {
    701:           /* Unaccepted values give information for next shift */
    702:           idxP=0;/* Number of values left from shift */
    703:           for (i=0;i<pep->nconv;i++) {
    704:             lambda = PetscRealPart(pep->eigr[i]);
    705:             if (sr->dir*(lambda - sPres->value) <0) idxP++;
    706:             else break;
    707:           }
    708:           /* Avoiding subtraction of eigenvalues (might be the same).*/
    709:           if (idxP>0) {
    710:             d_prev = PetscAbsReal(sPres->value - PetscRealPart(pep->eigr[0]))/(idxP+0.3);
    711:           } else {
    712:             d_prev = PetscAbsReal(sPres->value - PetscRealPart(pep->eigr[pep->nconv-1]))/(pep->nconv+0.3);
    713:           }
    714:           *newS = sPres->value + (sr->dir*d_prev*pep->nev)/2;
    715:           sr->dirch = PETSC_FALSE;
    716:         } else { /* No values found, no information for next shift */
    717:           PetscCheck(!sr->dirch,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"First shift renders no information");
    718:           sr->dirch = PETSC_TRUE;
    719:           *newS = sr->int1;
    720:         }
    721:       }
    722:     } else { /* Accepted values found */
    723:       sr->dirch = PETSC_FALSE;
    724:       sr->nleap = 0;
    725:       /* Average distance of values in previous subinterval */
    726:       s = sPres->neighb[0];
    727:       while (s && PetscAbs(s->inertia - sPres->inertia)==0) {
    728:         s = s->neighb[0];/* Looking for previous shifts with eigenvalues within */
    729:       }
    730:       if (s) {
    731:         d_prev = PetscAbsReal((sPres->value - s->value)/(sPres->inertia - s->inertia));
    732:       } else { /* First shift. Average distance obtained with values in this shift */
    733:         /* first shift might be too far from first wanted eigenvalue (no values found outside the interval)*/
    734:         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)) {
    735:           d_prev =  PetscAbsReal((PetscRealPart(sr->eigr[sr->indexEig-1]) - PetscRealPart(sr->eigr[0])))/(sPres->neigs+0.3);
    736:         } else {
    737:           d_prev = PetscAbsReal(PetscRealPart(sr->eigr[sr->indexEig-1]) - sPres->value)/(sPres->neigs+0.3);
    738:         }
    739:       }
    740:       /* Average distance is used for next shift by adding it to value on the right or to shift */
    741:       if (sr->dir*(PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1]) - sPres->value)>0) {
    742:         *newS = PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1])+ (sr->dir*d_prev*pep->nev)/2;
    743:       } else { /* Last accepted value is on the left of shift. Adding to shift */
    744:         *newS = sPres->value + (sr->dir*d_prev*pep->nev)/2;
    745:       }
    746:     }
    747:     /* End of interval can not be surpassed */
    748:     if (sr->dir*(sr->int1 - *newS) < 0) *newS = sr->int1;
    749:   }/* of neighb[side]==null */
    750:   PetscFunctionReturn(PETSC_SUCCESS);
    751: }

    753: /*
    754:   Function for sorting an array of real values
    755: */
    756: static PetscErrorCode sortRealEigenvalues(PetscScalar *r,PetscInt *perm,PetscInt nr,PetscBool prev,PetscInt dir)
    757: {
    758:   PetscReal re;
    759:   PetscInt  i,j,tmp;

    761:   PetscFunctionBegin;
    762:   if (!prev) for (i=0;i<nr;i++) perm[i] = i;
    763:   /* Insertion sort */
    764:   for (i=1;i<nr;i++) {
    765:     re = PetscRealPart(r[perm[i]]);
    766:     j = i-1;
    767:     while (j>=0 && dir*(re - PetscRealPart(r[perm[j]])) <= 0) {
    768:       tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
    769:     }
    770:   }
    771:   PetscFunctionReturn(PETSC_SUCCESS);
    772: }

    774: /* Stores the pairs obtained since the last shift in the global arrays */
    775: static PetscErrorCode PEPStoreEigenpairs(PEP pep)
    776: {
    777:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
    778:   PetscReal      lambda,err,*errest;
    779:   PetscInt       i,*aux,count=0,ndef,ld,nconv=pep->nconv,d=pep->nmat-1,idx;
    780:   PetscBool      iscayley,divide=PETSC_FALSE;
    781:   PEP_SR         sr = ctx->sr;
    782:   PEP_shift      sPres;
    783:   Vec            w,vomega;
    784:   Mat            MS;
    785:   BV             tV;
    786:   PetscScalar    *S,*eigr,*tS,*omega;

    788:   PetscFunctionBegin;
    789:   sPres = sr->sPres;
    790:   sPres->index = sr->indexEig;

    792:   if (nconv>sr->ndef0+sr->ndef1) {
    793:     /* Back-transform */
    794:     PetscCall(STBackTransform(pep->st,nconv,pep->eigr,pep->eigi));
    795:     for (i=0;i<nconv;i++) {
    796: #if defined(PETSC_USE_COMPLEX)
    797:       if (PetscImaginaryPart(pep->eigr[i])) pep->eigr[i] = sr->int0-sr->dir;
    798: #else
    799:       if (pep->eigi[i]) pep->eigr[i] = sr->int0-sr->dir;
    800: #endif
    801:     }
    802:     PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STCAYLEY,&iscayley));
    803:     /* Sort eigenvalues */
    804:     PetscCall(sortRealEigenvalues(pep->eigr,pep->perm,nconv,PETSC_FALSE,sr->dir));
    805:     PetscCall(VecCreateSeq(PETSC_COMM_SELF,nconv,&vomega));
    806:     PetscCall(BVGetSignature(ctx->V,vomega));
    807:     PetscCall(VecGetArray(vomega,&omega));
    808:     PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
    809:     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    810:     PetscCall(MatDenseGetArray(MS,&S));
    811:     /* Values stored in global array */
    812:     PetscCall(PetscCalloc4(nconv,&eigr,nconv,&errest,nconv*nconv*d,&tS,nconv,&aux));
    813:     ndef = sr->ndef0+sr->ndef1;
    814:     for (i=0;i<nconv;i++) {
    815:       lambda = PetscRealPart(pep->eigr[pep->perm[i]]);
    816:       err = pep->errest[pep->perm[i]];
    817:       if (sr->dir*(lambda - sPres->ext[0]) > 0 && (sr->dir)*(sPres->ext[1] - lambda) > 0) {/* Valid value */
    818:         PetscCheck(sr->indexEig+count-ndef<sr->numEigs,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Unexpected error in Spectrum Slicing");
    819:         PetscCall(PEPQSliceCheckEigenvalueType(pep,lambda,PetscRealPart(omega[pep->perm[i]]),PETSC_FALSE));
    820:         eigr[count] = lambda;
    821:         errest[count] = err;
    822:         if ((sr->dir*(sPres->value - lambda) > 0) && (sr->dir*(lambda - sPres->ext[0]) > 0)) sPres->nconv[0]++;
    823:         if ((sr->dir*(lambda - sPres->value) > 0) && (sr->dir*(sPres->ext[1] - lambda) > 0)) sPres->nconv[1]++;
    824:         PetscCall(PetscArraycpy(tS+count*(d*nconv),S+pep->perm[i]*(d*ld),nconv));
    825:         PetscCall(PetscArraycpy(tS+count*(d*nconv)+nconv,S+pep->perm[i]*(d*ld)+ld,nconv));
    826:         count++;
    827:       }
    828:     }
    829:     PetscCall(VecRestoreArray(vomega,&omega));
    830:     PetscCall(VecDestroy(&vomega));
    831:     for (i=0;i<count;i++) {
    832:       PetscCall(PetscArraycpy(S+i*(d*ld),tS+i*nconv*d,nconv));
    833:       PetscCall(PetscArraycpy(S+i*(d*ld)+ld,tS+i*nconv*d+nconv,nconv));
    834:     }
    835:     PetscCall(MatDenseRestoreArray(MS,&S));
    836:     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    837:     PetscCall(BVSetActiveColumns(ctx->V,0,count));
    838:     PetscCall(BVTensorCompress(ctx->V,count));
    839:     if (sr->sPres->nconv[0] && sr->sPres->nconv[1]) {
    840:       divide = PETSC_TRUE;
    841:       PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    842:       PetscCall(MatDenseGetArray(MS,&S));
    843:       PetscCall(PetscArrayzero(tS,nconv*nconv*d));
    844:       for (i=0;i<count;i++) {
    845:         PetscCall(PetscArraycpy(tS+i*nconv*d,S+i*(d*ld),count));
    846:         PetscCall(PetscArraycpy(tS+i*nconv*d+nconv,S+i*(d*ld)+ld,count));
    847:       }
    848:       PetscCall(MatDenseRestoreArray(MS,&S));
    849:       PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    850:       PetscCall(BVSetActiveColumns(pep->V,0,count));
    851:       PetscCall(BVDuplicateResize(pep->V,count,&tV));
    852:       PetscCall(BVCopy(pep->V,tV));
    853:     }
    854:     if (sr->sPres->nconv[0]) {
    855:       if (divide) {
    856:         PetscCall(BVSetActiveColumns(ctx->V,0,sr->sPres->nconv[0]));
    857:         PetscCall(BVTensorCompress(ctx->V,sr->sPres->nconv[0]));
    858:       }
    859:       for (i=0;i<sr->ndef0;i++) aux[i] = sr->idxDef0[i];
    860:       for (i=sr->ndef0;i<sr->sPres->nconv[0];i++) aux[i] = sr->indexEig+i-sr->ndef0;
    861:       PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    862:       PetscCall(MatDenseGetArray(MS,&S));
    863:       for (i=0;i<sr->sPres->nconv[0];i++) {
    864:         sr->eigr[aux[i]] = eigr[i];
    865:         sr->errest[aux[i]] = errest[i];
    866:         PetscCall(BVGetColumn(pep->V,i,&w));
    867:         PetscCall(BVInsertVec(sr->V,aux[i],w));
    868:         PetscCall(BVRestoreColumn(pep->V,i,&w));
    869:         idx = sr->ld*d*aux[i];
    870:         PetscCall(PetscArrayzero(sr->S+idx,sr->ld*d));
    871:         PetscCall(PetscArraycpy(sr->S+idx,S+i*(ld*d),sr->sPres->nconv[0]));
    872:         PetscCall(PetscArraycpy(sr->S+idx+sr->ld,S+i*(ld*d)+ld,sr->sPres->nconv[0]));
    873:         PetscCall(PetscFree(sr->qinfo[aux[i]].q));
    874:         PetscCall(PetscMalloc1(sr->sPres->nconv[0],&sr->qinfo[aux[i]].q));
    875:         PetscCall(PetscArraycpy(sr->qinfo[aux[i]].q,aux,sr->sPres->nconv[0]));
    876:         sr->qinfo[aux[i]].nq = sr->sPres->nconv[0];
    877:       }
    878:       PetscCall(MatDenseRestoreArray(MS,&S));
    879:       PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    880:     }

    882:     if (sr->sPres->nconv[1]) {
    883:       if (divide) {
    884:         PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    885:         PetscCall(MatDenseGetArray(MS,&S));
    886:         for (i=0;i<sr->sPres->nconv[1];i++) {
    887:           PetscCall(PetscArraycpy(S+i*(d*ld),tS+(sr->sPres->nconv[0]+i)*nconv*d,count));
    888:           PetscCall(PetscArraycpy(S+i*(d*ld)+ld,tS+(sr->sPres->nconv[0]+i)*nconv*d+nconv,count));
    889:         }
    890:         PetscCall(MatDenseRestoreArray(MS,&S));
    891:         PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    892:         PetscCall(BVSetActiveColumns(pep->V,0,count));
    893:         PetscCall(BVCopy(tV,pep->V));
    894:         PetscCall(BVSetActiveColumns(ctx->V,0,sr->sPres->nconv[1]));
    895:         PetscCall(BVTensorCompress(ctx->V,sr->sPres->nconv[1]));
    896:       }
    897:       for (i=0;i<sr->ndef1;i++) aux[i] = sr->idxDef1[i];
    898:       for (i=sr->ndef1;i<sr->sPres->nconv[1];i++) aux[i] = sr->indexEig+sr->sPres->nconv[0]-sr->ndef0+i-sr->ndef1;
    899:       PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    900:       PetscCall(MatDenseGetArray(MS,&S));
    901:       for (i=0;i<sr->sPres->nconv[1];i++) {
    902:         sr->eigr[aux[i]] = eigr[sr->sPres->nconv[0]+i];
    903:         sr->errest[aux[i]] = errest[sr->sPres->nconv[0]+i];
    904:         PetscCall(BVGetColumn(pep->V,i,&w));
    905:         PetscCall(BVInsertVec(sr->V,aux[i],w));
    906:         PetscCall(BVRestoreColumn(pep->V,i,&w));
    907:         idx = sr->ld*d*aux[i];
    908:         PetscCall(PetscArrayzero(sr->S+idx,sr->ld*d));
    909:         PetscCall(PetscArraycpy(sr->S+idx,S+i*(ld*d),sr->sPres->nconv[1]));
    910:         PetscCall(PetscArraycpy(sr->S+idx+sr->ld,S+i*(ld*d)+ld,sr->sPres->nconv[1]));
    911:         PetscCall(PetscFree(sr->qinfo[aux[i]].q));
    912:         PetscCall(PetscMalloc1(sr->sPres->nconv[1],&sr->qinfo[aux[i]].q));
    913:         PetscCall(PetscArraycpy(sr->qinfo[aux[i]].q,aux,sr->sPres->nconv[1]));
    914:         sr->qinfo[aux[i]].nq = sr->sPres->nconv[1];
    915:       }
    916:       PetscCall(MatDenseRestoreArray(MS,&S));
    917:       PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    918:     }
    919:     sPres->neigs = count-sr->ndef0-sr->ndef1;
    920:     sr->indexEig += sPres->neigs;
    921:     sPres->nconv[0]-= sr->ndef0;
    922:     sPres->nconv[1]-= sr->ndef1;
    923:     PetscCall(PetscFree4(eigr,errest,tS,aux));
    924:   } else {
    925:     sPres->neigs = 0;
    926:     sPres->nconv[0]= 0;
    927:     sPres->nconv[1]= 0;
    928:   }
    929:   /* Global ordering array updating */
    930:   PetscCall(sortRealEigenvalues(sr->eigr,sr->perm,sr->indexEig,PETSC_FALSE,sr->dir));
    931:   /* Check for completion */
    932:   sPres->comp[0] = PetscNot(sPres->nconv[0] < sPres->nsch[0]);
    933:   sPres->comp[1] = PetscNot(sPres->nconv[1] < sPres->nsch[1]);
    934:   PetscCheck(sPres->nconv[0]<=sPres->nsch[0] && sPres->nconv[1]<=sPres->nsch[1],PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia");
    935:   if (divide) PetscCall(BVDestroy(&tV));
    936:   PetscFunctionReturn(PETSC_SUCCESS);
    937: }

    939: static PetscErrorCode PEPLookForDeflation(PEP pep)
    940: {
    941:   PetscReal val;
    942:   PetscInt  i,count0=0,count1=0;
    943:   PEP_shift sPres;
    944:   PetscInt  ini,fin;
    945:   PEP_SR    sr;
    946:   PEP_STOAR *ctx=(PEP_STOAR*)pep->data;

    948:   PetscFunctionBegin;
    949:   sr = ctx->sr;
    950:   sPres = sr->sPres;

    952:   if (sPres->neighb[0]) ini = (sr->dir)*(sPres->neighb[0]->inertia - sr->inertia0);
    953:   else ini = 0;
    954:   fin = sr->indexEig;
    955:   /* Selection of ends for searching new values */
    956:   if (!sPres->neighb[0]) sPres->ext[0] = sr->int0;/* First shift */
    957:   else sPres->ext[0] = sPres->neighb[0]->value;
    958:   if (!sPres->neighb[1]) {
    959:     if (sr->hasEnd) sPres->ext[1] = sr->int1;
    960:     else sPres->ext[1] = (sr->dir > 0)?PETSC_MAX_REAL:PETSC_MIN_REAL;
    961:   } else sPres->ext[1] = sPres->neighb[1]->value;
    962:   /* Selection of values between right and left ends */
    963:   for (i=ini;i<fin;i++) {
    964:     val=PetscRealPart(sr->eigr[sr->perm[i]]);
    965:     /* Values to the right of left shift */
    966:     if (sr->dir*(val - sPres->ext[1]) < 0) {
    967:       if (sr->dir*(val - sPres->value) < 0) count0++;
    968:       else count1++;
    969:     } else break;
    970:   }
    971:   /* The number of values on each side are found */
    972:   if (sPres->neighb[0]) {
    973:     sPres->nsch[0] = (sr->dir)*(sPres->inertia - sPres->neighb[0]->inertia)-count0;
    974:     PetscCheck(sPres->nsch[0]>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia");
    975:   } else sPres->nsch[0] = 0;

    977:   if (sPres->neighb[1]) {
    978:     sPres->nsch[1] = (sr->dir)*(sPres->neighb[1]->inertia - sPres->inertia) - count1;
    979:     PetscCheck(sPres->nsch[1]>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia");
    980:   } else sPres->nsch[1] = (sr->dir)*(sr->inertia1 - sPres->inertia);

    982:   /* Completing vector of indexes for deflation */
    983:   for (i=0;i<count0;i++) sr->idxDef0[i] = sr->perm[ini+i];
    984:   sr->ndef0 = count0;
    985:   for (i=0;i<count1;i++) sr->idxDef1[i] = sr->perm[ini+count0+i];
    986:   sr->ndef1 = count1;
    987:   PetscFunctionReturn(PETSC_SUCCESS);
    988: }

    990: /*
    991:   Compute a run of Lanczos iterations
    992: */
    993: static PetscErrorCode PEPSTOARrun_QSlice(PEP pep,PetscReal *a,PetscReal *b,PetscReal *omega,PetscInt k,PetscInt *M,PetscBool *breakdown,PetscBool *symmlost,Vec *t_)
    994: {
    995:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
    996:   PetscInt       i,j,m=*M,l,lock;
    997:   PetscInt       lds,d,ld,offq,nqt,ldds;
    998:   Vec            v=t_[0],t=t_[1],q=t_[2];
    999:   PetscReal      norm,sym=0.0,fro=0.0,*f;
    1000:   PetscScalar    *y,*S,sigma;
    1001:   PetscBLASInt   j_,one=1;
    1002:   PetscBool      lindep;
    1003:   Mat            MS;

    1005:   PetscFunctionBegin;
    1006:   PetscCall(PetscMalloc1(*M,&y));
    1007:   PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
    1008:   PetscCall(BVTensorGetDegree(ctx->V,&d));
    1009:   PetscCall(BVGetActiveColumns(pep->V,&lock,&nqt));
    1010:   lds = d*ld;
    1011:   offq = ld;
    1012:   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));

    1014:   *breakdown = PETSC_FALSE; /* ----- */
    1015:   PetscCall(STGetShift(pep->st,&sigma));
    1016:   PetscCall(DSGetDimensions(pep->ds,NULL,&l,NULL,NULL));
    1017:   PetscCall(BVSetActiveColumns(ctx->V,0,m));
    1018:   PetscCall(BVSetActiveColumns(pep->V,0,nqt));
    1019:   for (j=k;j<m;j++) {
    1020:     /* apply operator */
    1021:     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1022:     PetscCall(MatDenseGetArray(MS,&S));
    1023:     PetscCall(BVGetColumn(pep->V,nqt,&t));
    1024:     PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+j*lds));
    1025:     PetscCall(MatMult(pep->A[1],v,q));
    1026:     PetscCall(MatMult(pep->A[2],v,t));
    1027:     PetscCall(VecAXPY(q,sigma*pep->sfactor,t));
    1028:     PetscCall(VecScale(q,pep->sfactor));
    1029:     PetscCall(BVMultVec(pep->V,1.0,0.0,v,S+offq+j*lds));
    1030:     PetscCall(MatMult(pep->A[2],v,t));
    1031:     PetscCall(VecAXPY(q,pep->sfactor*pep->sfactor,t));
    1032:     PetscCall(STMatSolve(pep->st,q,t));
    1033:     PetscCall(VecScale(t,-1.0));
    1034:     PetscCall(BVRestoreColumn(pep->V,nqt,&t));

    1036:     /* orthogonalize */
    1037:     PetscCall(BVOrthogonalizeColumn(pep->V,nqt,S+(j+1)*lds,&norm,&lindep));
    1038:     if (!lindep) {
    1039:       *(S+(j+1)*lds+nqt) = norm;
    1040:       PetscCall(BVScaleColumn(pep->V,nqt,1.0/norm));
    1041:       nqt++;
    1042:     }
    1043:     for (i=0;i<nqt;i++) *(S+(j+1)*lds+offq+i) = *(S+j*lds+i)+sigma*(*(S+(j+1)*lds+i));
    1044:     PetscCall(BVSetActiveColumns(pep->V,0,nqt));
    1045:     PetscCall(MatDenseRestoreArray(MS,&S));
    1046:     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));

    1048:     /* level-2 orthogonalization */
    1049:     PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep));
    1050:     a[j] = PetscRealPart(y[j]);
    1051:     omega[j+1] = (norm > 0)?1.0:-1.0;
    1052:     PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
    1053:     b[j] = PetscAbsReal(norm);

    1055:     /* check symmetry */
    1056:     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&f));
    1057:     if (j==k) {
    1058:       for (i=l;i<j-1;i++) y[i] = PetscAbsScalar(y[i])-PetscAbsReal(f[2*ldds+i]);
    1059:       for (i=0;i<l;i++) y[i] = 0.0;
    1060:     }
    1061:     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&f));
    1062:     if (j>0) y[j-1] = PetscAbsScalar(y[j-1])-PetscAbsReal(b[j-1]);
    1063:     PetscCall(PetscBLASIntCast(j,&j_));
    1064:     sym = SlepcAbs(BLASnrm2_(&j_,y,&one),sym);
    1065:     fro = SlepcAbs(fro,SlepcAbs(a[j],b[j]));
    1066:     if (j>0) fro = SlepcAbs(fro,b[j-1]);
    1067:     if (sym/fro>PetscMax(PETSC_SQRT_MACHINE_EPSILON,10*pep->tol)) {
    1068:       *symmlost = PETSC_TRUE;
    1069:       *M=j;
    1070:       break;
    1071:     }
    1072:   }
    1073:   PetscCall(BVSetActiveColumns(pep->V,lock,nqt));
    1074:   PetscCall(BVSetActiveColumns(ctx->V,0,*M));
    1075:   PetscCall(PetscFree(y));
    1076:   PetscFunctionReturn(PETSC_SUCCESS);
    1077: }

    1079: static PetscErrorCode PEPSTOAR_QSlice(PEP pep,Mat B)
    1080: {
    1081:   PEP_STOAR      *ctx = (PEP_STOAR*)pep->data;
    1082:   PetscInt       j,k,l,nv=0,ld,ldds,t,nq=0,idx;
    1083:   PetscInt       nconv=0,deg=pep->nmat-1,count0=0,count1=0;
    1084:   PetscScalar    *om,sigma,*back,*S,*pQ;
    1085:   PetscReal      beta,norm=1.0,*omega,*a,*b,eta,lambda;
    1086:   PetscBool      breakdown,symmlost=PETSC_FALSE,sinv,falselock=PETSC_TRUE;
    1087:   Mat            MS,MQ,D;
    1088:   Vec            v,vomega;
    1089:   PEP_SR         sr;
    1090:   BVOrthogType   otype;
    1091:   BVOrthogBlockType obtype;

    1093:   PetscFunctionBegin;
    1094:   /* Resize if needed for deflating vectors  */
    1095:   sr = ctx->sr;
    1096:   sigma = sr->sPres->value;
    1097:   k = sr->ndef0+sr->ndef1;
    1098:   pep->ncv = ctx->ncv+k;
    1099:   pep->nev = ctx->nev+k;
    1100:   PetscCall(PEPAllocateSolution(pep,3));
    1101:   PetscCall(BVDestroy(&ctx->V));
    1102:   PetscCall(BVCreateTensor(pep->V,pep->nmat-1,&ctx->V));
    1103:   PetscCall(BVGetOrthogonalization(pep->V,&otype,NULL,&eta,&obtype));
    1104:   PetscCall(BVSetOrthogonalization(ctx->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
    1105:   PetscCall(DSAllocate(pep->ds,pep->ncv+2));
    1106:   PetscCall(PetscMalloc1(pep->ncv,&back));
    1107:   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));
    1108:   PetscCall(BVSetMatrix(ctx->V,B,PETSC_TRUE));
    1109:   PetscCheck(ctx->lock,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"A locking variant is needed for spectrum slicing");
    1110:   /* undocumented option to use a cheaper locking instead of the true locking */
    1111:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-pep_stoar_falselocking",&falselock,NULL));
    1112:   PetscCall(PetscObjectTypeCompare((PetscObject)pep->st,STSINVERT,&sinv));
    1113:   PetscCall(RGPushScale(pep->rg,sinv?pep->sfactor:1.0/pep->sfactor));
    1114:   PetscCall(STScaleShift(pep->st,sinv?pep->sfactor:1.0/pep->sfactor));

    1116:   /* Get the starting Arnoldi vector */
    1117:   PetscCall(BVSetActiveColumns(pep->V,0,1));
    1118:   PetscCall(BVTensorBuildFirstColumn(ctx->V,pep->nini));
    1119:   PetscCall(BVSetActiveColumns(ctx->V,0,1));
    1120:   if (k) {
    1121:     /* Insert deflated vectors */
    1122:     PetscCall(BVSetActiveColumns(pep->V,0,0));
    1123:     idx = sr->ndef0?sr->idxDef0[0]:sr->idxDef1[0];
    1124:     for (j=0;j<k;j++) {
    1125:       PetscCall(BVGetColumn(pep->V,j,&v));
    1126:       PetscCall(BVCopyVec(sr->V,sr->qinfo[idx].q[j],v));
    1127:       PetscCall(BVRestoreColumn(pep->V,j,&v));
    1128:     }
    1129:     /* Update innerproduct matrix */
    1130:     PetscCall(BVSetActiveColumns(ctx->V,0,0));
    1131:     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1132:     PetscCall(BVSetActiveColumns(pep->V,0,k));
    1133:     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));

    1135:     PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
    1136:     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1137:     PetscCall(MatDenseGetArray(MS,&S));
    1138:     for (j=0;j<sr->ndef0;j++) {
    1139:       PetscCall(PetscArrayzero(S+j*ld*deg,ld*deg));
    1140:       PetscCall(PetscArraycpy(S+j*ld*deg,sr->S+sr->idxDef0[j]*sr->ld*deg,k));
    1141:       PetscCall(PetscArraycpy(S+j*ld*deg+ld,sr->S+sr->idxDef0[j]*sr->ld*deg+sr->ld,k));
    1142:       pep->eigr[j] = sr->eigr[sr->idxDef0[j]];
    1143:       pep->errest[j] = sr->errest[sr->idxDef0[j]];
    1144:     }
    1145:     for (j=0;j<sr->ndef1;j++) {
    1146:       PetscCall(PetscArrayzero(S+(j+sr->ndef0)*ld*deg,ld*deg));
    1147:       PetscCall(PetscArraycpy(S+(j+sr->ndef0)*ld*deg,sr->S+sr->idxDef1[j]*sr->ld*deg,k));
    1148:       PetscCall(PetscArraycpy(S+(j+sr->ndef0)*ld*deg+ld,sr->S+sr->idxDef1[j]*sr->ld*deg+sr->ld,k));
    1149:       pep->eigr[j+sr->ndef0] = sr->eigr[sr->idxDef1[j]];
    1150:       pep->errest[j+sr->ndef0] = sr->errest[sr->idxDef1[j]];
    1151:     }
    1152:     PetscCall(MatDenseRestoreArray(MS,&S));
    1153:     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    1154:     PetscCall(BVSetActiveColumns(ctx->V,0,k+1));
    1155:     PetscCall(VecCreateSeq(PETSC_COMM_SELF,k+1,&vomega));
    1156:     PetscCall(VecGetArray(vomega,&om));
    1157:     for (j=0;j<k;j++) {
    1158:       PetscCall(BVOrthogonalizeColumn(ctx->V,j,NULL,&norm,NULL));
    1159:       PetscCall(BVScaleColumn(ctx->V,j,1/norm));
    1160:       om[j] = (norm>=0.0)?1.0:-1.0;
    1161:     }
    1162:     PetscCall(BVTensorGetFactors(ctx->V,NULL,&MS));
    1163:     PetscCall(MatDenseGetArray(MS,&S));
    1164:     for (j=0;j<deg;j++) {
    1165:       PetscCall(BVSetRandomColumn(pep->V,k+j));
    1166:       PetscCall(BVOrthogonalizeColumn(pep->V,k+j,S+k*ld*deg+j*ld,&norm,NULL));
    1167:       PetscCall(BVScaleColumn(pep->V,k+j,1.0/norm));
    1168:       S[k*ld*deg+j*ld+k+j] = norm;
    1169:     }
    1170:     PetscCall(MatDenseRestoreArray(MS,&S));
    1171:     PetscCall(BVSetActiveColumns(pep->V,0,k+deg));
    1172:     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));
    1173:     PetscCall(BVOrthogonalizeColumn(ctx->V,k,NULL,&norm,NULL));
    1174:     PetscCall(BVScaleColumn(ctx->V,k,1.0/norm));
    1175:     om[k] = (norm>=0.0)?1.0:-1.0;
    1176:     PetscCall(VecRestoreArray(vomega,&om));
    1177:     PetscCall(BVSetSignature(ctx->V,vomega));
    1178:     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&a));
    1179:     PetscCall(VecGetArray(vomega,&om));
    1180:     for (j=0;j<k;j++) a[j] = PetscRealPart(om[j]/(pep->eigr[j]-sigma));
    1181:     PetscCall(VecRestoreArray(vomega,&om));
    1182:     PetscCall(VecDestroy(&vomega));
    1183:     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&a));
    1184:     PetscCall(DSGetArray(pep->ds,DS_MAT_Q,&pQ));
    1185:     PetscCall(PetscArrayzero(pQ,ldds*k));
    1186:     for (j=0;j<k;j++) pQ[j+j*ldds] = 1.0;
    1187:     PetscCall(DSRestoreArray(pep->ds,DS_MAT_Q,&pQ));
    1188:   }
    1189:   PetscCall(BVSetActiveColumns(ctx->V,0,k+1));
    1190:   PetscCall(DSSetDimensions(pep->ds,k+1,PETSC_DETERMINE,PETSC_DETERMINE));
    1191:   PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
    1192:   PetscCall(BVGetSignature(ctx->V,vomega));
    1193:   PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));

    1195:   PetscCall(PetscInfo(pep,"Start STOAR: sigma=%g in [%g,%g], for deflation: left=%" PetscInt_FMT " right=%" PetscInt_FMT ", searching: left=%" PetscInt_FMT " right=%" PetscInt_FMT "\n",(double)sr->sPres->value,(double)(sr->sPres->neighb[0]?sr->sPres->neighb[0]->value:sr->int0),(double)(sr->sPres->neighb[1]?sr->sPres->neighb[1]->value:sr->int1),sr->ndef0,sr->ndef1,sr->sPres->nsch[0],sr->sPres->nsch[1]));

    1197:   /* Restart loop */
    1198:   l = 0;
    1199:   pep->nconv = k;
    1200:   while (pep->reason == PEP_CONVERGED_ITERATING) {
    1201:     pep->its++;
    1202:     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&a));
    1203:     b = a+ldds;
    1204:     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_D,&omega));

    1206:     /* Compute an nv-step Lanczos factorization */
    1207:     nv = PetscMin(pep->nconv+pep->mpd,pep->ncv);
    1208:     PetscCall(PEPSTOARrun_QSlice(pep,a,b,omega,pep->nconv+l,&nv,&breakdown,&symmlost,pep->work));
    1209:     beta = b[nv-1];
    1210:     if (symmlost && nv==pep->nconv+l) {
    1211:       pep->reason = PEP_DIVERGED_SYMMETRY_LOST;
    1212:       pep->nconv = nconv;
    1213:       PetscCall(PetscInfo(pep,"Symmetry lost in STOAR sigma=%g nconv=%" PetscInt_FMT "\n",(double)sr->sPres->value,nconv));
    1214:       if (falselock || !ctx->lock) {
    1215:         PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
    1216:         PetscCall(BVTensorCompress(ctx->V,0));
    1217:       }
    1218:       break;
    1219:     }
    1220:     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_T,&a));
    1221:     PetscCall(DSRestoreArrayReal(pep->ds,DS_MAT_D,&omega));
    1222:     PetscCall(DSSetDimensions(pep->ds,nv,pep->nconv,pep->nconv+l));
    1223:     if (l==0) PetscCall(DSSetState(pep->ds,DS_STATE_INTERMEDIATE));
    1224:     else PetscCall(DSSetState(pep->ds,DS_STATE_RAW));

    1226:     /* Solve projected problem */
    1227:     PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
    1228:     PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
    1229:     PetscCall(DSUpdateExtraRow(pep->ds));
    1230:     PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));

    1232:     /* Check convergence */
    1233:     /* PetscCall(PEPSTOARpreKConvergence(pep,nv,&norm,pep->work));*/
    1234:     norm = 1.0;
    1235:     PetscCall(DSGetDimensions(pep->ds,NULL,NULL,NULL,&t));
    1236:     PetscCall(PEPKrylovConvergence(pep,PETSC_FALSE,pep->nconv,t-pep->nconv,PetscAbsReal(beta)*norm,&k));
    1237:     PetscCall((*pep->stopping)(pep,pep->its,pep->max_it,k,pep->nev,&pep->reason,pep->stoppingctx));
    1238:     for (j=0;j<k;j++) back[j] = pep->eigr[j];
    1239:     PetscCall(STBackTransform(pep->st,k,back,pep->eigi));
    1240:     count0=count1=0;
    1241:     for (j=0;j<k;j++) {
    1242:       lambda = PetscRealPart(back[j]);
    1243:       if ((sr->dir*(sr->sPres->value - lambda) > 0) && (sr->dir*(lambda - sr->sPres->ext[0]) > 0)) count0++;
    1244:       if ((sr->dir*(lambda - sr->sPres->value) > 0) && (sr->dir*(sr->sPres->ext[1] - lambda) > 0)) count1++;
    1245:     }
    1246:     if ((count0-sr->ndef0 >= sr->sPres->nsch[0]) && (count1-sr->ndef1 >= sr->sPres->nsch[1])) pep->reason = PEP_CONVERGED_TOL;
    1247:     /* Update l */
    1248:     if (pep->reason != PEP_CONVERGED_ITERATING || breakdown) l = 0;
    1249:     else {
    1250:       l = PetscMax(1,(PetscInt)((nv-k)/2));
    1251:       l = PetscMin(l,t);
    1252:       PetscCall(DSGetTruncateSize(pep->ds,k,t,&l));
    1253:       if (!breakdown) {
    1254:         /* Prepare the Rayleigh quotient for restart */
    1255:         PetscCall(DSTruncate(pep->ds,k+l,PETSC_FALSE));
    1256:       }
    1257:     }
    1258:     nconv = k;
    1259:     if (!ctx->lock && pep->reason == PEP_CONVERGED_ITERATING && !breakdown) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
    1260:     if (l) PetscCall(PetscInfo(pep,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));

    1262:     /* Update S */
    1263:     PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&MQ));
    1264:     PetscCall(BVMultInPlace(ctx->V,MQ,pep->nconv,k+l));
    1265:     PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&MQ));

    1267:     /* Copy last column of S */
    1268:     PetscCall(BVCopyColumn(ctx->V,nv,k+l));
    1269:     PetscCall(BVSetActiveColumns(ctx->V,0,k+l));
    1270:     if (k+l) {
    1271:       PetscCall(DSSetDimensions(pep->ds,k+l,PETSC_DETERMINE,PETSC_DETERMINE));
    1272:       PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
    1273:       PetscCall(BVSetSignature(ctx->V,vomega));
    1274:       PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
    1275:     }

    1277:     if (breakdown && pep->reason == PEP_CONVERGED_ITERATING) {
    1278:       /* stop if breakdown */
    1279:       PetscCall(PetscInfo(pep,"Breakdown TOAR method (it=%" PetscInt_FMT " norm=%g)\n",pep->its,(double)beta));
    1280:       pep->reason = PEP_DIVERGED_BREAKDOWN;
    1281:     }
    1282:     if (pep->reason != PEP_CONVERGED_ITERATING) l--;
    1283:     PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
    1284:     if (k+l+deg<=nq) {
    1285:       PetscCall(BVSetActiveColumns(ctx->V,pep->nconv,k+l+1));
    1286:       if (!falselock && ctx->lock) PetscCall(BVTensorCompress(ctx->V,k-pep->nconv));
    1287:       else PetscCall(BVTensorCompress(ctx->V,0));
    1288:     }
    1289:     pep->nconv = k;
    1290:     PetscCall(PEPMonitor(pep,pep->its,nconv,pep->eigr,pep->eigi,pep->errest,nv));
    1291:   }
    1292:   sr->itsKs += pep->its;
    1293:   if (pep->nconv>0) {
    1294:     PetscCall(BVSetActiveColumns(ctx->V,0,pep->nconv));
    1295:     PetscCall(BVGetActiveColumns(pep->V,NULL,&nq));
    1296:     PetscCall(BVSetActiveColumns(pep->V,0,nq));
    1297:     if (nq>pep->nconv) {
    1298:       PetscCall(BVTensorCompress(ctx->V,pep->nconv));
    1299:       PetscCall(BVSetActiveColumns(pep->V,0,pep->nconv));
    1300:     }
    1301:     for (j=0;j<pep->nconv;j++) {
    1302:       pep->eigr[j] *= pep->sfactor;
    1303:       pep->eigi[j] *= pep->sfactor;
    1304:     }
    1305:   }
    1306:   PetscCall(PetscInfo(pep,"Finished STOAR: nconv=%" PetscInt_FMT " (deflated=%" PetscInt_FMT ", left=%" PetscInt_FMT ", right=%" PetscInt_FMT ")\n",pep->nconv,sr->ndef0+sr->ndef1,count0-sr->ndef0,count1-sr->ndef1));
    1307:   PetscCall(STScaleShift(pep->st,sinv?1.0/pep->sfactor:pep->sfactor));
    1308:   PetscCall(RGPopScale(pep->rg));

    1310:   PetscCheck(pep->reason!=PEP_DIVERGED_SYMMETRY_LOST || nconv>=sr->ndef0+sr->ndef1,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Symmetry lost at sigma=%g",(double)sr->sPres->value);
    1311:   if (pep->reason == PEP_DIVERGED_SYMMETRY_LOST && nconv==sr->ndef0+sr->ndef1) {
    1312:     PetscCheck(++sr->symmlost<=10,PetscObjectComm((PetscObject)pep),PETSC_ERR_PLIB,"Symmetry lost at sigma=%g",(double)sr->sPres->value);
    1313:   } else sr->symmlost = 0;

    1315:   PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
    1316:   PetscCall(PetscFree(back));
    1317:   PetscFunctionReturn(PETSC_SUCCESS);
    1318: }

    1320: static PetscErrorCode PEPQSliceGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
    1321: {
    1322:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
    1323:   PEP_SR          sr=ctx->sr;
    1324:   PetscInt        i=0,j,tmpi;
    1325:   PetscReal       v,tmpr;
    1326:   PEP_shift       s;

    1328:   PetscFunctionBegin;
    1329:   PetscCheck(pep->state,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Must call PEPSetUp() first");
    1330:   PetscCheck(sr,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see PEPSetInterval()");
    1331:   if (!sr->s0) {  /* PEPSolve not called yet */
    1332:     *n = 2;
    1333:   } else {
    1334:     *n = 1;
    1335:     s = sr->s0;
    1336:     while (s) {
    1337:       (*n)++;
    1338:       s = s->neighb[1];
    1339:     }
    1340:   }
    1341:   PetscCall(PetscMalloc1(*n,shifts));
    1342:   PetscCall(PetscMalloc1(*n,inertias));
    1343:   if (!sr->s0) {  /* PEPSolve not called yet */
    1344:     (*shifts)[0]   = sr->int0;
    1345:     (*shifts)[1]   = sr->int1;
    1346:     (*inertias)[0] = sr->inertia0;
    1347:     (*inertias)[1] = sr->inertia1;
    1348:   } else {
    1349:     s = sr->s0;
    1350:     while (s) {
    1351:       (*shifts)[i]     = s->value;
    1352:       (*inertias)[i++] = s->inertia;
    1353:       s = s->neighb[1];
    1354:     }
    1355:     (*shifts)[i]   = sr->int1;
    1356:     (*inertias)[i] = sr->inertia1;
    1357:   }
    1358:   /* remove possible duplicate in last position */
    1359:   if ((*shifts)[(*n)-1]==(*shifts)[(*n)-2]) (*n)--;
    1360:   /* sort result */
    1361:   for (i=0;i<*n;i++) {
    1362:     v = (*shifts)[i];
    1363:     for (j=i+1;j<*n;j++) {
    1364:       if (v > (*shifts)[j]) {
    1365:         SlepcSwap((*shifts)[i],(*shifts)[j],tmpr);
    1366:         SlepcSwap((*inertias)[i],(*inertias)[j],tmpi);
    1367:         v = (*shifts)[i];
    1368:       }
    1369:     }
    1370:   }
    1371:   PetscFunctionReturn(PETSC_SUCCESS);
    1372: }

    1374: PetscErrorCode PEPSolve_STOAR_QSlice(PEP pep)
    1375: {
    1376:   PetscInt       i,j,ti,deg=pep->nmat-1;
    1377:   PetscReal      newS;
    1378:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
    1379:   PEP_SR         sr=ctx->sr;
    1380:   Mat            S,B;
    1381:   PetscScalar    *pS;

    1383:   PetscFunctionBegin;
    1384:   PetscCall(PetscCitationsRegister(citation,&cited));

    1386:   /* Only with eigenvalues present in the interval ...*/
    1387:   if (sr->numEigs==0) {
    1388:     pep->reason = PEP_CONVERGED_TOL;
    1389:     PetscFunctionReturn(PETSC_SUCCESS);
    1390:   }

    1392:   /* Inner product matrix */
    1393:   PetscCall(PEPSTOARSetUpInnerMatrix(pep,&B));

    1395:   /* Array of pending shifts */
    1396:   sr->maxPend = 100; /* Initial size */
    1397:   sr->nPend = 0;
    1398:   PetscCall(PetscMalloc1(sr->maxPend,&sr->pending));
    1399:   PetscCall(PEPCreateShift(pep,sr->int0,NULL,NULL));
    1400:   /* extract first shift */
    1401:   sr->sPrev = NULL;
    1402:   sr->sPres = sr->pending[--sr->nPend];
    1403:   sr->sPres->inertia = sr->inertia0;
    1404:   pep->target = sr->sPres->value;
    1405:   sr->s0 = sr->sPres;
    1406:   sr->indexEig = 0;

    1408:   for (i=0;i<sr->numEigs;i++) {
    1409:     sr->eigr[i]   = 0.0;
    1410:     sr->eigi[i]   = 0.0;
    1411:     sr->errest[i] = 0.0;
    1412:     sr->perm[i]   = i;
    1413:   }
    1414:   /* Vectors for deflation */
    1415:   PetscCall(PetscMalloc2(sr->numEigs,&sr->idxDef0,sr->numEigs,&sr->idxDef1));
    1416:   sr->indexEig = 0;
    1417:   while (sr->sPres) {
    1418:     /* Search for deflation */
    1419:     PetscCall(PEPLookForDeflation(pep));
    1420:     /* KrylovSchur */
    1421:     PetscCall(PEPSTOAR_QSlice(pep,B));

    1423:     PetscCall(PEPStoreEigenpairs(pep));
    1424:     /* Select new shift */
    1425:     if (!sr->sPres->comp[1]) {
    1426:       PetscCall(PEPGetNewShiftValue(pep,1,&newS));
    1427:       PetscCall(PEPCreateShift(pep,newS,sr->sPres,sr->sPres->neighb[1]));
    1428:     }
    1429:     if (!sr->sPres->comp[0]) {
    1430:       /* Completing earlier interval */
    1431:       PetscCall(PEPGetNewShiftValue(pep,0,&newS));
    1432:       PetscCall(PEPCreateShift(pep,newS,sr->sPres->neighb[0],sr->sPres));
    1433:     }
    1434:     /* Preparing for a new search of values */
    1435:     PetscCall(PEPExtractShift(pep));
    1436:   }

    1438:   /* Updating pep values prior to exit */
    1439:   PetscCall(PetscFree2(sr->idxDef0,sr->idxDef1));
    1440:   PetscCall(PetscFree(sr->pending));
    1441:   pep->nconv  = sr->indexEig;
    1442:   pep->reason = PEP_CONVERGED_TOL;
    1443:   pep->its    = sr->itsKs;
    1444:   pep->nev    = sr->indexEig;
    1445:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,pep->nconv,pep->nconv,NULL,&S));
    1446:   PetscCall(MatDenseGetArray(S,&pS));
    1447:   for (i=0;i<pep->nconv;i++) {
    1448:     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);
    1449:   }
    1450:   PetscCall(MatDenseRestoreArray(S,&pS));
    1451:   PetscCall(BVSetActiveColumns(sr->V,0,pep->nconv));
    1452:   PetscCall(BVMultInPlace(sr->V,S,0,pep->nconv));
    1453:   PetscCall(MatDestroy(&S));
    1454:   PetscCall(BVDestroy(&pep->V));
    1455:   pep->V = sr->V;
    1456:   PetscCall(PetscFree4(pep->eigr,pep->eigi,pep->errest,pep->perm));
    1457:   pep->eigr   = sr->eigr;
    1458:   pep->eigi   = sr->eigi;
    1459:   pep->perm   = sr->perm;
    1460:   pep->errest = sr->errest;
    1461:   if (sr->dir<0) {
    1462:     for (i=0;i<pep->nconv/2;i++) {
    1463:       ti = sr->perm[i]; sr->perm[i] = sr->perm[pep->nconv-1-i]; sr->perm[pep->nconv-1-i] = ti;
    1464:     }
    1465:   }
    1466:   PetscCall(PetscFree(ctx->inertias));
    1467:   PetscCall(PetscFree(ctx->shifts));
    1468:   PetscCall(MatDestroy(&B));
    1469:   PetscCall(PEPQSliceGetInertias(pep,&ctx->nshifts,&ctx->shifts,&ctx->inertias));
    1470:   PetscFunctionReturn(PETSC_SUCCESS);
    1471: }