Actual source code: ks-slice.c
slepc-3.21.1 2024-04-26
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 eigensolver: "krylovschur"
13: Method: Krylov-Schur with spectrum slicing for symmetric eigenproblems
15: References:
17: [1] R.G. Grimes et al., "A shifted block Lanczos algorithm for
18: solving sparse symmetric generalized eigenproblems", SIAM J.
19: Matrix Anal. Appl. 15(1):228-272, 1994.
21: [2] C. Campos and J.E. Roman, "Spectrum slicing strategies based
22: on restarted Lanczos methods", Numer. Algor. 60(2):279-295,
23: 2012.
24: */
26: #include <slepc/private/epsimpl.h>
27: #include "krylovschur.h"
29: static PetscBool cited = PETSC_FALSE;
30: static const char citation[] =
31: "@Article{slepc-slice,\n"
32: " author = \"C. Campos and J. E. Roman\",\n"
33: " title = \"Strategies for spectrum slicing based on restarted {Lanczos} methods\",\n"
34: " journal = \"Numer. Algorithms\",\n"
35: " volume = \"60\",\n"
36: " number = \"2\",\n"
37: " pages = \"279--295\",\n"
38: " year = \"2012,\"\n"
39: " doi = \"https://doi.org/10.1007/s11075-012-9564-z\"\n"
40: "}\n";
42: #define SLICE_PTOL PETSC_SQRT_MACHINE_EPSILON
44: static PetscErrorCode EPSSliceResetSR(EPS eps)
45: {
46: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
47: EPS_SR sr=ctx->sr;
48: EPS_shift s;
50: PetscFunctionBegin;
51: if (sr) {
52: if (ctx->npart>1) {
53: PetscCall(BVDestroy(&sr->V));
54: PetscCall(PetscFree4(sr->eigr,sr->eigi,sr->errest,sr->perm));
55: }
56: /* Reviewing list of shifts to free memory */
57: s = sr->s0;
58: if (s) {
59: while (s->neighb[1]) {
60: s = s->neighb[1];
61: PetscCall(PetscFree(s->neighb[0]));
62: }
63: PetscCall(PetscFree(s));
64: }
65: PetscCall(PetscFree(sr));
66: }
67: ctx->sr = NULL;
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: PetscErrorCode EPSReset_KrylovSchur_Slice(EPS eps)
72: {
73: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
75: PetscFunctionBegin;
76: if (!ctx->global) PetscFunctionReturn(PETSC_SUCCESS);
77: /* Reset auxiliary EPS */
78: PetscCall(EPSSliceResetSR(ctx->eps));
79: PetscCall(EPSReset(ctx->eps));
80: PetscCall(EPSSliceResetSR(eps));
81: PetscCall(PetscFree(ctx->inertias));
82: PetscCall(PetscFree(ctx->shifts));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: PetscErrorCode EPSDestroy_KrylovSchur_Slice(EPS eps)
87: {
88: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
90: PetscFunctionBegin;
91: if (!ctx->global) PetscFunctionReturn(PETSC_SUCCESS);
92: /* Destroy auxiliary EPS */
93: PetscCall(EPSReset_KrylovSchur_Slice(eps));
94: PetscCall(EPSDestroy(&ctx->eps));
95: if (ctx->npart>1) {
96: PetscCall(PetscSubcommDestroy(&ctx->subc));
97: if (ctx->commset) {
98: PetscCallMPI(MPI_Comm_free(&ctx->commrank));
99: ctx->commset = PETSC_FALSE;
100: }
101: PetscCall(ISDestroy(&ctx->isrow));
102: PetscCall(ISDestroy(&ctx->iscol));
103: PetscCall(MatDestroyMatrices(1,&ctx->submata));
104: PetscCall(MatDestroyMatrices(1,&ctx->submatb));
105: }
106: PetscCall(PetscFree(ctx->subintervals));
107: PetscCall(PetscFree(ctx->nconv_loc));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
111: /*
112: EPSSliceAllocateSolution - Allocate memory storage for common variables such
113: as eigenvalues and eigenvectors. The argument extra is used for methods
114: that require a working basis slightly larger than ncv.
115: */
116: static PetscErrorCode EPSSliceAllocateSolution(EPS eps,PetscInt extra)
117: {
118: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
119: PetscReal eta;
120: PetscInt k;
121: BVType type;
122: BVOrthogType orthog_type;
123: BVOrthogRefineType orthog_ref;
124: BVOrthogBlockType ob_type;
125: Mat matrix;
126: Vec t;
127: EPS_SR sr = ctx->sr;
129: PetscFunctionBegin;
130: /* allocate space for eigenvalues and friends */
131: k = PetscMax(1,sr->numEigs);
132: PetscCall(PetscFree4(sr->eigr,sr->eigi,sr->errest,sr->perm));
133: PetscCall(PetscMalloc4(k,&sr->eigr,k,&sr->eigi,k,&sr->errest,k,&sr->perm));
135: /* allocate sr->V and transfer options from eps->V */
136: PetscCall(BVDestroy(&sr->V));
137: PetscCall(BVCreate(PetscObjectComm((PetscObject)eps),&sr->V));
138: if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
139: if (!((PetscObject)eps->V)->type_name) PetscCall(BVSetType(sr->V,BVMAT));
140: else {
141: PetscCall(BVGetType(eps->V,&type));
142: PetscCall(BVSetType(sr->V,type));
143: }
144: PetscCall(STMatCreateVecsEmpty(eps->st,&t,NULL));
145: PetscCall(BVSetSizesFromVec(sr->V,t,k));
146: PetscCall(VecDestroy(&t));
147: PetscCall(EPS_SetInnerProduct(eps));
148: PetscCall(BVGetMatrix(eps->V,&matrix,NULL));
149: PetscCall(BVSetMatrix(sr->V,matrix,PETSC_FALSE));
150: PetscCall(BVGetOrthogonalization(eps->V,&orthog_type,&orthog_ref,&eta,&ob_type));
151: PetscCall(BVSetOrthogonalization(sr->V,orthog_type,orthog_ref,eta,ob_type));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode EPSSliceGetEPS(EPS eps)
156: {
157: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
158: BV V;
159: BVType type;
160: PetscReal eta;
161: BVOrthogType orthog_type;
162: BVOrthogRefineType orthog_ref;
163: BVOrthogBlockType ob_type;
164: PetscInt i;
165: PetscReal h,a,b;
166: PetscRandom rand;
167: EPS_SR sr=ctx->sr;
169: PetscFunctionBegin;
170: if (!ctx->eps) PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
172: /* Determine subintervals */
173: if (ctx->npart==1) {
174: a = eps->inta; b = eps->intb;
175: } else {
176: if (!ctx->subintset) { /* uniform distribution if no set by user */
177: PetscCheck(sr->hasEnd,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Global interval must be bounded for splitting it in uniform subintervals");
178: h = (eps->intb-eps->inta)/ctx->npart;
179: a = eps->inta+ctx->subc->color*h;
180: b = (ctx->subc->color==ctx->npart-1)?eps->intb:eps->inta+(ctx->subc->color+1)*h;
181: PetscCall(PetscFree(ctx->subintervals));
182: PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
183: for (i=0;i<ctx->npart;i++) ctx->subintervals[i] = eps->inta+h*i;
184: ctx->subintervals[ctx->npart] = eps->intb;
185: } else {
186: a = ctx->subintervals[ctx->subc->color];
187: b = ctx->subintervals[ctx->subc->color+1];
188: }
189: }
190: PetscCall(EPSSetInterval(ctx->eps,a,b));
191: PetscCall(EPSSetConvergenceTest(ctx->eps,eps->conv));
192: PetscCall(EPSSetDimensions(ctx->eps,ctx->nev,ctx->ncv,ctx->mpd));
193: PetscCall(EPSKrylovSchurSetLocking(ctx->eps,ctx->lock));
195: ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
196: ctx_local->detect = ctx->detect;
198: /* transfer options from eps->V */
199: PetscCall(EPSGetBV(ctx->eps,&V));
200: PetscCall(BVGetRandomContext(V,&rand)); /* make sure the random context is available when duplicating */
201: if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
202: if (!((PetscObject)eps->V)->type_name) PetscCall(BVSetType(V,BVMAT));
203: else {
204: PetscCall(BVGetType(eps->V,&type));
205: PetscCall(BVSetType(V,type));
206: }
207: PetscCall(BVGetOrthogonalization(eps->V,&orthog_type,&orthog_ref,&eta,&ob_type));
208: PetscCall(BVSetOrthogonalization(V,orthog_type,orthog_ref,eta,ob_type));
210: ctx->eps->which = eps->which;
211: ctx->eps->max_it = eps->max_it;
212: ctx->eps->tol = eps->tol;
213: ctx->eps->purify = eps->purify;
214: if (eps->tol==(PetscReal)PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL;
215: PetscCall(EPSSetProblemType(ctx->eps,eps->problem_type));
216: PetscCall(EPSSetUp(ctx->eps));
217: ctx->eps->nconv = 0;
218: ctx->eps->its = 0;
219: for (i=0;i<ctx->eps->ncv;i++) {
220: ctx->eps->eigr[i] = 0.0;
221: ctx->eps->eigi[i] = 0.0;
222: ctx->eps->errest[i] = 0.0;
223: }
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: static PetscErrorCode EPSSliceGetInertia(EPS eps,PetscReal shift,PetscInt *inertia,PetscInt *zeros)
228: {
229: KSP ksp,kspr;
230: PC pc;
231: Mat F;
232: PetscReal nzshift=shift;
233: PetscBool flg;
235: PetscFunctionBegin;
236: if (shift >= PETSC_MAX_REAL) { /* Right-open interval */
237: if (inertia) *inertia = eps->n;
238: } else if (shift <= PETSC_MIN_REAL) {
239: if (inertia) *inertia = 0;
240: if (zeros) *zeros = 0;
241: } else {
242: /* If the shift is zero, perturb it to a very small positive value.
243: The goal is that the nonzero pattern is the same in all cases and reuse
244: the symbolic factorizations */
245: nzshift = (shift==0.0)? 10.0/PETSC_MAX_REAL: shift;
246: PetscCall(STSetShift(eps->st,nzshift));
247: PetscCall(STGetKSP(eps->st,&ksp));
248: PetscCall(KSPGetPC(ksp,&pc));
249: PetscCall(PetscObjectTypeCompare((PetscObject)pc,PCREDUNDANT,&flg));
250: if (flg) {
251: PetscCall(PCRedundantGetKSP(pc,&kspr));
252: PetscCall(KSPGetPC(kspr,&pc));
253: }
254: PetscCall(PCFactorGetMatrix(pc,&F));
255: PetscCall(MatGetInertia(F,inertia,zeros,NULL));
256: }
257: if (inertia) PetscCall(PetscInfo(eps,"Computed inertia at shift %g: %" PetscInt_FMT "\n",(double)nzshift,*inertia));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: /*
262: Dummy backtransform operation
263: */
264: static PetscErrorCode EPSBackTransform_Skip(EPS eps)
265: {
266: PetscFunctionBegin;
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: PetscErrorCode EPSSetUp_KrylovSchur_Slice(EPS eps)
271: {
272: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*ctx_glob;
273: EPS_SR sr,sr_loc,sr_glob;
274: PetscInt nEigs,dssz=1,i,zeros=0,off=0,method,hiteig=0;
275: PetscMPIInt nproc,rank=0,aux;
276: PetscReal r;
277: MPI_Request req;
278: Mat A,B=NULL;
279: DSParallelType ptype;
280: MPI_Comm child;
282: PetscFunctionBegin;
283: if (ctx->global) {
284: EPSCheckHermitianDefiniteCondition(eps,PETSC_TRUE," with spectrum slicing");
285: EPSCheckSinvertCayleyCondition(eps,PETSC_TRUE," with spectrum slicing");
286: PetscCheck(eps->inta!=eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing all eigenvalues unless you provide a computational interval with EPSSetInterval()");
287: PetscCheck(eps->intb<PETSC_MAX_REAL || eps->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The defined computational interval should have at least one of their sides bounded");
288: PetscCheck(eps->nds==0,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectrum slicing not supported in combination with deflation space");
289: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_STOPPING,PETSC_TRUE," with spectrum slicing");
290: EPSCheckIgnoredCondition(eps,EPS_FEATURE_BALANCE,PETSC_TRUE," with spectrum slicing");
291: if (eps->tol==(PetscReal)PETSC_DEFAULT) {
292: #if defined(PETSC_USE_REAL_SINGLE)
293: eps->tol = SLEPC_DEFAULT_TOL;
294: #else
295: /* use tighter tolerance */
296: eps->tol = SLEPC_DEFAULT_TOL*1e-2;
297: #endif
298: }
299: if (eps->max_it==PETSC_DEFAULT) eps->max_it = 100;
300: if (ctx->nev==1) ctx->nev = PetscMin(40,eps->n); /* nev not set, use default value */
301: PetscCheck(eps->n<=10 || ctx->nev>=10,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"nev cannot be less than 10 in spectrum slicing runs");
302: }
303: eps->ops->backtransform = EPSBackTransform_Skip;
305: /* create spectrum slicing context and initialize it */
306: PetscCall(EPSSliceResetSR(eps));
307: PetscCall(PetscNew(&sr));
308: ctx->sr = sr;
309: sr->itsKs = 0;
310: sr->nleap = 0;
311: sr->nMAXCompl = eps->nev/4;
312: sr->iterCompl = eps->max_it/4;
313: sr->sPres = NULL;
314: sr->nS = 0;
316: if (ctx->npart==1 || ctx->global) {
317: /* check presence of ends and finding direction */
318: if ((eps->inta > PETSC_MIN_REAL && !(ctx->subintervals && ctx->subintervals[0]==ctx->subintervals[1])) || eps->intb >= PETSC_MAX_REAL) {
319: sr->int0 = eps->inta;
320: sr->int1 = eps->intb;
321: sr->dir = 1;
322: if (eps->intb >= PETSC_MAX_REAL) { /* Right-open interval */
323: sr->hasEnd = PETSC_FALSE;
324: } else sr->hasEnd = PETSC_TRUE;
325: } else {
326: sr->int0 = eps->intb;
327: sr->int1 = eps->inta;
328: sr->dir = -1;
329: sr->hasEnd = PetscNot(eps->inta <= PETSC_MIN_REAL);
330: }
331: }
332: if (ctx->global) {
333: PetscCall(EPSSetDimensions_Default(eps,ctx->nev,&ctx->ncv,&ctx->mpd));
334: /* create subintervals and initialize auxiliary eps for slicing runs */
335: PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
336: /* prevent computation of factorization in global eps */
337: PetscCall(STSetTransform(eps->st,PETSC_FALSE));
338: PetscCall(EPSSliceGetEPS(eps));
339: sr_loc = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
340: if (ctx->npart>1) {
341: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
342: if ((sr->dir>0&&ctx->subc->color==0)||(sr->dir<0&&ctx->subc->color==ctx->npart-1)) sr->inertia0 = sr_loc->inertia0;
343: PetscCallMPI(MPI_Comm_rank(child,&rank));
344: if (!rank) PetscCallMPI(MPI_Bcast(&sr->inertia0,1,MPIU_INT,(sr->dir>0)?0:ctx->npart-1,ctx->commrank));
345: PetscCallMPI(MPI_Bcast(&sr->inertia0,1,MPIU_INT,0,child));
346: PetscCall(PetscFree(ctx->nconv_loc));
347: PetscCall(PetscMalloc1(ctx->npart,&ctx->nconv_loc));
348: PetscCallMPI(MPI_Comm_size(((PetscObject)eps)->comm,&nproc));
349: if (sr->dir<0) off = 1;
350: if (nproc%ctx->npart==0) { /* subcommunicators with the same size */
351: PetscCall(PetscMPIIntCast(sr_loc->numEigs,&aux));
352: PetscCallMPI(MPI_Allgather(&aux,1,MPI_INT,ctx->nconv_loc,1,MPI_INT,ctx->commrank));
353: PetscCallMPI(MPI_Allgather(sr_loc->dir==sr->dir?&sr_loc->int0:&sr_loc->int1,1,MPIU_REAL,ctx->subintervals+off,1,MPIU_REAL,ctx->commrank));
354: } else {
355: PetscCallMPI(MPI_Comm_rank(child,&rank));
356: if (!rank) {
357: PetscCall(PetscMPIIntCast(sr_loc->numEigs,&aux));
358: PetscCallMPI(MPI_Allgather(&aux,1,MPI_INT,ctx->nconv_loc,1,MPI_INT,ctx->commrank));
359: PetscCallMPI(MPI_Allgather(sr_loc->dir==sr->dir?&sr_loc->int0:&sr_loc->int1,1,MPIU_REAL,ctx->subintervals+off,1,MPIU_REAL,ctx->commrank));
360: }
361: PetscCall(PetscMPIIntCast(ctx->npart,&aux));
362: PetscCallMPI(MPI_Bcast(ctx->nconv_loc,aux,MPI_INT,0,child));
363: PetscCallMPI(MPI_Bcast(ctx->subintervals+off,aux,MPIU_REAL,0,child));
364: }
365: nEigs = 0;
366: for (i=0;i<ctx->npart;i++) nEigs += ctx->nconv_loc[i];
367: } else {
368: nEigs = sr_loc->numEigs;
369: sr->inertia0 = sr_loc->inertia0;
370: sr->dir = sr_loc->dir;
371: }
372: sr->inertia1 = sr->inertia0+sr->dir*nEigs;
373: sr->numEigs = nEigs;
374: eps->nev = nEigs;
375: eps->ncv = nEigs;
376: eps->mpd = nEigs;
377: } else {
378: ctx_glob = (EPS_KRYLOVSCHUR*)ctx->eps->data;
379: sr_glob = ctx_glob->sr;
380: if (ctx->npart>1) {
381: sr->dir = sr_glob->dir;
382: sr->int0 = (sr->dir==1)?eps->inta:eps->intb;
383: sr->int1 = (sr->dir==1)?eps->intb:eps->inta;
384: if ((sr->dir>0&&ctx->subc->color==ctx->npart-1)||(sr->dir<0&&ctx->subc->color==0)) sr->hasEnd = sr_glob->hasEnd;
385: else sr->hasEnd = PETSC_TRUE;
386: }
387: /* sets first shift */
388: PetscCall(STSetShift(eps->st,(sr->int0==0.0)?10.0/PETSC_MAX_REAL:sr->int0));
389: PetscCall(STSetUp(eps->st));
391: /* compute inertia0 */
392: PetscCall(EPSSliceGetInertia(eps,sr->int0,&sr->inertia0,ctx->detect?&zeros:NULL));
393: /* undocumented option to control what to do when an eigenvalue is found:
394: - error out if it's the endpoint of the user-provided interval (or sub-interval)
395: - if it's an endpoint computed internally:
396: + if hiteig=0 error out
397: + else if hiteig=1 the subgroup that hit the eigenvalue does nothing
398: + otherwise the subgroup that hit the eigenvalue perturbs the shift and recomputes inertia
399: */
400: PetscCall(PetscOptionsGetInt(NULL,NULL,"-eps_krylovschur_hiteigenvalue",&hiteig,NULL));
401: if (zeros) { /* error in factorization */
402: PetscCheck(sr->int0!=ctx->eps->inta && sr->int0!=ctx->eps->intb,((PetscObject)eps)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in the interval endpoint");
403: PetscCheck(!ctx_glob->subintset || hiteig,((PetscObject)eps)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in an interval endpoint defined by user");
404: if (hiteig==1) { /* idle subgroup */
405: sr->inertia0 = -1;
406: } else { /* perturb shift */
407: sr->int0 *= (1.0+SLICE_PTOL);
408: PetscCall(EPSSliceGetInertia(eps,sr->int0,&sr->inertia0,&zeros));
409: PetscCheck(zeros==0,((PetscObject)eps)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)sr->int1);
410: }
411: }
412: if (ctx->npart>1) {
413: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
414: /* inertia1 is received from neighbour */
415: PetscCallMPI(MPI_Comm_rank(child,&rank));
416: if (!rank) {
417: if (sr->inertia0!=-1 && ((sr->dir>0 && ctx->subc->color>0) || (sr->dir<0 && ctx->subc->color<ctx->npart-1))) { /* send inertia0 to neighbour0 */
418: PetscCallMPI(MPI_Isend(&sr->inertia0,1,MPIU_INT,ctx->subc->color-sr->dir,0,ctx->commrank,&req));
419: PetscCallMPI(MPI_Isend(&sr->int0,1,MPIU_REAL,ctx->subc->color-sr->dir,0,ctx->commrank,&req));
420: }
421: if ((sr->dir>0 && ctx->subc->color<ctx->npart-1)|| (sr->dir<0 && ctx->subc->color>0)) { /* receive inertia1 from neighbour1 */
422: PetscCallMPI(MPI_Recv(&sr->inertia1,1,MPIU_INT,ctx->subc->color+sr->dir,0,ctx->commrank,MPI_STATUS_IGNORE));
423: PetscCallMPI(MPI_Recv(&sr->int1,1,MPIU_REAL,ctx->subc->color+sr->dir,0,ctx->commrank,MPI_STATUS_IGNORE));
424: }
425: if (sr->inertia0==-1 && !(sr->dir>0 && ctx->subc->color==ctx->npart-1) && !(sr->dir<0 && ctx->subc->color==0)) {
426: sr->inertia0 = sr->inertia1; sr->int0 = sr->int1;
427: PetscCallMPI(MPI_Isend(&sr->inertia0,1,MPIU_INT,ctx->subc->color-sr->dir,0,ctx->commrank,&req));
428: PetscCallMPI(MPI_Isend(&sr->int0,1,MPIU_REAL,ctx->subc->color-sr->dir,0,ctx->commrank,&req));
429: }
430: }
431: if ((sr->dir>0 && ctx->subc->color<ctx->npart-1)||(sr->dir<0 && ctx->subc->color>0)) {
432: PetscCallMPI(MPI_Bcast(&sr->inertia1,1,MPIU_INT,0,child));
433: PetscCallMPI(MPI_Bcast(&sr->int1,1,MPIU_REAL,0,child));
434: } else sr_glob->inertia1 = sr->inertia1;
435: }
437: /* last process in eps comm computes inertia1 */
438: if (ctx->npart==1 || ((sr->dir>0 && ctx->subc->color==ctx->npart-1) || (sr->dir<0 && ctx->subc->color==0))) {
439: PetscCall(EPSSliceGetInertia(eps,sr->int1,&sr->inertia1,ctx->detect?&zeros:NULL));
440: PetscCheck(zeros==0,((PetscObject)eps)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in an interval endpoint defined by user");
441: if (!rank && sr->inertia0==-1) {
442: sr->inertia0 = sr->inertia1; sr->int0 = sr->int1;
443: PetscCallMPI(MPI_Isend(&sr->inertia0,1,MPIU_INT,ctx->subc->color-sr->dir,0,ctx->commrank,&req));
444: PetscCallMPI(MPI_Isend(&sr->int0,1,MPIU_REAL,ctx->subc->color-sr->dir,0,ctx->commrank,&req));
445: }
446: if (sr->hasEnd) {
447: sr->dir = -sr->dir; r = sr->int0; sr->int0 = sr->int1; sr->int1 = r;
448: i = sr->inertia0; sr->inertia0 = sr->inertia1; sr->inertia1 = i;
449: }
450: }
452: /* number of eigenvalues in interval */
453: sr->numEigs = (sr->dir)*(sr->inertia1 - sr->inertia0);
454: if (ctx->npart>1) {
455: /* memory allocate for subinterval eigenpairs */
456: PetscCall(EPSSliceAllocateSolution(eps,1));
457: }
458: dssz = eps->ncv+1;
459: PetscCall(DSGetParallel(ctx->eps->ds,&ptype));
460: PetscCall(DSSetParallel(eps->ds,ptype));
461: PetscCall(DSGetMethod(ctx->eps->ds,&method));
462: PetscCall(DSSetMethod(eps->ds,method));
463: }
464: PetscCall(DSSetType(eps->ds,DSHEP));
465: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
466: PetscCall(DSAllocate(eps->ds,dssz));
467: /* keep state of subcomm matrices to check that the user does not modify them */
468: PetscCall(EPSGetOperators(eps,&A,&B));
469: PetscCall(PetscObjectStateGet((PetscObject)A,&ctx->Astate));
470: PetscCall(PetscObjectGetId((PetscObject)A,&ctx->Aid));
471: if (B) {
472: PetscCall(PetscObjectStateGet((PetscObject)B,&ctx->Bstate));
473: PetscCall(PetscObjectGetId((PetscObject)B,&ctx->Bid));
474: } else {
475: ctx->Bstate=0;
476: ctx->Bid=0;
477: }
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: static PetscErrorCode EPSSliceGatherEigenVectors(EPS eps)
482: {
483: Vec v,vg,v_loc;
484: IS is1,is2;
485: VecScatter vec_sc;
486: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
487: PetscInt nloc,m0,n0,i,si,idx,*idx1,*idx2,j;
488: PetscScalar *array;
489: EPS_SR sr_loc;
490: BV V_loc;
492: PetscFunctionBegin;
493: sr_loc = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
494: V_loc = sr_loc->V;
496: /* Gather parallel eigenvectors */
497: PetscCall(BVGetColumn(eps->V,0,&v));
498: PetscCall(VecGetOwnershipRange(v,&n0,&m0));
499: PetscCall(BVRestoreColumn(eps->V,0,&v));
500: PetscCall(BVGetColumn(ctx->eps->V,0,&v));
501: PetscCall(VecGetLocalSize(v,&nloc));
502: PetscCall(BVRestoreColumn(ctx->eps->V,0,&v));
503: PetscCall(PetscMalloc2(m0-n0,&idx1,m0-n0,&idx2));
504: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)eps),nloc,PETSC_DECIDE,&vg));
505: idx = -1;
506: for (si=0;si<ctx->npart;si++) {
507: j = 0;
508: for (i=n0;i<m0;i++) {
509: idx1[j] = i;
510: idx2[j++] = i+eps->n*si;
511: }
512: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)eps),(m0-n0),idx1,PETSC_COPY_VALUES,&is1));
513: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)eps),(m0-n0),idx2,PETSC_COPY_VALUES,&is2));
514: PetscCall(BVGetColumn(eps->V,0,&v));
515: PetscCall(VecScatterCreate(v,is1,vg,is2,&vec_sc));
516: PetscCall(BVRestoreColumn(eps->V,0,&v));
517: PetscCall(ISDestroy(&is1));
518: PetscCall(ISDestroy(&is2));
519: for (i=0;i<ctx->nconv_loc[si];i++) {
520: PetscCall(BVGetColumn(eps->V,++idx,&v));
521: if (ctx->subc->color==si) {
522: PetscCall(BVGetColumn(V_loc,i,&v_loc));
523: PetscCall(VecGetArray(v_loc,&array));
524: PetscCall(VecPlaceArray(vg,array));
525: }
526: PetscCall(VecScatterBegin(vec_sc,vg,v,INSERT_VALUES,SCATTER_REVERSE));
527: PetscCall(VecScatterEnd(vec_sc,vg,v,INSERT_VALUES,SCATTER_REVERSE));
528: if (ctx->subc->color==si) {
529: PetscCall(VecResetArray(vg));
530: PetscCall(VecRestoreArray(v_loc,&array));
531: PetscCall(BVRestoreColumn(V_loc,i,&v_loc));
532: }
533: PetscCall(BVRestoreColumn(eps->V,idx,&v));
534: }
535: PetscCall(VecScatterDestroy(&vec_sc));
536: }
537: PetscCall(PetscFree2(idx1,idx2));
538: PetscCall(VecDestroy(&vg));
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: /*
543: EPSComputeVectors_Slice - Recover Eigenvectors from subcomunicators
544: */
545: PetscErrorCode EPSComputeVectors_Slice(EPS eps)
546: {
547: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
549: PetscFunctionBegin;
550: if (ctx->global && ctx->npart>1) {
551: PetscCall(EPSComputeVectors(ctx->eps));
552: PetscCall(EPSSliceGatherEigenVectors(eps));
553: }
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: #define SWAP(a,b,t) do {t=a;a=b;b=t;} while (0)
559: static PetscErrorCode EPSSliceGetInertias(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
560: {
561: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
562: PetscInt i=0,j,tmpi;
563: PetscReal v,tmpr;
564: EPS_shift s;
566: PetscFunctionBegin;
567: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
568: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
569: if (!ctx->sr->s0) { /* EPSSolve not called yet */
570: *n = 2;
571: } else {
572: *n = 1;
573: s = ctx->sr->s0;
574: while (s) {
575: (*n)++;
576: s = s->neighb[1];
577: }
578: }
579: PetscCall(PetscMalloc1(*n,shifts));
580: PetscCall(PetscMalloc1(*n,inertias));
581: if (!ctx->sr->s0) { /* EPSSolve not called yet */
582: (*shifts)[0] = ctx->sr->int0;
583: (*shifts)[1] = ctx->sr->int1;
584: (*inertias)[0] = ctx->sr->inertia0;
585: (*inertias)[1] = ctx->sr->inertia1;
586: } else {
587: s = ctx->sr->s0;
588: while (s) {
589: (*shifts)[i] = s->value;
590: (*inertias)[i++] = s->inertia;
591: s = s->neighb[1];
592: }
593: (*shifts)[i] = ctx->sr->int1;
594: (*inertias)[i] = ctx->sr->inertia1;
595: }
596: /* remove possible duplicate in last position */
597: if ((*shifts)[(*n)-1]==(*shifts)[(*n)-2]) (*n)--;
598: /* sort result */
599: for (i=0;i<*n;i++) {
600: v = (*shifts)[i];
601: for (j=i+1;j<*n;j++) {
602: if (v > (*shifts)[j]) {
603: SWAP((*shifts)[i],(*shifts)[j],tmpr);
604: SWAP((*inertias)[i],(*inertias)[j],tmpi);
605: v = (*shifts)[i];
606: }
607: }
608: }
609: PetscFunctionReturn(PETSC_SUCCESS);
610: }
612: static PetscErrorCode EPSSliceGatherSolution(EPS eps)
613: {
614: PetscMPIInt rank,nproc,*disp,*ns_loc,aux;
615: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
616: PetscInt i,idx,j,*perm_loc,off=0,*inertias_loc,ns;
617: PetscScalar *eigr_loc;
618: EPS_SR sr_loc;
619: PetscReal *shifts_loc;
620: MPI_Comm child;
622: PetscFunctionBegin;
623: eps->nconv = 0;
624: for (i=0;i<ctx->npart;i++) eps->nconv += ctx->nconv_loc[i];
625: sr_loc = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
627: /* Gather the shifts used and the inertias computed */
628: PetscCall(EPSSliceGetInertias(ctx->eps,&ns,&shifts_loc,&inertias_loc));
629: if (ctx->sr->dir>0 && shifts_loc[ns-1]==sr_loc->int1 && ctx->subc->color<ctx->npart-1) ns--;
630: if (ctx->sr->dir<0 && shifts_loc[ns-1]==sr_loc->int0 && ctx->subc->color>0) {
631: ns--;
632: for (i=0;i<ns;i++) {
633: inertias_loc[i] = inertias_loc[i+1];
634: shifts_loc[i] = shifts_loc[i+1];
635: }
636: }
637: PetscCall(PetscMalloc1(ctx->npart,&ns_loc));
638: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
639: PetscCallMPI(MPI_Comm_rank(child,&rank));
640: PetscCall(PetscMPIIntCast(ns,&aux));
641: if (!rank) PetscCallMPI(MPI_Allgather(&aux,1,MPI_INT,ns_loc,1,MPI_INT,ctx->commrank));
642: PetscCall(PetscMPIIntCast(ctx->npart,&aux));
643: PetscCallMPI(MPI_Bcast(ns_loc,aux,MPI_INT,0,child));
644: ctx->nshifts = 0;
645: for (i=0;i<ctx->npart;i++) ctx->nshifts += ns_loc[i];
646: PetscCall(PetscFree(ctx->inertias));
647: PetscCall(PetscFree(ctx->shifts));
648: PetscCall(PetscMalloc1(ctx->nshifts,&ctx->inertias));
649: PetscCall(PetscMalloc1(ctx->nshifts,&ctx->shifts));
651: /* Gather eigenvalues (same ranks have fully set of eigenvalues)*/
652: eigr_loc = sr_loc->eigr;
653: perm_loc = sr_loc->perm;
654: PetscCallMPI(MPI_Comm_size(((PetscObject)eps)->comm,&nproc));
655: PetscCall(PetscMalloc1(ctx->npart,&disp));
656: disp[0] = 0;
657: for (i=1;i<ctx->npart;i++) disp[i] = disp[i-1]+ctx->nconv_loc[i-1];
658: if (nproc%ctx->npart==0) { /* subcommunicators with the same size */
659: PetscCall(PetscMPIIntCast(sr_loc->numEigs,&aux));
660: PetscCallMPI(MPI_Allgatherv(eigr_loc,aux,MPIU_SCALAR,eps->eigr,ctx->nconv_loc,disp,MPIU_SCALAR,ctx->commrank)); /* eigenvalues */
661: PetscCallMPI(MPI_Allgatherv(perm_loc,aux,MPIU_INT,eps->perm,ctx->nconv_loc,disp,MPIU_INT,ctx->commrank)); /* perm */
662: for (i=1;i<ctx->npart;i++) disp[i] = disp[i-1]+ns_loc[i-1];
663: PetscCall(PetscMPIIntCast(ns,&aux));
664: PetscCallMPI(MPI_Allgatherv(shifts_loc,aux,MPIU_REAL,ctx->shifts,ns_loc,disp,MPIU_REAL,ctx->commrank)); /* shifts */
665: PetscCallMPI(MPI_Allgatherv(inertias_loc,aux,MPIU_INT,ctx->inertias,ns_loc,disp,MPIU_INT,ctx->commrank)); /* inertias */
666: PetscCall(MPIU_Allreduce(&sr_loc->itsKs,&eps->its,1,MPIU_INT,MPI_SUM,ctx->commrank));
667: } else { /* subcommunicators with different size */
668: if (!rank) {
669: PetscCall(PetscMPIIntCast(sr_loc->numEigs,&aux));
670: PetscCallMPI(MPI_Allgatherv(eigr_loc,aux,MPIU_SCALAR,eps->eigr,ctx->nconv_loc,disp,MPIU_SCALAR,ctx->commrank)); /* eigenvalues */
671: PetscCallMPI(MPI_Allgatherv(perm_loc,aux,MPIU_INT,eps->perm,ctx->nconv_loc,disp,MPIU_INT,ctx->commrank)); /* perm */
672: for (i=1;i<ctx->npart;i++) disp[i] = disp[i-1]+ns_loc[i-1];
673: PetscCall(PetscMPIIntCast(ns,&aux));
674: PetscCallMPI(MPI_Allgatherv(shifts_loc,aux,MPIU_REAL,ctx->shifts,ns_loc,disp,MPIU_REAL,ctx->commrank)); /* shifts */
675: PetscCallMPI(MPI_Allgatherv(inertias_loc,aux,MPIU_INT,ctx->inertias,ns_loc,disp,MPIU_INT,ctx->commrank)); /* inertias */
676: PetscCall(MPIU_Allreduce(&sr_loc->itsKs,&eps->its,1,MPIU_INT,MPI_SUM,ctx->commrank));
677: }
678: PetscCall(PetscMPIIntCast(eps->nconv,&aux));
679: PetscCallMPI(MPI_Bcast(eps->eigr,aux,MPIU_SCALAR,0,child));
680: PetscCallMPI(MPI_Bcast(eps->perm,aux,MPIU_INT,0,child));
681: PetscCallMPI(MPI_Bcast(ctx->shifts,ctx->nshifts,MPIU_REAL,0,child));
682: PetscCall(PetscMPIIntCast(ctx->nshifts,&aux));
683: PetscCallMPI(MPI_Bcast(ctx->inertias,aux,MPIU_INT,0,child));
684: PetscCallMPI(MPI_Bcast(&eps->its,1,MPIU_INT,0,child));
685: }
686: /* Update global array eps->perm */
687: idx = ctx->nconv_loc[0];
688: for (i=1;i<ctx->npart;i++) {
689: off += ctx->nconv_loc[i-1];
690: for (j=0;j<ctx->nconv_loc[i];j++) eps->perm[idx++] += off;
691: }
693: /* Gather parallel eigenvectors */
694: PetscCall(PetscFree(ns_loc));
695: PetscCall(PetscFree(disp));
696: PetscCall(PetscFree(shifts_loc));
697: PetscCall(PetscFree(inertias_loc));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: /*
702: Fills the fields of a shift structure
703: */
704: static PetscErrorCode EPSCreateShift(EPS eps,PetscReal val,EPS_shift neighb0,EPS_shift neighb1)
705: {
706: EPS_shift s,*pending2;
707: PetscInt i;
708: EPS_SR sr;
709: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
711: PetscFunctionBegin;
712: sr = ctx->sr;
713: if ((neighb0 && val==neighb0->value) || (neighb1 && val==neighb1->value)) {
714: sr->nPend++;
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
717: PetscCall(PetscNew(&s));
718: s->value = val;
719: s->neighb[0] = neighb0;
720: if (neighb0) neighb0->neighb[1] = s;
721: s->neighb[1] = neighb1;
722: if (neighb1) neighb1->neighb[0] = s;
723: s->comp[0] = PETSC_FALSE;
724: s->comp[1] = PETSC_FALSE;
725: s->index = -1;
726: s->neigs = 0;
727: s->nconv[0] = s->nconv[1] = 0;
728: s->nsch[0] = s->nsch[1]=0;
729: /* Inserts in the stack of pending shifts */
730: /* If needed, the array is resized */
731: if (sr->nPend >= sr->maxPend) {
732: sr->maxPend *= 2;
733: PetscCall(PetscMalloc1(sr->maxPend,&pending2));
734: for (i=0;i<sr->nPend;i++) pending2[i] = sr->pending[i];
735: PetscCall(PetscFree(sr->pending));
736: sr->pending = pending2;
737: }
738: sr->pending[sr->nPend++]=s;
739: PetscFunctionReturn(PETSC_SUCCESS);
740: }
742: /* Prepare for Rational Krylov update */
743: static PetscErrorCode EPSPrepareRational(EPS eps)
744: {
745: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
746: PetscInt dir,i,k,ld,nv;
747: PetscScalar *A;
748: EPS_SR sr = ctx->sr;
749: Vec v;
751: PetscFunctionBegin;
752: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
753: dir = (sr->sPres->neighb[0] == sr->sPrev)?1:-1;
754: dir*=sr->dir;
755: k = 0;
756: for (i=0;i<sr->nS;i++) {
757: if (dir*PetscRealPart(sr->S[i])>0.0) {
758: sr->S[k] = sr->S[i];
759: sr->S[sr->nS+k] = sr->S[sr->nS+i];
760: PetscCall(BVGetColumn(sr->Vnext,k,&v));
761: PetscCall(BVCopyVec(eps->V,eps->nconv+i,v));
762: PetscCall(BVRestoreColumn(sr->Vnext,k,&v));
763: k++;
764: if (k>=sr->nS/2)break;
765: }
766: }
767: /* Copy to DS */
768: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&A));
769: PetscCall(PetscArrayzero(A,ld*ld));
770: for (i=0;i<k;i++) {
771: A[i*(1+ld)] = sr->S[i];
772: A[k+i*ld] = sr->S[sr->nS+i];
773: }
774: sr->nS = k;
775: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&A));
776: PetscCall(DSGetDimensions(eps->ds,&nv,NULL,NULL,NULL));
777: PetscCall(DSSetDimensions(eps->ds,nv,0,k));
778: /* Append u to V */
779: PetscCall(BVGetColumn(sr->Vnext,sr->nS,&v));
780: PetscCall(BVCopyVec(eps->V,sr->nv,v));
781: PetscCall(BVRestoreColumn(sr->Vnext,sr->nS,&v));
782: PetscFunctionReturn(PETSC_SUCCESS);
783: }
785: /* Provides next shift to be computed */
786: static PetscErrorCode EPSExtractShift(EPS eps)
787: {
788: PetscInt iner,zeros=0;
789: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
790: EPS_SR sr;
791: PetscReal newShift,diam,ptol;
792: EPS_shift sPres;
794: PetscFunctionBegin;
795: sr = ctx->sr;
796: if (sr->nPend > 0) {
797: if (sr->sPres==sr->pending[sr->nPend-1]) {
798: eps->reason = EPS_CONVERGED_ITERATING;
799: eps->its = 0;
800: sr->nPend--;
801: sr->sPres->rep = PETSC_TRUE;
802: PetscFunctionReturn(PETSC_SUCCESS);
803: }
804: sr->sPrev = sr->sPres;
805: sr->sPres = sr->pending[--sr->nPend];
806: sPres = sr->sPres;
807: PetscCall(EPSSliceGetInertia(eps,sPres->value,&iner,ctx->detect?&zeros:NULL));
808: if (zeros) {
809: diam = PetscMin(PetscAbsReal(sPres->neighb[0]->value-sPres->value),PetscAbsReal(sPres->neighb[1]->value-sPres->value));
810: ptol = PetscMin(SLICE_PTOL,diam/2);
811: newShift = sPres->value*(1.0+ptol);
812: if (sr->dir*(sPres->neighb[0] && newShift-sPres->neighb[0]->value) < 0) newShift = (sPres->value+sPres->neighb[0]->value)/2;
813: else if (sPres->neighb[1] && sr->dir*(sPres->neighb[1]->value-newShift) < 0) newShift = (sPres->value+sPres->neighb[1]->value)/2;
814: PetscCall(EPSSliceGetInertia(eps,newShift,&iner,&zeros));
815: PetscCheck(zeros==0,((PetscObject)eps)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)newShift);
816: sPres->value = newShift;
817: }
818: sr->sPres->inertia = iner;
819: eps->target = sr->sPres->value;
820: eps->reason = EPS_CONVERGED_ITERATING;
821: eps->its = 0;
822: } else sr->sPres = NULL;
823: PetscFunctionReturn(PETSC_SUCCESS);
824: }
826: /*
827: Symmetric KrylovSchur adapted to spectrum slicing:
828: Allows searching an specific amount of eigenvalues in the subintervals left and right.
829: Returns whether the search has succeeded
830: */
831: static PetscErrorCode EPSKrylovSchur_Slice(EPS eps)
832: {
833: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
834: PetscInt i,k,l,ld,nv,*iwork,j,count0,count1,iterCompl=0,n0,n1;
835: Mat U,Op,T;
836: PetscScalar *Q,*A;
837: PetscReal *a,*b,beta,lambda;
838: EPS_shift sPres;
839: PetscBool breakdown,complIterating,sch0,sch1;
840: EPS_SR sr = ctx->sr;
842: PetscFunctionBegin;
843: /* Spectrum slicing data */
844: sPres = sr->sPres;
845: complIterating =PETSC_FALSE;
846: sch1 = sch0 = PETSC_TRUE;
847: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
848: PetscCall(PetscMalloc1(2*ld,&iwork));
849: count0=0;count1=0; /* Found on both sides */
850: if (!sPres->rep && sr->nS > 0 && (sPres->neighb[0] == sr->sPrev || sPres->neighb[1] == sr->sPrev)) {
851: /* Rational Krylov */
852: PetscCall(DSTranslateRKS(eps->ds,sr->sPrev->value-sPres->value));
853: PetscCall(DSGetDimensions(eps->ds,NULL,NULL,&l,NULL));
854: PetscCall(DSSetDimensions(eps->ds,l+1,0,0));
855: PetscCall(BVSetActiveColumns(eps->V,0,l+1));
856: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
857: PetscCall(BVMultInPlace(eps->V,U,0,l+1));
858: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
859: } else {
860: /* Get the starting Lanczos vector */
861: PetscCall(EPSGetStartVector(eps,0,NULL));
862: l = 0;
863: }
864: /* Restart loop */
865: while (eps->reason == EPS_CONVERGED_ITERATING) {
866: eps->its++; sr->itsKs++;
867: /* Compute an nv-step Lanczos factorization */
868: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
869: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
870: PetscCall(DSGetMat(eps->ds,DS_MAT_T,&T));
871: PetscCall(STGetOperator(eps->st,&Op));
872: PetscCall(BVMatLanczos(eps->V,Op,T,eps->nconv+l,&nv,&beta,&breakdown));
873: PetscCall(STRestoreOperator(eps->st,&Op));
874: sr->nv = nv;
875: PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&T));
876: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
877: if (l==0) PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
878: else PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
879: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
881: /* Solve projected problem and compute residual norm estimates */
882: if (eps->its == 1 && l > 0) {/* After rational update */
883: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&A));
884: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
885: b = a + ld;
886: k = eps->nconv+l;
887: A[k*ld+k-1] = A[(k-1)*ld+k];
888: A[k*ld+k] = a[k];
889: for (j=k+1; j< nv; j++) {
890: A[j*ld+j] = a[j];
891: A[j*ld+j-1] = b[j-1] ;
892: A[(j-1)*ld+j] = b[j-1];
893: }
894: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&A));
895: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
896: PetscCall(DSSolve(eps->ds,eps->eigr,NULL));
897: PetscCall(DSSort(eps->ds,eps->eigr,NULL,NULL,NULL,NULL));
898: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
899: } else { /* Restart */
900: PetscCall(DSSolve(eps->ds,eps->eigr,NULL));
901: PetscCall(DSSort(eps->ds,eps->eigr,NULL,NULL,NULL,NULL));
902: }
903: PetscCall(DSSynchronize(eps->ds,eps->eigr,NULL));
905: /* Residual */
906: PetscCall(EPSKrylovConvergence(eps,PETSC_TRUE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
907: /* Checking values obtained for completing */
908: for (i=0;i<k;i++) {
909: sr->back[i]=eps->eigr[i];
910: }
911: PetscCall(STBackTransform(eps->st,k,sr->back,eps->eigi));
912: count0=count1=0;
913: for (i=0;i<k;i++) {
914: lambda = PetscRealPart(sr->back[i]);
915: if ((sr->dir*(sPres->value - lambda) > 0) && (sr->dir*(lambda - sPres->ext[0]) > 0)) count0++;
916: if ((sr->dir*(lambda - sPres->value) > 0) && (sr->dir*(sPres->ext[1] - lambda) > 0)) count1++;
917: }
918: if (k>eps->nev && eps->ncv-k<5) eps->reason = EPS_CONVERGED_TOL;
919: else {
920: /* Checks completion */
921: if ((!sch0||count0 >= sPres->nsch[0]) && (!sch1 ||count1 >= sPres->nsch[1])) {
922: eps->reason = EPS_CONVERGED_TOL;
923: } else {
924: if (!complIterating && eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
925: if (complIterating) {
926: if (--iterCompl <= 0) eps->reason = EPS_DIVERGED_ITS;
927: } else if (k >= eps->nev) {
928: n0 = sPres->nsch[0]-count0;
929: n1 = sPres->nsch[1]-count1;
930: if (sr->iterCompl>0 && ((n0>0 && n0<= sr->nMAXCompl)||(n1>0&&n1<=sr->nMAXCompl))) {
931: /* Iterating for completion*/
932: complIterating = PETSC_TRUE;
933: if (n0 >sr->nMAXCompl)sch0 = PETSC_FALSE;
934: if (n1 >sr->nMAXCompl)sch1 = PETSC_FALSE;
935: iterCompl = sr->iterCompl;
936: } else eps->reason = EPS_CONVERGED_TOL;
937: }
938: }
939: }
940: /* Update l */
941: if (eps->reason == EPS_CONVERGED_ITERATING) l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
942: else l = nv-k;
943: if (breakdown) l=0;
944: if (!ctx->lock && l>0 && eps->reason == EPS_CONVERGED_ITERATING) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
946: if (eps->reason == EPS_CONVERGED_ITERATING) {
947: if (breakdown) {
948: /* Start a new Lanczos factorization */
949: PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
950: PetscCall(EPSGetStartVector(eps,k,&breakdown));
951: if (breakdown) {
952: eps->reason = EPS_DIVERGED_BREAKDOWN;
953: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
954: }
955: } else {
956: /* Prepare the Rayleigh quotient for restart */
957: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
958: PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Q));
959: b = a + ld;
960: for (i=k;i<k+l;i++) {
961: a[i] = PetscRealPart(eps->eigr[i]);
962: b[i] = PetscRealPart(Q[nv-1+i*ld]*beta);
963: }
964: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
965: PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Q));
966: }
967: }
968: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
969: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
970: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
971: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
973: /* Normalize u and append it to V */
974: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(eps->V,nv,k+l));
975: eps->nconv = k;
976: if (eps->reason != EPS_CONVERGED_ITERATING) {
977: /* Store approximated values for next shift */
978: PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Q));
979: sr->nS = l;
980: for (i=0;i<l;i++) {
981: sr->S[i] = eps->eigr[i+k];/* Diagonal elements */
982: sr->S[i+l] = Q[nv-1+(i+k)*ld]*beta; /* Out of diagonal elements */
983: }
984: PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Q));
985: }
986: }
987: /* Check for completion */
988: for (i=0;i< eps->nconv; i++) {
989: if (sr->dir*PetscRealPart(eps->eigr[i])>0) sPres->nconv[1]++;
990: else sPres->nconv[0]++;
991: }
992: sPres->comp[0] = PetscNot(count0 < sPres->nsch[0]);
993: sPres->comp[1] = PetscNot(count1 < sPres->nsch[1]);
994: PetscCall(PetscInfo(eps,"Lanczos: %" PetscInt_FMT " evals in [%g,%g]%s and %" PetscInt_FMT " evals in [%g,%g]%s\n",count0,(double)(sr->dir==1?sPres->ext[0]:sPres->value),(double)(sr->dir==1?sPres->value:sPres->ext[0]),sPres->comp[0]?"*":"",count1,(double)(sr->dir==1?sPres->value:sPres->ext[1]),(double)(sr->dir==1?sPres->ext[1]:sPres->value),sPres->comp[1]?"*":""));
995: PetscCheck(count0<=sPres->nsch[0] && count1<=sPres->nsch[1],PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia%s",ctx->detect?"":", consider using EPSKrylovSchurSetDetectZeros()");
996: PetscCall(PetscFree(iwork));
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: /*
1001: Obtains value of subsequent shift
1002: */
1003: static PetscErrorCode EPSGetNewShiftValue(EPS eps,PetscInt side,PetscReal *newS)
1004: {
1005: PetscReal lambda,d_prev;
1006: PetscInt i,idxP;
1007: EPS_SR sr;
1008: EPS_shift sPres,s;
1009: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1011: PetscFunctionBegin;
1012: sr = ctx->sr;
1013: sPres = sr->sPres;
1014: if (sPres->neighb[side]) {
1015: /* Completing a previous interval */
1016: *newS = (sPres->value + sPres->neighb[side]->value)/2;
1017: if (PetscAbsReal(sPres->value - *newS)/PetscAbsReal(sPres->value)<=100*PETSC_SQRT_MACHINE_EPSILON) *newS = sPres->value;
1018: } else { /* (Only for side=1). Creating a new interval. */
1019: if (sPres->neigs==0) {/* No value has been accepted*/
1020: if (sPres->neighb[0]) {
1021: /* Multiplying by 10 the previous distance */
1022: *newS = sPres->value + 10*sr->dir*PetscAbsReal(sPres->value - sPres->neighb[0]->value);
1023: sr->nleap++;
1024: /* Stops when the interval is open and no values are found in the last 5 shifts (there might be infinite eigenvalues) */
1025: PetscCheck(sr->hasEnd || sr->nleap<=5,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unable to compute the wanted eigenvalues with open interval");
1026: } else { /* First shift */
1027: PetscCheck(eps->nconv!=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"First shift renders no information");
1028: /* Unaccepted values give information for next shift */
1029: idxP=0;/* Number of values left from shift */
1030: for (i=0;i<eps->nconv;i++) {
1031: lambda = PetscRealPart(eps->eigr[i]);
1032: if (sr->dir*(lambda - sPres->value) <0) idxP++;
1033: else break;
1034: }
1035: /* Avoiding subtraction of eigenvalues (might be the same).*/
1036: if (idxP>0) {
1037: d_prev = PetscAbsReal(sPres->value - PetscRealPart(eps->eigr[0]))/(idxP+0.3);
1038: } else {
1039: d_prev = PetscAbsReal(sPres->value - PetscRealPart(eps->eigr[eps->nconv-1]))/(eps->nconv+0.3);
1040: }
1041: *newS = sPres->value + (sr->dir*d_prev*eps->nev)/2;
1042: }
1043: } else { /* Accepted values found */
1044: sr->nleap = 0;
1045: /* Average distance of values in previous subinterval */
1046: s = sPres->neighb[0];
1047: while (s && PetscAbs(s->inertia - sPres->inertia)==0) {
1048: s = s->neighb[0];/* Looking for previous shifts with eigenvalues within */
1049: }
1050: if (s) {
1051: d_prev = PetscAbsReal((sPres->value - s->value)/(sPres->inertia - s->inertia));
1052: } else { /* First shift. Average distance obtained with values in this shift */
1053: /* first shift might be too far from first wanted eigenvalue (no values found outside the interval)*/
1054: 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(eps->tol)) {
1055: d_prev = PetscAbsReal((PetscRealPart(sr->eigr[sr->indexEig-1]) - PetscRealPart(sr->eigr[0])))/(sPres->neigs+0.3);
1056: } else {
1057: d_prev = PetscAbsReal(PetscRealPart(sr->eigr[sr->indexEig-1]) - sPres->value)/(sPres->neigs+0.3);
1058: }
1059: }
1060: /* Average distance is used for next shift by adding it to value on the right or to shift */
1061: if (sr->dir*(PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1]) - sPres->value)>0) {
1062: *newS = PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1])+ (sr->dir*d_prev*eps->nev)/2;
1063: } else { /* Last accepted value is on the left of shift. Adding to shift */
1064: *newS = sPres->value + (sr->dir*d_prev*eps->nev)/2;
1065: }
1066: }
1067: /* End of interval can not be surpassed */
1068: if (sr->dir*(sr->int1 - *newS) < 0) *newS = sr->int1;
1069: }/* of neighb[side]==null */
1070: PetscFunctionReturn(PETSC_SUCCESS);
1071: }
1073: /*
1074: Function for sorting an array of real values
1075: */
1076: static PetscErrorCode sortRealEigenvalues(PetscScalar *r,PetscInt *perm,PetscInt nr,PetscBool prev,PetscInt dir)
1077: {
1078: PetscReal re;
1079: PetscInt i,j,tmp;
1081: PetscFunctionBegin;
1082: if (!prev) for (i=0;i<nr;i++) perm[i] = i;
1083: /* Insertion sort */
1084: for (i=1;i<nr;i++) {
1085: re = PetscRealPart(r[perm[i]]);
1086: j = i-1;
1087: while (j>=0 && dir*(re - PetscRealPart(r[perm[j]])) <= 0) {
1088: tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
1089: }
1090: }
1091: PetscFunctionReturn(PETSC_SUCCESS);
1092: }
1094: /* Stores the pairs obtained since the last shift in the global arrays */
1095: static PetscErrorCode EPSStoreEigenpairs(EPS eps)
1096: {
1097: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1098: PetscReal lambda,err,norm;
1099: PetscInt i,count;
1100: PetscBool iscayley;
1101: EPS_SR sr = ctx->sr;
1102: EPS_shift sPres;
1103: Vec v,w;
1105: PetscFunctionBegin;
1106: sPres = sr->sPres;
1107: sPres->index = sr->indexEig;
1108: count = sr->indexEig;
1109: /* Back-transform */
1110: PetscCall(STBackTransform(eps->st,eps->nconv,eps->eigr,eps->eigi));
1111: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley));
1112: /* Sort eigenvalues */
1113: PetscCall(sortRealEigenvalues(eps->eigr,eps->perm,eps->nconv,PETSC_FALSE,sr->dir));
1114: /* Values stored in global array */
1115: for (i=0;i<eps->nconv;i++) {
1116: lambda = PetscRealPart(eps->eigr[eps->perm[i]]);
1117: err = eps->errest[eps->perm[i]];
1119: if (sr->dir*(lambda - sPres->ext[0]) > 0 && (sr->dir)*(sPres->ext[1] - lambda) > 0) {/* Valid value */
1120: PetscCheck(count<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error in Spectrum Slicing");
1121: sr->eigr[count] = lambda;
1122: sr->errest[count] = err;
1123: /* Explicit purification */
1124: PetscCall(BVGetColumn(eps->V,eps->perm[i],&w));
1125: if (eps->purify) {
1126: PetscCall(BVGetColumn(sr->V,count,&v));
1127: PetscCall(STApply(eps->st,w,v));
1128: PetscCall(BVRestoreColumn(sr->V,count,&v));
1129: } else PetscCall(BVInsertVec(sr->V,count,w));
1130: PetscCall(BVRestoreColumn(eps->V,eps->perm[i],&w));
1131: PetscCall(BVNormColumn(sr->V,count,NORM_2,&norm));
1132: PetscCall(BVScaleColumn(sr->V,count,1.0/norm));
1133: count++;
1134: }
1135: }
1136: sPres->neigs = count - sr->indexEig;
1137: sr->indexEig = count;
1138: /* Global ordering array updating */
1139: PetscCall(sortRealEigenvalues(sr->eigr,sr->perm,count,PETSC_TRUE,sr->dir));
1140: PetscFunctionReturn(PETSC_SUCCESS);
1141: }
1143: static PetscErrorCode EPSLookForDeflation(EPS eps)
1144: {
1145: PetscReal val;
1146: PetscInt i,count0=0,count1=0;
1147: EPS_shift sPres;
1148: PetscInt ini,fin,k,idx0,idx1;
1149: EPS_SR sr;
1150: Vec v;
1151: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1153: PetscFunctionBegin;
1154: sr = ctx->sr;
1155: sPres = sr->sPres;
1157: if (sPres->neighb[0]) ini = (sr->dir)*(sPres->neighb[0]->inertia - sr->inertia0);
1158: else ini = 0;
1159: fin = sr->indexEig;
1160: /* Selection of ends for searching new values */
1161: if (!sPres->neighb[0]) sPres->ext[0] = sr->int0;/* First shift */
1162: else sPres->ext[0] = sPres->neighb[0]->value;
1163: if (!sPres->neighb[1]) {
1164: if (sr->hasEnd) sPres->ext[1] = sr->int1;
1165: else sPres->ext[1] = (sr->dir > 0)?PETSC_MAX_REAL:PETSC_MIN_REAL;
1166: } else sPres->ext[1] = sPres->neighb[1]->value;
1167: /* Selection of values between right and left ends */
1168: for (i=ini;i<fin;i++) {
1169: val=PetscRealPart(sr->eigr[sr->perm[i]]);
1170: /* Values to the right of left shift */
1171: if (sr->dir*(val - sPres->ext[1]) < 0) {
1172: if (sr->dir*(val - sPres->value) < 0) count0++;
1173: else count1++;
1174: } else break;
1175: }
1176: /* The number of values on each side are found */
1177: if (sPres->neighb[0]) {
1178: sPres->nsch[0] = (sr->dir)*(sPres->inertia - sPres->neighb[0]->inertia)-count0;
1179: PetscCheck(sPres->nsch[0]>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia%s",ctx->detect?"":", consider using EPSKrylovSchurSetDetectZeros()");
1180: } else sPres->nsch[0] = 0;
1182: if (sPres->neighb[1]) {
1183: sPres->nsch[1] = (sr->dir)*(sPres->neighb[1]->inertia - sPres->inertia) - count1;
1184: PetscCheck(sPres->nsch[1]>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Mismatch between number of values found and information from inertia%s",ctx->detect?"":", consider using EPSKrylovSchurSetDetectZeros()");
1185: } else sPres->nsch[1] = (sr->dir)*(sr->inertia1 - sPres->inertia);
1187: /* Completing vector of indexes for deflation */
1188: idx0 = ini;
1189: idx1 = ini+count0+count1;
1190: k=0;
1191: for (i=idx0;i<idx1;i++) sr->idxDef[k++]=sr->perm[i];
1192: PetscCall(BVDuplicateResize(eps->V,k+eps->ncv+1,&sr->Vnext));
1193: PetscCall(BVSetNumConstraints(sr->Vnext,k));
1194: for (i=0;i<k;i++) {
1195: PetscCall(BVGetColumn(sr->Vnext,-i-1,&v));
1196: PetscCall(BVCopyVec(sr->V,sr->idxDef[i],v));
1197: PetscCall(BVRestoreColumn(sr->Vnext,-i-1,&v));
1198: }
1200: /* For rational Krylov */
1201: if (sr->nS>0 && (sr->sPrev == sr->sPres->neighb[0] || sr->sPrev == sr->sPres->neighb[1])) PetscCall(EPSPrepareRational(eps));
1202: eps->nconv = 0;
1203: /* Get rid of temporary Vnext */
1204: PetscCall(BVDestroy(&eps->V));
1205: eps->V = sr->Vnext;
1206: sr->Vnext = NULL;
1207: PetscFunctionReturn(PETSC_SUCCESS);
1208: }
1210: PetscErrorCode EPSSolve_KrylovSchur_Slice(EPS eps)
1211: {
1212: PetscInt i,lds,ti;
1213: PetscReal newS;
1214: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1215: EPS_SR sr=ctx->sr;
1216: Mat A,B=NULL;
1217: PetscObjectState Astate,Bstate=0;
1218: PetscObjectId Aid,Bid=0;
1220: PetscFunctionBegin;
1221: PetscCall(PetscCitationsRegister(citation,&cited));
1222: if (ctx->global) {
1223: PetscCall(EPSSolve_KrylovSchur_Slice(ctx->eps));
1224: ctx->eps->state = EPS_STATE_SOLVED;
1225: eps->reason = EPS_CONVERGED_TOL;
1226: if (ctx->npart>1) {
1227: /* Gather solution from subsolvers */
1228: PetscCall(EPSSliceGatherSolution(eps));
1229: } else {
1230: eps->nconv = sr->numEigs;
1231: eps->its = ctx->eps->its;
1232: PetscCall(PetscFree(ctx->inertias));
1233: PetscCall(PetscFree(ctx->shifts));
1234: PetscCall(EPSSliceGetInertias(ctx->eps,&ctx->nshifts,&ctx->shifts,&ctx->inertias));
1235: }
1236: } else {
1237: if (ctx->npart==1) {
1238: sr->eigr = ctx->eps->eigr;
1239: sr->eigi = ctx->eps->eigi;
1240: sr->perm = ctx->eps->perm;
1241: sr->errest = ctx->eps->errest;
1242: sr->V = ctx->eps->V;
1243: }
1244: /* Check that the user did not modify subcomm matrices */
1245: PetscCall(EPSGetOperators(eps,&A,&B));
1246: PetscCall(PetscObjectStateGet((PetscObject)A,&Astate));
1247: PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1248: if (B) {
1249: PetscCall(PetscObjectStateGet((PetscObject)B,&Bstate));
1250: PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1251: }
1252: PetscCheck(Astate==ctx->Astate && (!B || Bstate==ctx->Bstate) && Aid==ctx->Aid && (!B || Bid==ctx->Bid),PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Subcomm matrices have been modified by user");
1253: /* Only with eigenvalues present in the interval ...*/
1254: if (sr->numEigs==0) {
1255: eps->reason = EPS_CONVERGED_TOL;
1256: PetscFunctionReturn(PETSC_SUCCESS);
1257: }
1258: /* Array of pending shifts */
1259: sr->maxPend = 100; /* Initial size */
1260: sr->nPend = 0;
1261: PetscCall(PetscMalloc1(sr->maxPend,&sr->pending));
1262: PetscCall(EPSCreateShift(eps,sr->int0,NULL,NULL));
1263: /* extract first shift */
1264: sr->sPrev = NULL;
1265: sr->sPres = sr->pending[--sr->nPend];
1266: sr->sPres->inertia = sr->inertia0;
1267: eps->target = sr->sPres->value;
1268: sr->s0 = sr->sPres;
1269: sr->indexEig = 0;
1270: /* Memory reservation for auxiliary variables */
1271: lds = PetscMin(eps->mpd,eps->ncv);
1272: PetscCall(PetscCalloc1(lds*lds,&sr->S));
1273: PetscCall(PetscMalloc1(eps->ncv,&sr->back));
1274: for (i=0;i<sr->numEigs;i++) {
1275: sr->eigr[i] = 0.0;
1276: sr->eigi[i] = 0.0;
1277: sr->errest[i] = 0.0;
1278: sr->perm[i] = i;
1279: }
1280: /* Vectors for deflation */
1281: PetscCall(PetscMalloc1(sr->numEigs,&sr->idxDef));
1282: sr->indexEig = 0;
1283: /* Main loop */
1284: while (sr->sPres) {
1285: /* Search for deflation */
1286: PetscCall(EPSLookForDeflation(eps));
1287: /* KrylovSchur */
1288: PetscCall(EPSKrylovSchur_Slice(eps));
1290: PetscCall(EPSStoreEigenpairs(eps));
1291: /* Select new shift */
1292: if (!sr->sPres->comp[1]) {
1293: PetscCall(EPSGetNewShiftValue(eps,1,&newS));
1294: PetscCall(EPSCreateShift(eps,newS,sr->sPres,sr->sPres->neighb[1]));
1295: }
1296: if (!sr->sPres->comp[0]) {
1297: /* Completing earlier interval */
1298: PetscCall(EPSGetNewShiftValue(eps,0,&newS));
1299: PetscCall(EPSCreateShift(eps,newS,sr->sPres->neighb[0],sr->sPres));
1300: }
1301: /* Preparing for a new search of values */
1302: PetscCall(EPSExtractShift(eps));
1303: }
1305: /* Updating eps values prior to exit */
1306: PetscCall(PetscFree(sr->S));
1307: PetscCall(PetscFree(sr->idxDef));
1308: PetscCall(PetscFree(sr->pending));
1309: PetscCall(PetscFree(sr->back));
1310: PetscCall(BVDuplicateResize(eps->V,eps->ncv+1,&sr->Vnext));
1311: PetscCall(BVSetNumConstraints(sr->Vnext,0));
1312: PetscCall(BVDestroy(&eps->V));
1313: eps->V = sr->Vnext;
1314: eps->nconv = sr->indexEig;
1315: eps->reason = EPS_CONVERGED_TOL;
1316: eps->its = sr->itsKs;
1317: eps->nds = 0;
1318: if (sr->dir<0) {
1319: for (i=0;i<eps->nconv/2;i++) {
1320: ti = sr->perm[i]; sr->perm[i] = sr->perm[eps->nconv-1-i]; sr->perm[eps->nconv-1-i] = ti;
1321: }
1322: }
1323: }
1324: PetscFunctionReturn(PETSC_SUCCESS);
1325: }