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