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: }