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