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