Actual source code: ks-slice.c
slepc-3.22.2 2024-12-02
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_DETERMINE) 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_DETERMINE) {
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_DETERMINE) 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) {
345: PetscCall(PetscMPIIntCast((sr->dir>0)?0:ctx->npart-1,&aux));
346: PetscCallMPI(MPI_Bcast(&sr->inertia0,1,MPIU_INT,aux,ctx->commrank));
347: }
348: PetscCallMPI(MPI_Bcast(&sr->inertia0,1,MPIU_INT,0,child));
349: PetscCall(PetscFree(ctx->nconv_loc));
350: PetscCall(PetscMalloc1(ctx->npart,&ctx->nconv_loc));
351: PetscCallMPI(MPI_Comm_size(((PetscObject)eps)->comm,&nproc));
352: if (sr->dir<0) off = 1;
353: if (nproc%ctx->npart==0) { /* subcommunicators with the same size */
354: PetscCall(PetscMPIIntCast(sr_loc->numEigs,&aux));
355: PetscCallMPI(MPI_Allgather(&aux,1,MPI_INT,ctx->nconv_loc,1,MPI_INT,ctx->commrank));
356: 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));
357: } else {
358: PetscCallMPI(MPI_Comm_rank(child,&rank));
359: if (!rank) {
360: PetscCall(PetscMPIIntCast(sr_loc->numEigs,&aux));
361: PetscCallMPI(MPI_Allgather(&aux,1,MPI_INT,ctx->nconv_loc,1,MPI_INT,ctx->commrank));
362: 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));
363: }
364: PetscCall(PetscMPIIntCast(ctx->npart,&aux));
365: PetscCallMPI(MPI_Bcast(ctx->nconv_loc,aux,MPI_INT,0,child));
366: PetscCallMPI(MPI_Bcast(ctx->subintervals+off,aux,MPIU_REAL,0,child));
367: }
368: nEigs = 0;
369: for (i=0;i<ctx->npart;i++) nEigs += ctx->nconv_loc[i];
370: } else {
371: nEigs = sr_loc->numEigs;
372: sr->inertia0 = sr_loc->inertia0;
373: sr->dir = sr_loc->dir;
374: }
375: sr->inertia1 = sr->inertia0+sr->dir*nEigs;
376: sr->numEigs = nEigs;
377: eps->nev = nEigs;
378: eps->ncv = nEigs;
379: eps->mpd = nEigs;
380: } else {
381: ctx_glob = (EPS_KRYLOVSCHUR*)ctx->eps->data;
382: sr_glob = ctx_glob->sr;
383: if (ctx->npart>1) {
384: sr->dir = sr_glob->dir;
385: sr->int0 = (sr->dir==1)?eps->inta:eps->intb;
386: sr->int1 = (sr->dir==1)?eps->intb:eps->inta;
387: if ((sr->dir>0&&ctx->subc->color==ctx->npart-1)||(sr->dir<0&&ctx->subc->color==0)) sr->hasEnd = sr_glob->hasEnd;
388: else sr->hasEnd = PETSC_TRUE;
389: }
390: /* sets first shift */
391: PetscCall(STSetShift(eps->st,(sr->int0==0.0)?10.0/PETSC_MAX_REAL:sr->int0));
392: PetscCall(STSetUp(eps->st));
394: /* compute inertia0 */
395: PetscCall(EPSSliceGetInertia(eps,sr->int0,&sr->inertia0,ctx->detect?&zeros:NULL));
396: /* undocumented option to control what to do when an eigenvalue is found:
397: - error out if it's the endpoint of the user-provided interval (or sub-interval)
398: - if it's an endpoint computed internally:
399: + if hiteig=0 error out
400: + else if hiteig=1 the subgroup that hit the eigenvalue does nothing
401: + otherwise the subgroup that hit the eigenvalue perturbs the shift and recomputes inertia
402: */
403: PetscCall(PetscOptionsGetInt(NULL,NULL,"-eps_krylovschur_hiteigenvalue",&hiteig,NULL));
404: if (zeros) { /* error in factorization */
405: 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");
406: 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");
407: if (hiteig==1) { /* idle subgroup */
408: sr->inertia0 = -1;
409: } else { /* perturb shift */
410: sr->int0 *= (1.0+SLICE_PTOL);
411: PetscCall(EPSSliceGetInertia(eps,sr->int0,&sr->inertia0,&zeros));
412: PetscCheck(zeros==0,((PetscObject)eps)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)sr->int1);
413: }
414: }
415: if (ctx->npart>1) {
416: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
417: /* inertia1 is received from neighbour */
418: PetscCallMPI(MPI_Comm_rank(child,&rank));
419: if (!rank) {
420: if (sr->inertia0!=-1 && ((sr->dir>0 && ctx->subc->color>0) || (sr->dir<0 && ctx->subc->color<ctx->npart-1))) { /* send inertia0 to neighbour0 */
421: PetscCall(PetscMPIIntCast(ctx->subc->color-sr->dir,&aux));
422: PetscCallMPI(MPI_Isend(&sr->inertia0,1,MPIU_INT,aux,0,ctx->commrank,&req));
423: PetscCallMPI(MPI_Isend(&sr->int0,1,MPIU_REAL,aux,0,ctx->commrank,&req));
424: }
425: if ((sr->dir>0 && ctx->subc->color<ctx->npart-1)|| (sr->dir<0 && ctx->subc->color>0)) { /* receive inertia1 from neighbour1 */
426: PetscCall(PetscMPIIntCast(ctx->subc->color+sr->dir,&aux));
427: PetscCallMPI(MPI_Recv(&sr->inertia1,1,MPIU_INT,aux,0,ctx->commrank,MPI_STATUS_IGNORE));
428: PetscCallMPI(MPI_Recv(&sr->int1,1,MPIU_REAL,aux,0,ctx->commrank,MPI_STATUS_IGNORE));
429: }
430: if (sr->inertia0==-1 && !(sr->dir>0 && ctx->subc->color==ctx->npart-1) && !(sr->dir<0 && ctx->subc->color==0)) {
431: sr->inertia0 = sr->inertia1; sr->int0 = sr->int1;
432: PetscCall(PetscMPIIntCast(ctx->subc->color-sr->dir,&aux));
433: PetscCallMPI(MPI_Isend(&sr->inertia0,1,MPIU_INT,aux,0,ctx->commrank,&req));
434: PetscCallMPI(MPI_Isend(&sr->int0,1,MPIU_REAL,aux,0,ctx->commrank,&req));
435: }
436: }
437: if ((sr->dir>0 && ctx->subc->color<ctx->npart-1)||(sr->dir<0 && ctx->subc->color>0)) {
438: PetscCallMPI(MPI_Bcast(&sr->inertia1,1,MPIU_INT,0,child));
439: PetscCallMPI(MPI_Bcast(&sr->int1,1,MPIU_REAL,0,child));
440: } else sr_glob->inertia1 = sr->inertia1;
441: }
443: /* last process in eps comm computes inertia1 */
444: if (ctx->npart==1 || ((sr->dir>0 && ctx->subc->color==ctx->npart-1) || (sr->dir<0 && ctx->subc->color==0))) {
445: PetscCall(EPSSliceGetInertia(eps,sr->int1,&sr->inertia1,ctx->detect?&zeros:NULL));
446: PetscCheck(zeros==0,((PetscObject)eps)->comm,PETSC_ERR_USER,"Found singular matrix for the transformed problem in an interval endpoint defined by user");
447: if (!rank && sr->inertia0==-1) {
448: sr->inertia0 = sr->inertia1; sr->int0 = sr->int1;
449: PetscCall(PetscMPIIntCast(ctx->subc->color-sr->dir,&aux));
450: PetscCallMPI(MPI_Isend(&sr->inertia0,1,MPIU_INT,aux,0,ctx->commrank,&req));
451: PetscCallMPI(MPI_Isend(&sr->int0,1,MPIU_REAL,aux,0,ctx->commrank,&req));
452: }
453: if (sr->hasEnd) {
454: sr->dir = -sr->dir; r = sr->int0; sr->int0 = sr->int1; sr->int1 = r;
455: i = sr->inertia0; sr->inertia0 = sr->inertia1; sr->inertia1 = i;
456: }
457: }
459: /* number of eigenvalues in interval */
460: sr->numEigs = (sr->dir)*(sr->inertia1 - sr->inertia0);
461: if (ctx->npart>1) {
462: /* memory allocate for subinterval eigenpairs */
463: PetscCall(EPSSliceAllocateSolution(eps,1));
464: }
465: dssz = eps->ncv+1;
466: PetscCall(DSGetParallel(ctx->eps->ds,&ptype));
467: PetscCall(DSSetParallel(eps->ds,ptype));
468: PetscCall(DSGetMethod(ctx->eps->ds,&method));
469: PetscCall(DSSetMethod(eps->ds,method));
470: }
471: PetscCall(DSSetType(eps->ds,DSHEP));
472: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
473: PetscCall(DSAllocate(eps->ds,dssz));
474: /* keep state of subcomm matrices to check that the user does not modify them */
475: PetscCall(EPSGetOperators(eps,&A,&B));
476: PetscCall(MatGetState(A,&ctx->Astate));
477: PetscCall(PetscObjectGetId((PetscObject)A,&ctx->Aid));
478: if (B) {
479: PetscCall(MatGetState(B,&ctx->Bstate));
480: PetscCall(PetscObjectGetId((PetscObject)B,&ctx->Bid));
481: } else {
482: ctx->Bstate=0;
483: ctx->Bid=0;
484: }
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: static PetscErrorCode EPSSliceGatherEigenVectors(EPS eps)
489: {
490: Vec v,vg,v_loc;
491: IS is1,is2;
492: VecScatter vec_sc;
493: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
494: PetscInt nloc,m0,n0,i,si,idx,*idx1,*idx2,j;
495: PetscScalar *array;
496: EPS_SR sr_loc;
497: BV V_loc;
499: PetscFunctionBegin;
500: sr_loc = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
501: V_loc = sr_loc->V;
503: /* Gather parallel eigenvectors */
504: PetscCall(BVGetColumn(eps->V,0,&v));
505: PetscCall(VecGetOwnershipRange(v,&n0,&m0));
506: PetscCall(BVRestoreColumn(eps->V,0,&v));
507: PetscCall(BVGetColumn(ctx->eps->V,0,&v));
508: PetscCall(VecGetLocalSize(v,&nloc));
509: PetscCall(BVRestoreColumn(ctx->eps->V,0,&v));
510: PetscCall(PetscMalloc2(m0-n0,&idx1,m0-n0,&idx2));
511: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)eps),nloc,PETSC_DECIDE,&vg));
512: idx = -1;
513: for (si=0;si<ctx->npart;si++) {
514: j = 0;
515: for (i=n0;i<m0;i++) {
516: idx1[j] = i;
517: idx2[j++] = i+eps->n*si;
518: }
519: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)eps),(m0-n0),idx1,PETSC_COPY_VALUES,&is1));
520: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)eps),(m0-n0),idx2,PETSC_COPY_VALUES,&is2));
521: PetscCall(BVGetColumn(eps->V,0,&v));
522: PetscCall(VecScatterCreate(v,is1,vg,is2,&vec_sc));
523: PetscCall(BVRestoreColumn(eps->V,0,&v));
524: PetscCall(ISDestroy(&is1));
525: PetscCall(ISDestroy(&is2));
526: for (i=0;i<ctx->nconv_loc[si];i++) {
527: PetscCall(BVGetColumn(eps->V,++idx,&v));
528: if (ctx->subc->color==si) {
529: PetscCall(BVGetColumn(V_loc,i,&v_loc));
530: PetscCall(VecGetArray(v_loc,&array));
531: PetscCall(VecPlaceArray(vg,array));
532: }
533: PetscCall(VecScatterBegin(vec_sc,vg,v,INSERT_VALUES,SCATTER_REVERSE));
534: PetscCall(VecScatterEnd(vec_sc,vg,v,INSERT_VALUES,SCATTER_REVERSE));
535: if (ctx->subc->color==si) {
536: PetscCall(VecResetArray(vg));
537: PetscCall(VecRestoreArray(v_loc,&array));
538: PetscCall(BVRestoreColumn(V_loc,i,&v_loc));
539: }
540: PetscCall(BVRestoreColumn(eps->V,idx,&v));
541: }
542: PetscCall(VecScatterDestroy(&vec_sc));
543: }
544: PetscCall(PetscFree2(idx1,idx2));
545: PetscCall(VecDestroy(&vg));
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: /*
550: EPSComputeVectors_Slice - Recover Eigenvectors from subcomunicators
551: */
552: PetscErrorCode EPSComputeVectors_Slice(EPS eps)
553: {
554: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
556: PetscFunctionBegin;
557: if (ctx->global && ctx->npart>1) {
558: PetscCall(EPSComputeVectors(ctx->eps));
559: PetscCall(EPSSliceGatherEigenVectors(eps));
560: }
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: static PetscErrorCode EPSSliceGetInertias(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
565: {
566: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
567: PetscInt i=0,j,tmpi;
568: PetscReal v,tmpr;
569: EPS_shift s;
571: PetscFunctionBegin;
572: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
573: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
574: if (!ctx->sr->s0) { /* EPSSolve not called yet */
575: *n = 2;
576: } else {
577: *n = 1;
578: s = ctx->sr->s0;
579: while (s) {
580: (*n)++;
581: s = s->neighb[1];
582: }
583: }
584: PetscCall(PetscMalloc1(*n,shifts));
585: PetscCall(PetscMalloc1(*n,inertias));
586: if (!ctx->sr->s0) { /* EPSSolve not called yet */
587: (*shifts)[0] = ctx->sr->int0;
588: (*shifts)[1] = ctx->sr->int1;
589: (*inertias)[0] = ctx->sr->inertia0;
590: (*inertias)[1] = ctx->sr->inertia1;
591: } else {
592: s = ctx->sr->s0;
593: while (s) {
594: (*shifts)[i] = s->value;
595: (*inertias)[i++] = s->inertia;
596: s = s->neighb[1];
597: }
598: (*shifts)[i] = ctx->sr->int1;
599: (*inertias)[i] = ctx->sr->inertia1;
600: }
601: /* remove possible duplicate in last position */
602: if ((*shifts)[(*n)-1]==(*shifts)[(*n)-2]) (*n)--;
603: /* sort result */
604: for (i=0;i<*n;i++) {
605: v = (*shifts)[i];
606: for (j=i+1;j<*n;j++) {
607: if (v > (*shifts)[j]) {
608: SlepcSwap((*shifts)[i],(*shifts)[j],tmpr);
609: SlepcSwap((*inertias)[i],(*inertias)[j],tmpi);
610: v = (*shifts)[i];
611: }
612: }
613: }
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: static PetscErrorCode EPSSliceGatherSolution(EPS eps)
618: {
619: PetscMPIInt rank,nproc,*disp,*ns_loc,aux;
620: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
621: PetscInt i,idx,j,*perm_loc,off=0,*inertias_loc,ns;
622: PetscScalar *eigr_loc;
623: EPS_SR sr_loc;
624: PetscReal *shifts_loc;
625: MPI_Comm child;
627: PetscFunctionBegin;
628: eps->nconv = 0;
629: for (i=0;i<ctx->npart;i++) eps->nconv += ctx->nconv_loc[i];
630: sr_loc = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
632: /* Gather the shifts used and the inertias computed */
633: PetscCall(EPSSliceGetInertias(ctx->eps,&ns,&shifts_loc,&inertias_loc));
634: if (ctx->sr->dir>0 && shifts_loc[ns-1]==sr_loc->int1 && ctx->subc->color<ctx->npart-1) ns--;
635: if (ctx->sr->dir<0 && shifts_loc[ns-1]==sr_loc->int0 && ctx->subc->color>0) {
636: ns--;
637: for (i=0;i<ns;i++) {
638: inertias_loc[i] = inertias_loc[i+1];
639: shifts_loc[i] = shifts_loc[i+1];
640: }
641: }
642: PetscCall(PetscMalloc1(ctx->npart,&ns_loc));
643: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
644: PetscCallMPI(MPI_Comm_rank(child,&rank));
645: PetscCall(PetscMPIIntCast(ns,&aux));
646: if (!rank) PetscCallMPI(MPI_Allgather(&aux,1,MPI_INT,ns_loc,1,MPI_INT,ctx->commrank));
647: PetscCall(PetscMPIIntCast(ctx->npart,&aux));
648: PetscCallMPI(MPI_Bcast(ns_loc,aux,MPI_INT,0,child));
649: ctx->nshifts = 0;
650: for (i=0;i<ctx->npart;i++) ctx->nshifts += ns_loc[i];
651: PetscCall(PetscFree(ctx->inertias));
652: PetscCall(PetscFree(ctx->shifts));
653: PetscCall(PetscMalloc1(ctx->nshifts,&ctx->inertias));
654: PetscCall(PetscMalloc1(ctx->nshifts,&ctx->shifts));
656: /* Gather eigenvalues (same ranks have fully set of eigenvalues)*/
657: eigr_loc = sr_loc->eigr;
658: perm_loc = sr_loc->perm;
659: PetscCallMPI(MPI_Comm_size(((PetscObject)eps)->comm,&nproc));
660: PetscCall(PetscMalloc1(ctx->npart,&disp));
661: disp[0] = 0;
662: for (i=1;i<ctx->npart;i++) disp[i] = disp[i-1]+ctx->nconv_loc[i-1];
663: if (nproc%ctx->npart==0) { /* subcommunicators with the same size */
664: PetscCall(PetscMPIIntCast(sr_loc->numEigs,&aux));
665: PetscCallMPI(MPI_Allgatherv(eigr_loc,aux,MPIU_SCALAR,eps->eigr,ctx->nconv_loc,disp,MPIU_SCALAR,ctx->commrank)); /* eigenvalues */
666: PetscCallMPI(MPI_Allgatherv(perm_loc,aux,MPIU_INT,eps->perm,ctx->nconv_loc,disp,MPIU_INT,ctx->commrank)); /* perm */
667: for (i=1;i<ctx->npart;i++) disp[i] = disp[i-1]+ns_loc[i-1];
668: PetscCall(PetscMPIIntCast(ns,&aux));
669: PetscCallMPI(MPI_Allgatherv(shifts_loc,aux,MPIU_REAL,ctx->shifts,ns_loc,disp,MPIU_REAL,ctx->commrank)); /* shifts */
670: PetscCallMPI(MPI_Allgatherv(inertias_loc,aux,MPIU_INT,ctx->inertias,ns_loc,disp,MPIU_INT,ctx->commrank)); /* inertias */
671: PetscCallMPI(MPIU_Allreduce(&sr_loc->itsKs,&eps->its,1,MPIU_INT,MPI_SUM,ctx->commrank));
672: } else { /* subcommunicators with different size */
673: if (!rank) {
674: PetscCall(PetscMPIIntCast(sr_loc->numEigs,&aux));
675: PetscCallMPI(MPI_Allgatherv(eigr_loc,aux,MPIU_SCALAR,eps->eigr,ctx->nconv_loc,disp,MPIU_SCALAR,ctx->commrank)); /* eigenvalues */
676: PetscCallMPI(MPI_Allgatherv(perm_loc,aux,MPIU_INT,eps->perm,ctx->nconv_loc,disp,MPIU_INT,ctx->commrank)); /* perm */
677: for (i=1;i<ctx->npart;i++) disp[i] = disp[i-1]+ns_loc[i-1];
678: PetscCall(PetscMPIIntCast(ns,&aux));
679: PetscCallMPI(MPI_Allgatherv(shifts_loc,aux,MPIU_REAL,ctx->shifts,ns_loc,disp,MPIU_REAL,ctx->commrank)); /* shifts */
680: PetscCallMPI(MPI_Allgatherv(inertias_loc,aux,MPIU_INT,ctx->inertias,ns_loc,disp,MPIU_INT,ctx->commrank)); /* inertias */
681: PetscCallMPI(MPIU_Allreduce(&sr_loc->itsKs,&eps->its,1,MPIU_INT,MPI_SUM,ctx->commrank));
682: }
683: PetscCall(PetscMPIIntCast(eps->nconv,&aux));
684: PetscCallMPI(MPI_Bcast(eps->eigr,aux,MPIU_SCALAR,0,child));
685: PetscCallMPI(MPI_Bcast(eps->perm,aux,MPIU_INT,0,child));
686: PetscCall(PetscMPIIntCast(ctx->nshifts,&aux));
687: PetscCallMPI(MPI_Bcast(ctx->shifts,aux,MPIU_REAL,0,child));
688: PetscCallMPI(MPI_Bcast(ctx->inertias,aux,MPIU_INT,0,child));
689: PetscCallMPI(MPI_Bcast(&eps->its,1,MPIU_INT,0,child));
690: }
691: /* Update global array eps->perm */
692: idx = ctx->nconv_loc[0];
693: for (i=1;i<ctx->npart;i++) {
694: off += ctx->nconv_loc[i-1];
695: for (j=0;j<ctx->nconv_loc[i];j++) eps->perm[idx++] += off;
696: }
698: /* Gather parallel eigenvectors */
699: PetscCall(PetscFree(ns_loc));
700: PetscCall(PetscFree(disp));
701: PetscCall(PetscFree(shifts_loc));
702: PetscCall(PetscFree(inertias_loc));
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: /*
707: Fills the fields of a shift structure
708: */
709: static PetscErrorCode EPSCreateShift(EPS eps,PetscReal val,EPS_shift neighb0,EPS_shift neighb1)
710: {
711: EPS_shift s,*pending2;
712: PetscInt i;
713: EPS_SR sr;
714: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
716: PetscFunctionBegin;
717: sr = ctx->sr;
718: if ((neighb0 && val==neighb0->value) || (neighb1 && val==neighb1->value)) {
719: sr->nPend++;
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
722: PetscCall(PetscNew(&s));
723: s->value = val;
724: s->neighb[0] = neighb0;
725: if (neighb0) neighb0->neighb[1] = s;
726: s->neighb[1] = neighb1;
727: if (neighb1) neighb1->neighb[0] = s;
728: s->comp[0] = PETSC_FALSE;
729: s->comp[1] = PETSC_FALSE;
730: s->index = -1;
731: s->neigs = 0;
732: s->nconv[0] = s->nconv[1] = 0;
733: s->nsch[0] = s->nsch[1]=0;
734: /* Inserts in the stack of pending shifts */
735: /* If needed, the array is resized */
736: if (sr->nPend >= sr->maxPend) {
737: sr->maxPend *= 2;
738: PetscCall(PetscMalloc1(sr->maxPend,&pending2));
739: for (i=0;i<sr->nPend;i++) pending2[i] = sr->pending[i];
740: PetscCall(PetscFree(sr->pending));
741: sr->pending = pending2;
742: }
743: sr->pending[sr->nPend++]=s;
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: /* Prepare for Rational Krylov update */
748: static PetscErrorCode EPSPrepareRational(EPS eps)
749: {
750: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
751: PetscInt dir,i,k,ld,nv;
752: PetscScalar *A;
753: EPS_SR sr = ctx->sr;
754: Vec v;
756: PetscFunctionBegin;
757: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
758: dir = (sr->sPres->neighb[0] == sr->sPrev)?1:-1;
759: dir*=sr->dir;
760: k = 0;
761: for (i=0;i<sr->nS;i++) {
762: if (dir*PetscRealPart(sr->S[i])>0.0) {
763: sr->S[k] = sr->S[i];
764: sr->S[sr->nS+k] = sr->S[sr->nS+i];
765: PetscCall(BVGetColumn(sr->Vnext,k,&v));
766: PetscCall(BVCopyVec(eps->V,eps->nconv+i,v));
767: PetscCall(BVRestoreColumn(sr->Vnext,k,&v));
768: k++;
769: if (k>=sr->nS/2) break;
770: }
771: }
772: /* Copy to DS */
773: PetscCall(DSSetCompact(eps->ds,PETSC_FALSE)); /* make sure DS_MAT_A is allocated */
774: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&A));
775: PetscCall(PetscArrayzero(A,ld*ld));
776: for (i=0;i<k;i++) {
777: A[i*(1+ld)] = sr->S[i];
778: A[k+i*ld] = sr->S[sr->nS+i];
779: }
780: sr->nS = k;
781: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&A));
782: PetscCall(DSGetDimensions(eps->ds,&nv,NULL,NULL,NULL));
783: PetscCall(DSSetDimensions(eps->ds,nv,0,k));
784: /* Append u to V */
785: PetscCall(BVGetColumn(sr->Vnext,sr->nS,&v));
786: PetscCall(BVCopyVec(eps->V,sr->nv,v));
787: PetscCall(BVRestoreColumn(sr->Vnext,sr->nS,&v));
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: /* Provides next shift to be computed */
792: static PetscErrorCode EPSExtractShift(EPS eps)
793: {
794: PetscInt iner,zeros=0;
795: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
796: EPS_SR sr;
797: PetscReal newShift,diam,ptol;
798: EPS_shift sPres;
800: PetscFunctionBegin;
801: sr = ctx->sr;
802: if (sr->nPend > 0) {
803: if (sr->sPres==sr->pending[sr->nPend-1]) {
804: eps->reason = EPS_CONVERGED_ITERATING;
805: eps->its = 0;
806: sr->nPend--;
807: sr->sPres->rep = PETSC_TRUE;
808: PetscFunctionReturn(PETSC_SUCCESS);
809: }
810: sr->sPrev = sr->sPres;
811: sr->sPres = sr->pending[--sr->nPend];
812: sPres = sr->sPres;
813: PetscCall(EPSSliceGetInertia(eps,sPres->value,&iner,ctx->detect?&zeros:NULL));
814: if (zeros) {
815: diam = PetscMin(PetscAbsReal(sPres->neighb[0]->value-sPres->value),PetscAbsReal(sPres->neighb[1]->value-sPres->value));
816: ptol = PetscMin(SLICE_PTOL,diam/2);
817: newShift = sPres->value*(1.0+ptol);
818: if (sr->dir*(sPres->neighb[0] && newShift-sPres->neighb[0]->value) < 0) newShift = (sPres->value+sPres->neighb[0]->value)/2;
819: else if (sPres->neighb[1] && sr->dir*(sPres->neighb[1]->value-newShift) < 0) newShift = (sPres->value+sPres->neighb[1]->value)/2;
820: PetscCall(EPSSliceGetInertia(eps,newShift,&iner,&zeros));
821: PetscCheck(zeros==0,((PetscObject)eps)->comm,PETSC_ERR_CONV_FAILED,"Inertia computation fails in %g",(double)newShift);
822: sPres->value = newShift;
823: }
824: sr->sPres->inertia = iner;
825: eps->target = sr->sPres->value;
826: eps->reason = EPS_CONVERGED_ITERATING;
827: eps->its = 0;
828: } else sr->sPres = NULL;
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: /*
833: Symmetric KrylovSchur adapted to spectrum slicing:
834: Allows searching an specific amount of eigenvalues in the subintervals left and right.
835: Returns whether the search has succeeded
836: */
837: static PetscErrorCode EPSKrylovSchur_Slice(EPS eps)
838: {
839: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
840: PetscInt i,k,l,ld,nv,*iwork,j,count0,count1,iterCompl=0,n0,n1;
841: Mat U,Op,T;
842: PetscScalar *Q,*A;
843: PetscReal *a,*b,beta,lambda;
844: EPS_shift sPres;
845: PetscBool breakdown,complIterating,sch0,sch1;
846: EPS_SR sr = ctx->sr;
848: PetscFunctionBegin;
849: /* Spectrum slicing data */
850: sPres = sr->sPres;
851: complIterating =PETSC_FALSE;
852: sch1 = sch0 = PETSC_TRUE;
853: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
854: PetscCall(PetscMalloc1(2*ld,&iwork));
855: count0=0;count1=0; /* Found on both sides */
856: if (!sPres->rep && sr->nS > 0 && (sPres->neighb[0] == sr->sPrev || sPres->neighb[1] == sr->sPrev)) {
857: /* Rational Krylov */
858: PetscCall(DSTranslateRKS(eps->ds,sr->sPrev->value-sPres->value));
859: PetscCall(DSGetDimensions(eps->ds,NULL,NULL,&l,NULL));
860: PetscCall(DSSetDimensions(eps->ds,l+1,0,0));
861: PetscCall(BVSetActiveColumns(eps->V,0,l+1));
862: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
863: PetscCall(BVMultInPlace(eps->V,U,0,l+1));
864: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
865: } else {
866: /* Get the starting Lanczos vector */
867: PetscCall(EPSGetStartVector(eps,0,NULL));
868: l = 0;
869: }
870: /* Restart loop */
871: while (eps->reason == EPS_CONVERGED_ITERATING) {
872: eps->its++; sr->itsKs++;
873: /* Compute an nv-step Lanczos factorization */
874: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
875: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
876: PetscCall(DSGetMat(eps->ds,DS_MAT_T,&T));
877: PetscCall(STGetOperator(eps->st,&Op));
878: PetscCall(BVMatLanczos(eps->V,Op,T,eps->nconv+l,&nv,&beta,&breakdown));
879: PetscCall(STRestoreOperator(eps->st,&Op));
880: sr->nv = nv;
881: PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&T));
882: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
883: if (l==0) PetscCall(DSSetState(eps->ds,DS_STATE_INTERMEDIATE));
884: else PetscCall(DSSetState(eps->ds,DS_STATE_RAW));
885: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
887: /* Solve projected problem and compute residual norm estimates */
888: if (eps->its == 1 && l > 0) { /* After rational update, DS_MAT_A is available */
889: PetscCall(DSGetArray(eps->ds,DS_MAT_A,&A));
890: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
891: b = a + ld;
892: k = eps->nconv+l;
893: A[k*ld+k-1] = A[(k-1)*ld+k];
894: A[k*ld+k] = a[k];
895: for (j=k+1; j< nv; j++) {
896: A[j*ld+j] = a[j];
897: A[j*ld+j-1] = b[j-1] ;
898: A[(j-1)*ld+j] = b[j-1];
899: }
900: PetscCall(DSRestoreArray(eps->ds,DS_MAT_A,&A));
901: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
902: PetscCall(DSSolve(eps->ds,eps->eigr,NULL));
903: PetscCall(DSSort(eps->ds,eps->eigr,NULL,NULL,NULL,NULL));
904: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
905: } else { /* Restart */
906: PetscCall(DSSolve(eps->ds,eps->eigr,NULL));
907: PetscCall(DSSort(eps->ds,eps->eigr,NULL,NULL,NULL,NULL));
908: }
909: PetscCall(DSSynchronize(eps->ds,eps->eigr,NULL));
911: /* Residual */
912: PetscCall(EPSKrylovConvergence(eps,PETSC_TRUE,eps->nconv,nv-eps->nconv,beta,0.0,1.0,&k));
913: /* Checking values obtained for completing */
914: for (i=0;i<k;i++) {
915: sr->back[i]=eps->eigr[i];
916: }
917: PetscCall(STBackTransform(eps->st,k,sr->back,eps->eigi));
918: count0=count1=0;
919: for (i=0;i<k;i++) {
920: lambda = PetscRealPart(sr->back[i]);
921: if ((sr->dir*(sPres->value - lambda) > 0) && (sr->dir*(lambda - sPres->ext[0]) > 0)) count0++;
922: if ((sr->dir*(lambda - sPres->value) > 0) && (sr->dir*(sPres->ext[1] - lambda) > 0)) count1++;
923: }
924: if (k>eps->nev && eps->ncv-k<5) eps->reason = EPS_CONVERGED_TOL;
925: else {
926: /* Checks completion */
927: if ((!sch0||count0 >= sPres->nsch[0]) && (!sch1 ||count1 >= sPres->nsch[1])) {
928: eps->reason = EPS_CONVERGED_TOL;
929: } else {
930: if (!complIterating && eps->its >= eps->max_it) eps->reason = EPS_DIVERGED_ITS;
931: if (complIterating) {
932: if (--iterCompl <= 0) eps->reason = EPS_DIVERGED_ITS;
933: } else if (k >= eps->nev) {
934: n0 = sPres->nsch[0]-count0;
935: n1 = sPres->nsch[1]-count1;
936: if (sr->iterCompl>0 && ((n0>0 && n0<= sr->nMAXCompl)||(n1>0&&n1<=sr->nMAXCompl))) {
937: /* Iterating for completion*/
938: complIterating = PETSC_TRUE;
939: if (n0 >sr->nMAXCompl)sch0 = PETSC_FALSE;
940: if (n1 >sr->nMAXCompl)sch1 = PETSC_FALSE;
941: iterCompl = sr->iterCompl;
942: } else eps->reason = EPS_CONVERGED_TOL;
943: }
944: }
945: }
946: /* Update l */
947: if (eps->reason == EPS_CONVERGED_ITERATING) l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
948: else l = nv-k;
949: if (breakdown) l=0;
950: if (!ctx->lock && l>0 && eps->reason == EPS_CONVERGED_ITERATING) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
952: if (eps->reason == EPS_CONVERGED_ITERATING) {
953: if (breakdown) {
954: /* Start a new Lanczos factorization */
955: PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
956: PetscCall(EPSGetStartVector(eps,k,&breakdown));
957: if (breakdown) {
958: eps->reason = EPS_DIVERGED_BREAKDOWN;
959: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
960: }
961: } else {
962: /* Prepare the Rayleigh quotient for restart */
963: PetscCall(DSGetArrayReal(eps->ds,DS_MAT_T,&a));
964: PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Q));
965: b = a + ld;
966: for (i=k;i<k+l;i++) {
967: a[i] = PetscRealPart(eps->eigr[i]);
968: b[i] = PetscRealPart(Q[nv-1+i*ld]*beta);
969: }
970: PetscCall(DSRestoreArrayReal(eps->ds,DS_MAT_T,&a));
971: PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Q));
972: }
973: }
974: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
975: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
976: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
977: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
979: /* Normalize u and append it to V */
980: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(eps->V,nv,k+l));
981: eps->nconv = k;
982: if (eps->reason != EPS_CONVERGED_ITERATING) {
983: /* Store approximated values for next shift */
984: PetscCall(DSGetArray(eps->ds,DS_MAT_Q,&Q));
985: sr->nS = l;
986: for (i=0;i<l;i++) {
987: sr->S[i] = eps->eigr[i+k];/* Diagonal elements */
988: sr->S[i+l] = Q[nv-1+(i+k)*ld]*beta; /* Out of diagonal elements */
989: }
990: PetscCall(DSRestoreArray(eps->ds,DS_MAT_Q,&Q));
991: }
992: }
993: /* Check for completion */
994: for (i=0;i< eps->nconv; i++) {
995: if (sr->dir*PetscRealPart(eps->eigr[i])>0) sPres->nconv[1]++;
996: else sPres->nconv[0]++;
997: }
998: sPres->comp[0] = PetscNot(count0 < sPres->nsch[0]);
999: sPres->comp[1] = PetscNot(count1 < sPres->nsch[1]);
1000: 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]?"*":""));
1001: 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()");
1002: PetscCall(PetscFree(iwork));
1003: PetscFunctionReturn(PETSC_SUCCESS);
1004: }
1006: /*
1007: Obtains value of subsequent shift
1008: */
1009: static PetscErrorCode EPSGetNewShiftValue(EPS eps,PetscInt side,PetscReal *newS)
1010: {
1011: PetscReal lambda,d_prev;
1012: PetscInt i,idxP;
1013: EPS_SR sr;
1014: EPS_shift sPres,s;
1015: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1017: PetscFunctionBegin;
1018: sr = ctx->sr;
1019: sPres = sr->sPres;
1020: if (sPres->neighb[side]) {
1021: /* Completing a previous interval */
1022: *newS = (sPres->value + sPres->neighb[side]->value)/2;
1023: if (PetscAbsReal(sPres->value - *newS)/PetscAbsReal(sPres->value)<=100*PETSC_SQRT_MACHINE_EPSILON) *newS = sPres->value;
1024: } else { /* (Only for side=1). Creating a new interval. */
1025: if (sPres->neigs==0) {/* No value has been accepted*/
1026: if (sPres->neighb[0]) {
1027: /* Multiplying by 10 the previous distance */
1028: *newS = sPres->value + 10*sr->dir*PetscAbsReal(sPres->value - sPres->neighb[0]->value);
1029: sr->nleap++;
1030: /* Stops when the interval is open and no values are found in the last 5 shifts (there might be infinite eigenvalues) */
1031: PetscCheck(sr->hasEnd || sr->nleap<=5,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unable to compute the wanted eigenvalues with open interval");
1032: } else { /* First shift */
1033: PetscCheck(eps->nconv!=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"First shift renders no information");
1034: /* Unaccepted values give information for next shift */
1035: idxP=0;/* Number of values left from shift */
1036: for (i=0;i<eps->nconv;i++) {
1037: lambda = PetscRealPart(eps->eigr[i]);
1038: if (sr->dir*(lambda - sPres->value) <0) idxP++;
1039: else break;
1040: }
1041: /* Avoiding subtraction of eigenvalues (might be the same).*/
1042: if (idxP>0) {
1043: d_prev = PetscAbsReal(sPres->value - PetscRealPart(eps->eigr[0]))/(idxP+0.3);
1044: } else {
1045: d_prev = PetscAbsReal(sPres->value - PetscRealPart(eps->eigr[eps->nconv-1]))/(eps->nconv+0.3);
1046: }
1047: *newS = sPres->value + (sr->dir*d_prev*eps->nev)/2;
1048: }
1049: } else { /* Accepted values found */
1050: sr->nleap = 0;
1051: /* Average distance of values in previous subinterval */
1052: s = sPres->neighb[0];
1053: while (s && PetscAbs(s->inertia - sPres->inertia)==0) {
1054: s = s->neighb[0];/* Looking for previous shifts with eigenvalues within */
1055: }
1056: if (s) {
1057: d_prev = PetscAbsReal((sPres->value - s->value)/(sPres->inertia - s->inertia));
1058: } else { /* First shift. Average distance obtained with values in this shift */
1059: /* first shift might be too far from first wanted eigenvalue (no values found outside the interval)*/
1060: 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)) {
1061: d_prev = PetscAbsReal((PetscRealPart(sr->eigr[sr->indexEig-1]) - PetscRealPart(sr->eigr[0])))/(sPres->neigs+0.3);
1062: } else {
1063: d_prev = PetscAbsReal(PetscRealPart(sr->eigr[sr->indexEig-1]) - sPres->value)/(sPres->neigs+0.3);
1064: }
1065: }
1066: /* Average distance is used for next shift by adding it to value on the right or to shift */
1067: if (sr->dir*(PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1]) - sPres->value)>0) {
1068: *newS = PetscRealPart(sr->eigr[sPres->index + sPres->neigs -1])+ (sr->dir*d_prev*eps->nev)/2;
1069: } else { /* Last accepted value is on the left of shift. Adding to shift */
1070: *newS = sPres->value + (sr->dir*d_prev*eps->nev)/2;
1071: }
1072: }
1073: /* End of interval can not be surpassed */
1074: if (sr->dir*(sr->int1 - *newS) < 0) *newS = sr->int1;
1075: }/* of neighb[side]==null */
1076: PetscFunctionReturn(PETSC_SUCCESS);
1077: }
1079: /*
1080: Function for sorting an array of real values
1081: */
1082: static PetscErrorCode sortRealEigenvalues(PetscScalar *r,PetscInt *perm,PetscInt nr,PetscBool prev,PetscInt dir)
1083: {
1084: PetscReal re;
1085: PetscInt i,j,tmp;
1087: PetscFunctionBegin;
1088: if (!prev) for (i=0;i<nr;i++) perm[i] = i;
1089: /* Insertion sort */
1090: for (i=1;i<nr;i++) {
1091: re = PetscRealPart(r[perm[i]]);
1092: j = i-1;
1093: while (j>=0 && dir*(re - PetscRealPart(r[perm[j]])) <= 0) {
1094: tmp = perm[j]; perm[j] = perm[j+1]; perm[j+1] = tmp; j--;
1095: }
1096: }
1097: PetscFunctionReturn(PETSC_SUCCESS);
1098: }
1100: /* Stores the pairs obtained since the last shift in the global arrays */
1101: static PetscErrorCode EPSStoreEigenpairs(EPS eps)
1102: {
1103: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1104: PetscReal lambda,err,norm;
1105: PetscInt i,count;
1106: PetscBool iscayley;
1107: EPS_SR sr = ctx->sr;
1108: EPS_shift sPres;
1109: Vec v,w;
1111: PetscFunctionBegin;
1112: sPres = sr->sPres;
1113: sPres->index = sr->indexEig;
1114: count = sr->indexEig;
1115: /* Back-transform */
1116: PetscCall(STBackTransform(eps->st,eps->nconv,eps->eigr,eps->eigi));
1117: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley));
1118: /* Sort eigenvalues */
1119: PetscCall(sortRealEigenvalues(eps->eigr,eps->perm,eps->nconv,PETSC_FALSE,sr->dir));
1120: /* Values stored in global array */
1121: for (i=0;i<eps->nconv;i++) {
1122: lambda = PetscRealPart(eps->eigr[eps->perm[i]]);
1123: err = eps->errest[eps->perm[i]];
1125: if (sr->dir*(lambda - sPres->ext[0]) > 0 && (sr->dir)*(sPres->ext[1] - lambda) > 0) {/* Valid value */
1126: PetscCheck(count<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error in Spectrum Slicing");
1127: sr->eigr[count] = lambda;
1128: sr->errest[count] = err;
1129: /* Explicit purification */
1130: PetscCall(BVGetColumn(eps->V,eps->perm[i],&w));
1131: if (eps->purify) {
1132: PetscCall(BVGetColumn(sr->V,count,&v));
1133: PetscCall(STApply(eps->st,w,v));
1134: PetscCall(BVRestoreColumn(sr->V,count,&v));
1135: } else PetscCall(BVInsertVec(sr->V,count,w));
1136: PetscCall(BVRestoreColumn(eps->V,eps->perm[i],&w));
1137: PetscCall(BVNormColumn(sr->V,count,NORM_2,&norm));
1138: PetscCall(BVScaleColumn(sr->V,count,1.0/norm));
1139: count++;
1140: }
1141: }
1142: sPres->neigs = count - sr->indexEig;
1143: sr->indexEig = count;
1144: /* Global ordering array updating */
1145: PetscCall(sortRealEigenvalues(sr->eigr,sr->perm,count,PETSC_TRUE,sr->dir));
1146: PetscFunctionReturn(PETSC_SUCCESS);
1147: }
1149: static PetscErrorCode EPSLookForDeflation(EPS eps)
1150: {
1151: PetscReal val;
1152: PetscInt i,count0=0,count1=0;
1153: EPS_shift sPres;
1154: PetscInt ini,fin,k,idx0,idx1;
1155: EPS_SR sr;
1156: Vec v;
1157: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1159: PetscFunctionBegin;
1160: sr = ctx->sr;
1161: sPres = sr->sPres;
1163: if (sPres->neighb[0]) ini = (sr->dir)*(sPres->neighb[0]->inertia - sr->inertia0);
1164: else ini = 0;
1165: fin = sr->indexEig;
1166: /* Selection of ends for searching new values */
1167: if (!sPres->neighb[0]) sPres->ext[0] = sr->int0;/* First shift */
1168: else sPres->ext[0] = sPres->neighb[0]->value;
1169: if (!sPres->neighb[1]) {
1170: if (sr->hasEnd) sPres->ext[1] = sr->int1;
1171: else sPres->ext[1] = (sr->dir > 0)?PETSC_MAX_REAL:PETSC_MIN_REAL;
1172: } else sPres->ext[1] = sPres->neighb[1]->value;
1173: /* Selection of values between right and left ends */
1174: for (i=ini;i<fin;i++) {
1175: val=PetscRealPart(sr->eigr[sr->perm[i]]);
1176: /* Values to the right of left shift */
1177: if (sr->dir*(val - sPres->ext[1]) < 0) {
1178: if (sr->dir*(val - sPres->value) < 0) count0++;
1179: else count1++;
1180: } else break;
1181: }
1182: /* The number of values on each side are found */
1183: if (sPres->neighb[0]) {
1184: sPres->nsch[0] = (sr->dir)*(sPres->inertia - sPres->neighb[0]->inertia)-count0;
1185: 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()");
1186: } else sPres->nsch[0] = 0;
1188: if (sPres->neighb[1]) {
1189: sPres->nsch[1] = (sr->dir)*(sPres->neighb[1]->inertia - sPres->inertia) - count1;
1190: 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()");
1191: } else sPres->nsch[1] = (sr->dir)*(sr->inertia1 - sPres->inertia);
1193: /* Completing vector of indexes for deflation */
1194: idx0 = ini;
1195: idx1 = ini+count0+count1;
1196: k=0;
1197: for (i=idx0;i<idx1;i++) sr->idxDef[k++]=sr->perm[i];
1198: PetscCall(BVDuplicateResize(eps->V,k+eps->ncv+1,&sr->Vnext));
1199: PetscCall(BVSetNumConstraints(sr->Vnext,k));
1200: for (i=0;i<k;i++) {
1201: PetscCall(BVGetColumn(sr->Vnext,-i-1,&v));
1202: PetscCall(BVCopyVec(sr->V,sr->idxDef[i],v));
1203: PetscCall(BVRestoreColumn(sr->Vnext,-i-1,&v));
1204: }
1206: /* For rational Krylov */
1207: if (!sr->sPres->rep && sr->nS>0 && (sr->sPrev == sr->sPres->neighb[0] || sr->sPrev == sr->sPres->neighb[1])) PetscCall(EPSPrepareRational(eps));
1208: eps->nconv = 0;
1209: /* Get rid of temporary Vnext */
1210: PetscCall(BVDestroy(&eps->V));
1211: eps->V = sr->Vnext;
1212: sr->Vnext = NULL;
1213: PetscFunctionReturn(PETSC_SUCCESS);
1214: }
1216: PetscErrorCode EPSSolve_KrylovSchur_Slice(EPS eps)
1217: {
1218: PetscInt i,lds,ti;
1219: PetscReal newS;
1220: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1221: EPS_SR sr=ctx->sr;
1222: Mat A,B=NULL;
1223: PetscObjectState Astate,Bstate=0;
1224: PetscObjectId Aid,Bid=0;
1226: PetscFunctionBegin;
1227: PetscCall(PetscCitationsRegister(citation,&cited));
1228: if (ctx->global) {
1229: PetscCall(EPSSolve_KrylovSchur_Slice(ctx->eps));
1230: ctx->eps->state = EPS_STATE_SOLVED;
1231: eps->reason = EPS_CONVERGED_TOL;
1232: if (ctx->npart>1) {
1233: /* Gather solution from subsolvers */
1234: PetscCall(EPSSliceGatherSolution(eps));
1235: } else {
1236: eps->nconv = sr->numEigs;
1237: eps->its = ctx->eps->its;
1238: PetscCall(PetscFree(ctx->inertias));
1239: PetscCall(PetscFree(ctx->shifts));
1240: PetscCall(EPSSliceGetInertias(ctx->eps,&ctx->nshifts,&ctx->shifts,&ctx->inertias));
1241: }
1242: } else {
1243: if (ctx->npart==1) {
1244: sr->eigr = ctx->eps->eigr;
1245: sr->eigi = ctx->eps->eigi;
1246: sr->perm = ctx->eps->perm;
1247: sr->errest = ctx->eps->errest;
1248: sr->V = ctx->eps->V;
1249: }
1250: /* Check that the user did not modify subcomm matrices */
1251: PetscCall(EPSGetOperators(eps,&A,&B));
1252: PetscCall(MatGetState(A,&Astate));
1253: PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1254: if (B) {
1255: PetscCall(MatGetState(B,&Bstate));
1256: PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1257: }
1258: 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");
1259: /* Only with eigenvalues present in the interval ...*/
1260: if (sr->numEigs==0) {
1261: eps->reason = EPS_CONVERGED_TOL;
1262: PetscFunctionReturn(PETSC_SUCCESS);
1263: }
1264: /* Array of pending shifts */
1265: sr->maxPend = 100; /* Initial size */
1266: sr->nPend = 0;
1267: PetscCall(PetscMalloc1(sr->maxPend,&sr->pending));
1268: PetscCall(EPSCreateShift(eps,sr->int0,NULL,NULL));
1269: /* extract first shift */
1270: sr->sPrev = NULL;
1271: sr->sPres = sr->pending[--sr->nPend];
1272: sr->sPres->inertia = sr->inertia0;
1273: eps->target = sr->sPres->value;
1274: sr->s0 = sr->sPres;
1275: sr->indexEig = 0;
1276: /* Memory reservation for auxiliary variables */
1277: lds = PetscMin(eps->mpd,eps->ncv);
1278: PetscCall(PetscCalloc1(lds*lds,&sr->S));
1279: PetscCall(PetscMalloc1(eps->ncv,&sr->back));
1280: for (i=0;i<sr->numEigs;i++) {
1281: sr->eigr[i] = 0.0;
1282: sr->eigi[i] = 0.0;
1283: sr->errest[i] = 0.0;
1284: sr->perm[i] = i;
1285: }
1286: /* Vectors for deflation */
1287: PetscCall(PetscMalloc1(sr->numEigs,&sr->idxDef));
1288: sr->indexEig = 0;
1289: /* Main loop */
1290: while (sr->sPres) {
1291: /* Search for deflation */
1292: PetscCall(EPSLookForDeflation(eps));
1293: /* KrylovSchur */
1294: PetscCall(EPSKrylovSchur_Slice(eps));
1296: PetscCall(EPSStoreEigenpairs(eps));
1297: /* Select new shift */
1298: if (!sr->sPres->comp[1]) {
1299: PetscCall(EPSGetNewShiftValue(eps,1,&newS));
1300: PetscCall(EPSCreateShift(eps,newS,sr->sPres,sr->sPres->neighb[1]));
1301: }
1302: if (!sr->sPres->comp[0]) {
1303: /* Completing earlier interval */
1304: PetscCall(EPSGetNewShiftValue(eps,0,&newS));
1305: PetscCall(EPSCreateShift(eps,newS,sr->sPres->neighb[0],sr->sPres));
1306: }
1307: /* Preparing for a new search of values */
1308: PetscCall(EPSExtractShift(eps));
1309: }
1311: /* Updating eps values prior to exit */
1312: PetscCall(PetscFree(sr->S));
1313: PetscCall(PetscFree(sr->idxDef));
1314: PetscCall(PetscFree(sr->pending));
1315: PetscCall(PetscFree(sr->back));
1316: PetscCall(BVDuplicateResize(eps->V,eps->ncv+1,&sr->Vnext));
1317: PetscCall(BVSetNumConstraints(sr->Vnext,0));
1318: PetscCall(BVDestroy(&eps->V));
1319: eps->V = sr->Vnext;
1320: eps->nconv = sr->indexEig;
1321: eps->reason = EPS_CONVERGED_TOL;
1322: eps->its = sr->itsKs;
1323: eps->nds = 0;
1324: if (sr->dir<0) {
1325: for (i=0;i<eps->nconv/2;i++) {
1326: ti = sr->perm[i]; sr->perm[i] = sr->perm[eps->nconv-1-i]; sr->perm[eps->nconv-1-i] = ti;
1327: }
1328: }
1329: }
1330: PetscFunctionReturn(PETSC_SUCCESS);
1331: }