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