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