Actual source code: qslice.c

slepc-main 2024-11-09
Report Typos and Errors
  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(PCFactorGetMatrix(pc,&F));
152:   PetscCall(MatGetInertia(F,inertia,zeros,NULL));
153:   PetscFunctionReturn(PETSC_SUCCESS);
154: }

156: static PetscErrorCode PEPQSliceGetInertia(PEP pep,PetscReal shift,PetscInt *inertia,PetscInt *zeros,PetscInt correction)
157: {
158:   KSP            ksp;
159:   Mat            P;
160:   PetscReal      nzshift=0.0,dot;
161:   PetscRandom    rand;
162:   PetscInt       nconv;
163:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
164:   PEP_SR         sr=ctx->sr;

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

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

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

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

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

304: static inline PetscErrorCode PEPQSliceEvaluateQEP(PEP pep,PetscScalar x,Mat M,MatStructure str)
305: {
306:   PetscFunctionBegin;
307:   PetscCall(MatCopy(pep->A[0],M,SAME_NONZERO_PATTERN));
308:   PetscCall(MatAXPY(M,x,pep->A[1],str));
309:   PetscCall(MatAXPY(M,x*x,pep->A[2],str));
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

313: /*@
314:    PEPCheckDefiniteQEP - Determines if a symmetric/Hermitian quadratic eigenvalue problem
315:    is definite or not.

317:    Collective

319:    Input Parameter:
320: .  pep  - eigensolver context

322:    Output Parameters:
323: +  xi - first computed parameter
324: .  mu - second computed parameter
325: .  definite - flag indicating that the problem is definite
326: -  hyperbolic - flag indicating that the problem is hyperbolic

328:    Notes:
329:    This function is intended for quadratic eigenvalue problems, Q(lambda)=A*lambda^2+B*lambda+C,
330:    with symmetric (or Hermitian) coefficient matrices A,B,C.

332:    On output, the flag 'definite' may have the values -1 (meaning that the QEP is not
333:    definite), 1 (if the problem is definite), or 0 if the algorithm was not able to
334:    determine whether the problem is definite or not.

336:    If definite=1, the output flag 'hyperbolic' informs in a similar way about whether the
337:    problem is hyperbolic or not.

339:    If definite=1, the computed values xi and mu satisfy Q(xi)<0 and Q(mu)>0, as
340:    obtained via the method proposed in [Niendorf and Voss, LAA 2010]. Furthermore, if
341:    hyperbolic=1 then only xi is computed.

343:    Level: advanced

345: .seealso: PEPSetProblemType()
346: @*/
347: PetscErrorCode PEPCheckDefiniteQEP(PEP pep,PetscReal *xi,PetscReal *mu,PetscInt *definite,PetscInt *hyperbolic)
348: {
349:   PetscRandom    rand;
350:   Vec            u,w;
351:   PetscReal      d=0.0,s=0.0,sp,mut=0.0,omg=0.0,omgp;
352:   PetscInt       k,its=10,hyp=0,check=0,nconv,inertia,n;
353:   Mat            M=NULL;
354:   MatStructure   str;
355:   EPS            eps;
356:   PetscBool      transform,ptypehyp;

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

466: /*
467:    Dummy backtransform operation
468:  */
469: static PetscErrorCode PEPBackTransform_Skip(PEP pep)
470: {
471:   PetscFunctionBegin;
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: PetscErrorCode PEPSetUp_STOAR_QSlice(PEP pep)
476: {
477:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
478:   PEP_SR         sr;
479:   PetscInt       ld,i,zeros=0;
480:   SlepcSC        sc;
481:   PetscReal      r;

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

501:   /* create spectrum slicing context and initialize it */
502:   PetscCall(PEPQSliceResetSR(pep));
503:   PetscCall(PetscNew(&sr));
504:   ctx->sr   = sr;
505:   sr->itsKs = 0;
506:   sr->nleap = 0;
507:   sr->sPres = NULL;

509:   if (pep->solvematcoeffs) PetscCall(PetscFree(pep->solvematcoeffs));
510:   PetscCall(PetscMalloc1(pep->nmat,&pep->solvematcoeffs));
511:   if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
512:   PetscCall(STSetTransform(pep->st,PETSC_FALSE));
513:   PetscCall(STSetUp(pep->st));

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

517:   /* check presence of ends and finding direction */
518:   if (pep->inta > PETSC_MIN_REAL || pep->intb >= PETSC_MAX_REAL) {
519:     sr->int0 = pep->inta;
520:     sr->int1 = pep->intb;
521:     sr->dir = 1;
522:     if (pep->intb >= PETSC_MAX_REAL) { /* Right-open interval */
523:       sr->hasEnd = PETSC_FALSE;
524:     } else sr->hasEnd = PETSC_TRUE;
525:   } else {
526:     sr->int0 = pep->intb;
527:     sr->int1 = pep->inta;
528:     sr->dir = -1;
529:     sr->hasEnd = PetscNot(pep->inta <= PETSC_MIN_REAL);
530:   }

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

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

554:   if (sr->hasEnd) {
555:     sr->dir = -sr->dir; r = sr->int0; sr->int0 = sr->int1; sr->int1 = r;
556:     i = sr->inertia0; sr->inertia0 = sr->inertia1; sr->inertia1 = i;
557:   }

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

581: /*
582:    Fills the fields of a shift structure
583: */
584: static PetscErrorCode PEPCreateShift(PEP pep,PetscReal val,PEP_shift neighb0,PEP_shift neighb1)
585: {
586:   PEP_shift      s,*pending2;
587:   PetscInt       i;
588:   PEP_SR         sr;
589:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;

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

618: /* Provides next shift to be computed */
619: static PetscErrorCode PEPExtractShift(PEP pep)
620: {
621:   PetscInt       iner,zeros=0;
622:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
623:   PEP_SR         sr;
624:   PetscReal      newShift,aux;
625:   PEP_shift      sPres;

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

665: /*
666:   Obtains value of subsequent shift
667: */
668: static PetscErrorCode PEPGetNewShiftValue(PEP pep,PetscInt side,PetscReal *newS)
669: {
670:   PetscReal lambda,d_prev;
671:   PetscInt  i,idxP;
672:   PEP_SR    sr;
673:   PEP_shift sPres,s;
674:   PEP_STOAR *ctx=(PEP_STOAR*)pep->data;

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

747: /*
748:   Function for sorting an array of real values
749: */
750: static PetscErrorCode sortRealEigenvalues(PetscScalar *r,PetscInt *perm,PetscInt nr,PetscBool prev,PetscInt dir)
751: {
752:   PetscReal re;
753:   PetscInt  i,j,tmp;

755:   PetscFunctionBegin;
756:   if (!prev) for (i=0;i<nr;i++) perm[i] = i;
757:   /* Insertion sort */
758:   for (i=1;i<nr;i++) {
759:     re = PetscRealPart(r[perm[i]]);
760:     j = i-1;
761:     while (j>=0 && dir*(re - PetscRealPart(r[perm[j]])) <= 0) {
762:       tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
763:     }
764:   }
765:   PetscFunctionReturn(PETSC_SUCCESS);
766: }

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

782:   PetscFunctionBegin;
783:   sPres = sr->sPres;
784:   sPres->index = sr->indexEig;

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

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

933: static PetscErrorCode PEPLookForDeflation(PEP pep)
934: {
935:   PetscReal val;
936:   PetscInt  i,count0=0,count1=0;
937:   PEP_shift sPres;
938:   PetscInt  ini,fin;
939:   PEP_SR    sr;
940:   PEP_STOAR *ctx=(PEP_STOAR*)pep->data;

942:   PetscFunctionBegin;
943:   sr = ctx->sr;
944:   sPres = sr->sPres;

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

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

976:   /* Completing vector of indexes for deflation */
977:   for (i=0;i<count0;i++) sr->idxDef0[i] = sr->perm[ini+i];
978:   sr->ndef0 = count0;
979:   for (i=0;i<count1;i++) sr->idxDef1[i] = sr->perm[ini+count0+i];
980:   sr->ndef1 = count1;
981:   PetscFunctionReturn(PETSC_SUCCESS);
982: }

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

999:   PetscFunctionBegin;
1000:   PetscCall(PetscMalloc1(*M,&y));
1001:   PetscCall(BVGetSizes(pep->V,NULL,NULL,&ld));
1002:   PetscCall(BVTensorGetDegree(ctx->V,&d));
1003:   PetscCall(BVGetActiveColumns(pep->V,&lock,&nqt));
1004:   lds = d*ld;
1005:   offq = ld;
1006:   PetscCall(DSGetLeadingDimension(pep->ds,&ldds));

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

1030:     /* orthogonalize */
1031:     PetscCall(BVOrthogonalizeColumn(pep->V,nqt,S+(j+1)*lds,&norm,&lindep));
1032:     if (!lindep) {
1033:       *(S+(j+1)*lds+nqt) = norm;
1034:       PetscCall(BVScaleColumn(pep->V,nqt,1.0/norm));
1035:       nqt++;
1036:     }
1037:     for (i=0;i<nqt;i++) *(S+(j+1)*lds+offq+i) = *(S+j*lds+i)+sigma*(*(S+(j+1)*lds+i));
1038:     PetscCall(BVSetActiveColumns(pep->V,0,nqt));
1039:     PetscCall(MatDenseRestoreArray(MS,&S));
1040:     PetscCall(BVTensorRestoreFactors(ctx->V,NULL,&MS));

1042:     /* level-2 orthogonalization */
1043:     PetscCall(BVOrthogonalizeColumn(ctx->V,j+1,y,&norm,&lindep));
1044:     a[j] = PetscRealPart(y[j]);
1045:     omega[j+1] = (norm > 0)?1.0:-1.0;
1046:     PetscCall(BVScaleColumn(ctx->V,j+1,1.0/norm));
1047:     b[j] = PetscAbsReal(norm);

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

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

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

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

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

1189:   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]));

1191:   /* Restart loop */
1192:   l = 0;
1193:   pep->nconv = k;
1194:   while (pep->reason == PEP_CONVERGED_ITERATING) {
1195:     pep->its++;
1196:     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_T,&a));
1197:     b = a+ldds;
1198:     PetscCall(DSGetArrayReal(pep->ds,DS_MAT_D,&omega));

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

1220:     /* Solve projected problem */
1221:     PetscCall(DSSolve(pep->ds,pep->eigr,pep->eigi));
1222:     PetscCall(DSSort(pep->ds,pep->eigr,pep->eigi,NULL,NULL,NULL));
1223:     PetscCall(DSUpdateExtraRow(pep->ds));
1224:     PetscCall(DSSynchronize(pep->ds,pep->eigr,pep->eigi));

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

1256:     /* Update S */
1257:     PetscCall(DSGetMat(pep->ds,DS_MAT_Q,&MQ));
1258:     PetscCall(BVMultInPlace(ctx->V,MQ,pep->nconv,k+l));
1259:     PetscCall(DSRestoreMat(pep->ds,DS_MAT_Q,&MQ));

1261:     /* Copy last column of S */
1262:     PetscCall(BVCopyColumn(ctx->V,nv,k+l));
1263:     PetscCall(BVSetActiveColumns(ctx->V,0,k+l));
1264:     if (k+l) {
1265:       PetscCall(DSSetDimensions(pep->ds,k+l,PETSC_DETERMINE,PETSC_DETERMINE));
1266:       PetscCall(DSGetMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
1267:       PetscCall(BVSetSignature(ctx->V,vomega));
1268:       PetscCall(DSRestoreMatAndColumn(pep->ds,DS_MAT_D,0,&D,&vomega));
1269:     }

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

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

1309:   PetscCall(DSTruncate(pep->ds,pep->nconv,PETSC_TRUE));
1310:   PetscCall(PetscFree(back));
1311:   PetscFunctionReturn(PETSC_SUCCESS);
1312: }

1314: static PetscErrorCode PEPQSliceGetInertias(PEP pep,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
1315: {
1316:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
1317:   PEP_SR          sr=ctx->sr;
1318:   PetscInt        i=0,j,tmpi;
1319:   PetscReal       v,tmpr;
1320:   PEP_shift       s;

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

1368: PetscErrorCode PEPSolve_STOAR_QSlice(PEP pep)
1369: {
1370:   PetscInt       i,j,ti,deg=pep->nmat-1;
1371:   PetscReal      newS;
1372:   PEP_STOAR      *ctx=(PEP_STOAR*)pep->data;
1373:   PEP_SR         sr=ctx->sr;
1374:   Mat            S,B;
1375:   PetscScalar    *pS;

1377:   PetscFunctionBegin;
1378:   PetscCall(PetscCitationsRegister(citation,&cited));

1380:   /* Only with eigenvalues present in the interval ...*/
1381:   if (sr->numEigs==0) {
1382:     pep->reason = PEP_CONVERGED_TOL;
1383:     PetscFunctionReturn(PETSC_SUCCESS);
1384:   }

1386:   /* Inner product matrix */
1387:   PetscCall(PEPSTOARSetUpInnerMatrix(pep,&B));

1389:   /* Array of pending shifts */
1390:   sr->maxPend = 100; /* Initial size */
1391:   sr->nPend = 0;
1392:   PetscCall(PetscMalloc1(sr->maxPend,&sr->pending));
1393:   PetscCall(PEPCreateShift(pep,sr->int0,NULL,NULL));
1394:   /* extract first shift */
1395:   sr->sPrev = NULL;
1396:   sr->sPres = sr->pending[--sr->nPend];
1397:   sr->sPres->inertia = sr->inertia0;
1398:   pep->target = sr->sPres->value;
1399:   sr->s0 = sr->sPres;
1400:   sr->indexEig = 0;

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

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

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