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
14 :
15 : Algorithm:
16 :
17 : Single-vector Krylov-Schur method for non-symmetric problems,
18 : including harmonic extraction.
19 :
20 : References:
21 :
22 : [1] "Krylov-Schur Methods in SLEPc", SLEPc Technical Report STR-7,
23 : available at https://slepc.upv.es.
24 :
25 : [2] G.W. Stewart, "A Krylov-Schur Algorithm for Large Eigenproblems",
26 : SIAM J. Matrix Anal. App. 23(3):601-614, 2001.
27 :
28 : [3] "Practical Implementation of Harmonic Krylov-Schur", SLEPc Technical
29 : Report STR-9, available at https://slepc.upv.es.
30 : */
31 :
32 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
33 : #include "krylovschur.h"
34 :
35 20 : PetscErrorCode EPSGetArbitraryValues(EPS eps,PetscScalar *rr,PetscScalar *ri)
36 : {
37 20 : PetscInt i,newi,ld,n,l;
38 20 : Vec xr=eps->work[0],xi=eps->work[1];
39 20 : PetscScalar re,im,*Zr,*Zi,*X;
40 :
41 20 : PetscFunctionBegin;
42 20 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
43 20 : PetscCall(DSGetDimensions(eps->ds,&n,&l,NULL,NULL));
44 340 : for (i=l;i<n;i++) {
45 320 : re = eps->eigr[i];
46 320 : im = eps->eigi[i];
47 320 : PetscCall(STBackTransform(eps->st,1,&re,&im));
48 320 : newi = i;
49 320 : PetscCall(DSVectors(eps->ds,DS_MAT_X,&newi,NULL));
50 320 : PetscCall(DSGetArray(eps->ds,DS_MAT_X,&X));
51 320 : Zr = X+i*ld;
52 320 : if (newi==i+1) Zi = X+newi*ld;
53 : else Zi = NULL;
54 320 : PetscCall(EPSComputeRitzVector(eps,Zr,Zi,eps->V,xr,xi));
55 320 : PetscCall(DSRestoreArray(eps->ds,DS_MAT_X,&X));
56 320 : PetscCall((*eps->arbitrary)(re,im,xr,xi,rr+i,ri+i,eps->arbitraryctx));
57 : }
58 20 : PetscFunctionReturn(PETSC_SUCCESS);
59 : }
60 :
61 8 : static PetscErrorCode EPSSetUp_KrylovSchur_Filter(EPS eps)
62 : {
63 8 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
64 8 : PetscBool estimaterange=PETSC_TRUE;
65 8 : PetscReal rleft,rright;
66 8 : Mat A;
67 :
68 8 : PetscFunctionBegin;
69 8 : EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
70 8 : EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
71 8 : 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");
72 8 : EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_THRESHOLD,PETSC_TRUE," with polynomial filter");
73 8 : if (eps->tol==(PetscReal)PETSC_DETERMINE) eps->tol = SLEPC_DEFAULT_TOL*1e-2; /* use tighter tolerance */
74 8 : PetscCall(STFilterSetInterval(eps->st,eps->inta,eps->intb));
75 8 : if (!ctx->estimatedrange) {
76 6 : PetscCall(STFilterGetRange(eps->st,&rleft,&rright));
77 6 : estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
78 : }
79 6 : if (estimaterange) { /* user did not set a range */
80 6 : PetscCall(STGetMatrix(eps->st,0,&A));
81 6 : PetscCall(MatEstimateSpectralRange_EPS(A,&rleft,&rright));
82 6 : PetscCall(PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright));
83 6 : PetscCall(STFilterSetRange(eps->st,rleft,rright));
84 6 : ctx->estimatedrange = PETSC_TRUE;
85 : }
86 8 : if (eps->ncv==PETSC_DETERMINE && eps->nev==0) eps->nev = 40; /* user did not provide nev estimation */
87 8 : PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
88 8 : PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
89 8 : if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
90 8 : PetscFunctionReturn(PETSC_SUCCESS);
91 : }
92 :
93 526 : static PetscErrorCode EPSSetUp_KrylovSchur(EPS eps)
94 : {
95 526 : PetscReal eta;
96 526 : PetscBool isfilt=PETSC_FALSE;
97 526 : BVOrthogType otype;
98 526 : BVOrthogBlockType obtype;
99 526 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
100 526 : enum { EPS_KS_DEFAULT,EPS_KS_SYMM,EPS_KS_SLICE,EPS_KS_FILTER,EPS_KS_INDEF,EPS_KS_TWOSIDED } variant;
101 :
102 526 : PetscFunctionBegin;
103 526 : if (eps->which==EPS_ALL) { /* default values in case of spectrum slicing or polynomial filter */
104 44 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
105 44 : if (isfilt) PetscCall(EPSSetUp_KrylovSchur_Filter(eps));
106 36 : else PetscCall(EPSSetUp_KrylovSchur_Slice(eps));
107 482 : } else if (eps->isstructured) {
108 20 : PetscCall(EPSSetUp_KrylovSchur_BSE(eps));
109 20 : PetscFunctionReturn(PETSC_SUCCESS);
110 : } else {
111 462 : PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
112 462 : PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
113 687 : if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);
114 462 : if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
115 : }
116 506 : PetscCheck(ctx->lock || eps->mpd>=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
117 :
118 506 : EPSCheckDefiniteCondition(eps,eps->arbitrary," with arbitrary selection of eigenpairs");
119 :
120 506 : PetscCheck(eps->extraction==EPS_RITZ || eps->extraction==EPS_HARMONIC,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
121 :
122 506 : if (!ctx->keep) ctx->keep = 0.5;
123 :
124 506 : PetscCall(EPSAllocateSolution(eps,1));
125 506 : PetscCall(EPS_SetInnerProduct(eps));
126 506 : if (eps->arbitrary) PetscCall(EPSSetWorkVecs(eps,2));
127 504 : else if (eps->ishermitian && !eps->ispositive) PetscCall(EPSSetWorkVecs(eps,1));
128 :
129 : /* dispatch solve method */
130 506 : if (eps->ishermitian) {
131 261 : if (eps->which==EPS_ALL) {
132 44 : EPSCheckDefiniteCondition(eps,eps->which==EPS_ALL," with spectrum slicing");
133 44 : variant = isfilt? EPS_KS_FILTER: EPS_KS_SLICE;
134 217 : } else if (eps->isgeneralized && !eps->ispositive) {
135 : variant = EPS_KS_INDEF;
136 : } else {
137 201 : switch (eps->extraction) {
138 : case EPS_RITZ: variant = EPS_KS_SYMM; break;
139 : case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
140 0 : default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
141 : }
142 : }
143 245 : } else if (eps->twosided) {
144 : variant = EPS_KS_TWOSIDED;
145 : } else {
146 232 : switch (eps->extraction) {
147 : case EPS_RITZ: variant = EPS_KS_DEFAULT; break;
148 : case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
149 0 : default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
150 : }
151 : }
152 65 : switch (variant) {
153 237 : case EPS_KS_DEFAULT:
154 237 : eps->ops->solve = EPSSolve_KrylovSchur_Default;
155 237 : eps->ops->computevectors = EPSComputeVectors_Schur;
156 237 : PetscCall(DSSetType(eps->ds,DSNHEP));
157 237 : PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
158 237 : PetscCall(DSAllocate(eps->ds,eps->ncv+1));
159 : break;
160 204 : case EPS_KS_SYMM:
161 : case EPS_KS_FILTER:
162 204 : eps->ops->solve = EPSSolve_KrylovSchur_Default;
163 204 : eps->ops->computevectors = EPSComputeVectors_Hermitian;
164 204 : PetscCall(DSSetType(eps->ds,DSHEP));
165 204 : PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
166 204 : PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
167 204 : PetscCall(DSAllocate(eps->ds,eps->ncv+1));
168 : break;
169 : case EPS_KS_SLICE:
170 36 : eps->ops->solve = EPSSolve_KrylovSchur_Slice;
171 36 : eps->ops->computevectors = EPSComputeVectors_Slice;
172 36 : break;
173 : case EPS_KS_INDEF:
174 16 : eps->ops->solve = EPSSolve_KrylovSchur_Indefinite;
175 16 : eps->ops->computevectors = EPSComputeVectors_Indefinite;
176 16 : PetscCall(DSSetType(eps->ds,DSGHIEP));
177 16 : PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
178 16 : PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
179 16 : PetscCall(DSAllocate(eps->ds,eps->ncv+1));
180 : /* force reorthogonalization for pseudo-Lanczos */
181 16 : PetscCall(BVGetOrthogonalization(eps->V,&otype,NULL,&eta,&obtype));
182 16 : PetscCall(BVSetOrthogonalization(eps->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
183 : break;
184 : case EPS_KS_TWOSIDED:
185 13 : eps->ops->solve = EPSSolve_KrylovSchur_TwoSided;
186 13 : eps->ops->computevectors = EPSComputeVectors_Schur;
187 13 : PetscCall(DSSetType(eps->ds,DSNHEPTS));
188 13 : PetscCall(DSAllocate(eps->ds,eps->ncv+1));
189 13 : PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
190 : break;
191 506 : default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
192 : }
193 506 : PetscFunctionReturn(PETSC_SUCCESS);
194 : }
195 :
196 526 : static PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
197 : {
198 526 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
199 526 : SlepcSC sc;
200 526 : PetscBool isfilt;
201 :
202 526 : PetscFunctionBegin;
203 526 : PetscCall(EPSSetUpSort_Default(eps));
204 526 : if (eps->which==EPS_ALL) {
205 44 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
206 44 : if (isfilt) {
207 8 : PetscCall(DSGetSlepcSC(eps->ds,&sc));
208 8 : sc->rg = NULL;
209 8 : sc->comparison = SlepcCompareLargestReal;
210 8 : sc->comparisonctx = NULL;
211 8 : sc->map = NULL;
212 8 : sc->mapobj = NULL;
213 : } else {
214 36 : if (!ctx->global && ctx->sr->numEigs>0) {
215 18 : PetscCall(DSGetSlepcSC(eps->ds,&sc));
216 18 : sc->rg = NULL;
217 18 : sc->comparison = SlepcCompareLargestMagnitude;
218 18 : sc->comparisonctx = NULL;
219 18 : sc->map = NULL;
220 18 : sc->mapobj = NULL;
221 : }
222 : }
223 : }
224 526 : PetscFunctionReturn(PETSC_SUCCESS);
225 : }
226 :
227 441 : PetscErrorCode EPSSolve_KrylovSchur_Default(EPS eps)
228 : {
229 441 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
230 441 : PetscInt i,j,*pj,k,l,nv,ld,nconv;
231 441 : Mat U,Op,H,T;
232 441 : PetscScalar *g;
233 441 : PetscReal beta,gamma=1.0;
234 441 : PetscBool breakdown,harmonic,hermitian;
235 :
236 441 : PetscFunctionBegin;
237 441 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
238 441 : harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
239 441 : hermitian = (eps->ishermitian && !harmonic)?PETSC_TRUE:PETSC_FALSE;
240 441 : if (harmonic) PetscCall(PetscMalloc1(ld,&g));
241 441 : if (eps->arbitrary) pj = &j;
242 439 : else pj = NULL;
243 :
244 : /* Get the starting Arnoldi vector */
245 441 : PetscCall(EPSGetStartVector(eps,0,NULL));
246 441 : l = 0;
247 :
248 : /* Restart loop */
249 441 : while (eps->reason == EPS_CONVERGED_ITERATING) {
250 3389 : eps->its++;
251 :
252 : /* Compute an nv-step Arnoldi factorization */
253 3389 : nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
254 3389 : PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
255 3389 : PetscCall(STGetOperator(eps->st,&Op));
256 3389 : if (hermitian) {
257 2345 : PetscCall(DSGetMat(eps->ds,DS_MAT_T,&T));
258 2345 : PetscCall(BVMatLanczos(eps->V,Op,T,eps->nconv+l,&nv,&beta,&breakdown));
259 2345 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&T));
260 : } else {
261 1044 : PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
262 1044 : PetscCall(BVMatArnoldi(eps->V,Op,H,eps->nconv+l,&nv,&beta,&breakdown));
263 1044 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
264 : }
265 3389 : PetscCall(STRestoreOperator(eps->st,&Op));
266 3389 : PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
267 3389 : PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
268 3389 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
269 :
270 : /* Compute translation of Krylov decomposition if harmonic extraction used */
271 3389 : if (PetscUnlikely(harmonic)) PetscCall(DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,g,&gamma));
272 :
273 : /* Solve projected problem */
274 3389 : PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
275 3389 : if (PetscUnlikely(eps->arbitrary)) {
276 20 : PetscCall(EPSGetArbitraryValues(eps,eps->rr,eps->ri));
277 20 : j=1;
278 : }
279 3389 : PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj));
280 3389 : PetscCall(DSUpdateExtraRow(eps->ds));
281 3389 : PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
282 :
283 : /* Check convergence */
284 3389 : PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
285 3389 : EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
286 3389 : PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
287 3389 : nconv = k;
288 :
289 : /* Update l */
290 3389 : if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
291 : else {
292 2948 : l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
293 2948 : if (!hermitian) PetscCall(DSGetTruncateSize(eps->ds,k,nv,&l));
294 : }
295 3389 : if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
296 3389 : if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
297 :
298 3389 : if (eps->reason == EPS_CONVERGED_ITERATING) {
299 2948 : if (PetscUnlikely(breakdown || k==nv)) {
300 : /* Start a new Arnoldi factorization */
301 0 : PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
302 0 : if (k<eps->nev) {
303 0 : PetscCall(EPSGetStartVector(eps,k,&breakdown));
304 0 : if (breakdown) {
305 0 : eps->reason = EPS_DIVERGED_BREAKDOWN;
306 0 : PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
307 : }
308 : }
309 : } else {
310 : /* Undo translation of Krylov decomposition */
311 2948 : if (PetscUnlikely(harmonic)) {
312 24 : PetscCall(DSSetDimensions(eps->ds,nv,k,l));
313 24 : PetscCall(DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma));
314 : /* gamma u^ = u - U*g~ */
315 24 : PetscCall(BVSetActiveColumns(eps->V,0,nv));
316 24 : PetscCall(BVMultColumn(eps->V,-1.0,1.0,nv,g));
317 24 : PetscCall(BVScaleColumn(eps->V,nv,1.0/gamma));
318 24 : PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
319 24 : PetscCall(DSSetDimensions(eps->ds,nv,k,nv));
320 : }
321 : /* Prepare the Rayleigh quotient for restart */
322 2948 : PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
323 : }
324 : }
325 : /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
326 3389 : PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
327 3389 : PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
328 3389 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
329 :
330 3389 : if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
331 2948 : PetscCall(BVCopyColumn(eps->V,nv,k+l)); /* copy restart vector from the last column */
332 2948 : if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
333 3 : eps->ncv = eps->mpd+k;
334 3 : PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
335 27 : for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
336 3 : PetscCall(DSReallocate(eps->ds,eps->ncv+1));
337 3 : PetscCall(DSGetLeadingDimension(eps->ds,&ld));
338 : }
339 : }
340 :
341 3389 : eps->nconv = k;
342 3830 : PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
343 : }
344 :
345 441 : if (harmonic) PetscCall(PetscFree(g));
346 441 : PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
347 441 : PetscFunctionReturn(PETSC_SUCCESS);
348 : }
349 :
350 6 : static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
351 : {
352 6 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
353 :
354 6 : PetscFunctionBegin;
355 6 : if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
356 : else {
357 6 : PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument %g must be in the range [0.1,0.9]",(double)keep);
358 6 : ctx->keep = keep;
359 : }
360 6 : PetscFunctionReturn(PETSC_SUCCESS);
361 : }
362 :
363 : /*@
364 : EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
365 : method, in particular the proportion of basis vectors that must be kept
366 : after restart.
367 :
368 : Logically Collective
369 :
370 : Input Parameters:
371 : + eps - the eigenproblem solver context
372 : - keep - the number of vectors to be kept at restart
373 :
374 : Options Database Key:
375 : . -eps_krylovschur_restart - Sets the restart parameter
376 :
377 : Notes:
378 : Allowed values are in the range [0.1,0.9]. The default is 0.5.
379 :
380 : Level: advanced
381 :
382 : .seealso: EPSKrylovSchurGetRestart()
383 : @*/
384 6 : PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
385 : {
386 6 : PetscFunctionBegin;
387 6 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
388 18 : PetscValidLogicalCollectiveReal(eps,keep,2);
389 6 : PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
390 6 : PetscFunctionReturn(PETSC_SUCCESS);
391 : }
392 :
393 6 : static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
394 : {
395 6 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
396 :
397 6 : PetscFunctionBegin;
398 6 : *keep = ctx->keep;
399 6 : PetscFunctionReturn(PETSC_SUCCESS);
400 : }
401 :
402 : /*@
403 : EPSKrylovSchurGetRestart - Gets the restart parameter used in the
404 : Krylov-Schur method.
405 :
406 : Not Collective
407 :
408 : Input Parameter:
409 : . eps - the eigenproblem solver context
410 :
411 : Output Parameter:
412 : . keep - the restart parameter
413 :
414 : Level: advanced
415 :
416 : .seealso: EPSKrylovSchurSetRestart()
417 : @*/
418 6 : PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
419 : {
420 6 : PetscFunctionBegin;
421 6 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
422 6 : PetscAssertPointer(keep,2);
423 6 : PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
424 6 : PetscFunctionReturn(PETSC_SUCCESS);
425 : }
426 :
427 39 : static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
428 : {
429 39 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
430 :
431 39 : PetscFunctionBegin;
432 39 : ctx->lock = lock;
433 39 : PetscFunctionReturn(PETSC_SUCCESS);
434 : }
435 :
436 : /*@
437 : EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
438 : the Krylov-Schur method.
439 :
440 : Logically Collective
441 :
442 : Input Parameters:
443 : + eps - the eigenproblem solver context
444 : - lock - true if the locking variant must be selected
445 :
446 : Options Database Key:
447 : . -eps_krylovschur_locking - Sets the locking flag
448 :
449 : Notes:
450 : The default is to lock converged eigenpairs when the method restarts.
451 : This behaviour can be changed so that all directions are kept in the
452 : working subspace even if already converged to working accuracy (the
453 : non-locking variant).
454 :
455 : Level: advanced
456 :
457 : .seealso: EPSKrylovSchurGetLocking()
458 : @*/
459 39 : PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
460 : {
461 39 : PetscFunctionBegin;
462 39 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
463 117 : PetscValidLogicalCollectiveBool(eps,lock,2);
464 39 : PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
465 39 : PetscFunctionReturn(PETSC_SUCCESS);
466 : }
467 :
468 6 : static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
469 : {
470 6 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
471 :
472 6 : PetscFunctionBegin;
473 6 : *lock = ctx->lock;
474 6 : PetscFunctionReturn(PETSC_SUCCESS);
475 : }
476 :
477 : /*@
478 : EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
479 : method.
480 :
481 : Not Collective
482 :
483 : Input Parameter:
484 : . eps - the eigenproblem solver context
485 :
486 : Output Parameter:
487 : . lock - the locking flag
488 :
489 : Level: advanced
490 :
491 : .seealso: EPSKrylovSchurSetLocking()
492 : @*/
493 6 : PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
494 : {
495 6 : PetscFunctionBegin;
496 6 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
497 6 : PetscAssertPointer(lock,2);
498 6 : PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
499 6 : PetscFunctionReturn(PETSC_SUCCESS);
500 : }
501 :
502 5 : static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
503 : {
504 5 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
505 5 : PetscMPIInt size;
506 5 : PetscInt newnpart;
507 :
508 5 : PetscFunctionBegin;
509 5 : if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
510 : newnpart = 1;
511 : } else {
512 5 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
513 5 : PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
514 : newnpart = npart;
515 : }
516 5 : if (ctx->npart!=newnpart) {
517 5 : if (ctx->npart>1) {
518 0 : PetscCall(PetscSubcommDestroy(&ctx->subc));
519 0 : if (ctx->commset) {
520 0 : PetscCallMPI(MPI_Comm_free(&ctx->commrank));
521 0 : ctx->commset = PETSC_FALSE;
522 : }
523 : }
524 5 : PetscCall(EPSDestroy(&ctx->eps));
525 5 : ctx->npart = newnpart;
526 5 : eps->state = EPS_STATE_INITIAL;
527 : }
528 5 : PetscFunctionReturn(PETSC_SUCCESS);
529 : }
530 :
531 : /*@
532 : EPSKrylovSchurSetPartitions - Sets the number of partitions for the
533 : case of doing spectrum slicing for a computational interval with the
534 : communicator split in several sub-communicators.
535 :
536 : Logically Collective
537 :
538 : Input Parameters:
539 : + eps - the eigenproblem solver context
540 : - npart - number of partitions
541 :
542 : Options Database Key:
543 : . -eps_krylovschur_partitions <npart> - Sets the number of partitions
544 :
545 : Notes:
546 : By default, npart=1 so all processes in the communicator participate in
547 : the processing of the whole interval. If npart>1 then the interval is
548 : divided into npart subintervals, each of them being processed by a
549 : subset of processes.
550 :
551 : The interval is split proportionally unless the separation points are
552 : specified with EPSKrylovSchurSetSubintervals().
553 :
554 : Level: advanced
555 :
556 : .seealso: EPSKrylovSchurSetSubintervals(), EPSSetInterval()
557 : @*/
558 5 : PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
559 : {
560 5 : PetscFunctionBegin;
561 5 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
562 15 : PetscValidLogicalCollectiveInt(eps,npart,2);
563 5 : PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
564 5 : PetscFunctionReturn(PETSC_SUCCESS);
565 : }
566 :
567 2 : static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
568 : {
569 2 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
570 :
571 2 : PetscFunctionBegin;
572 2 : *npart = ctx->npart;
573 2 : PetscFunctionReturn(PETSC_SUCCESS);
574 : }
575 :
576 : /*@
577 : EPSKrylovSchurGetPartitions - Gets the number of partitions of the
578 : communicator in case of spectrum slicing.
579 :
580 : Not Collective
581 :
582 : Input Parameter:
583 : . eps - the eigenproblem solver context
584 :
585 : Output Parameter:
586 : . npart - number of partitions
587 :
588 : Level: advanced
589 :
590 : .seealso: EPSKrylovSchurSetPartitions()
591 : @*/
592 2 : PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
593 : {
594 2 : PetscFunctionBegin;
595 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
596 2 : PetscAssertPointer(npart,2);
597 2 : PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
598 2 : PetscFunctionReturn(PETSC_SUCCESS);
599 : }
600 :
601 3 : static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
602 : {
603 3 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
604 :
605 3 : PetscFunctionBegin;
606 3 : ctx->detect = detect;
607 3 : eps->state = EPS_STATE_INITIAL;
608 3 : PetscFunctionReturn(PETSC_SUCCESS);
609 : }
610 :
611 : /*@
612 : EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
613 : zeros during the factorizations throughout the spectrum slicing computation.
614 :
615 : Logically Collective
616 :
617 : Input Parameters:
618 : + eps - the eigenproblem solver context
619 : - detect - check for zeros
620 :
621 : Options Database Key:
622 : . -eps_krylovschur_detect_zeros - Check for zeros; this takes an optional
623 : bool value (0/1/no/yes/true/false)
624 :
625 : Notes:
626 : A zero in the factorization indicates that a shift coincides with an eigenvalue.
627 :
628 : This flag is turned off by default, and may be necessary in some cases,
629 : especially when several partitions are being used. This feature currently
630 : requires an external package for factorizations with support for zero
631 : detection, e.g. MUMPS.
632 :
633 : Level: advanced
634 :
635 : .seealso: EPSKrylovSchurSetPartitions(), EPSSetInterval()
636 : @*/
637 3 : PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
638 : {
639 3 : PetscFunctionBegin;
640 3 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
641 9 : PetscValidLogicalCollectiveBool(eps,detect,2);
642 3 : PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
643 3 : PetscFunctionReturn(PETSC_SUCCESS);
644 : }
645 :
646 6 : static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
647 : {
648 6 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
649 :
650 6 : PetscFunctionBegin;
651 6 : *detect = ctx->detect;
652 6 : PetscFunctionReturn(PETSC_SUCCESS);
653 : }
654 :
655 : /*@
656 : EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
657 : in spectrum slicing.
658 :
659 : Not Collective
660 :
661 : Input Parameter:
662 : . eps - the eigenproblem solver context
663 :
664 : Output Parameter:
665 : . detect - whether zeros detection is enforced during factorizations
666 :
667 : Level: advanced
668 :
669 : .seealso: EPSKrylovSchurSetDetectZeros()
670 : @*/
671 6 : PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
672 : {
673 6 : PetscFunctionBegin;
674 6 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
675 6 : PetscAssertPointer(detect,2);
676 6 : PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
677 6 : PetscFunctionReturn(PETSC_SUCCESS);
678 : }
679 :
680 3 : static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
681 : {
682 3 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
683 :
684 3 : PetscFunctionBegin;
685 3 : if (nev != PETSC_CURRENT) {
686 3 : PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
687 3 : ctx->nev = nev;
688 : }
689 3 : if (ncv == PETSC_DETERMINE) {
690 0 : ctx->ncv = PETSC_DETERMINE;
691 3 : } else if (ncv != PETSC_CURRENT) {
692 3 : PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
693 3 : ctx->ncv = ncv;
694 : }
695 3 : if (mpd == PETSC_DETERMINE) {
696 0 : ctx->mpd = PETSC_DETERMINE;
697 3 : } else if (mpd != PETSC_CURRENT) {
698 3 : PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
699 3 : ctx->mpd = mpd;
700 : }
701 3 : eps->state = EPS_STATE_INITIAL;
702 3 : PetscFunctionReturn(PETSC_SUCCESS);
703 : }
704 :
705 : /*@
706 : EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
707 : step in case of doing spectrum slicing for a computational interval.
708 : The meaning of the parameters is the same as in EPSSetDimensions().
709 :
710 : Logically Collective
711 :
712 : Input Parameters:
713 : + eps - the eigenproblem solver context
714 : . nev - number of eigenvalues to compute
715 : . ncv - the maximum dimension of the subspace to be used by the subsolve
716 : - mpd - the maximum dimension allowed for the projected problem
717 :
718 : Options Database Key:
719 : + -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
720 : . -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
721 : - -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension
722 :
723 : Note:
724 : Use PETSC_DETERMINE for ncv and mpd to assign a default value. For any
725 : of the arguments, use PETSC_CURRENT to preserve the current value.
726 :
727 : Level: advanced
728 :
729 : .seealso: EPSKrylovSchurGetDimensions(), EPSSetDimensions(), EPSSetInterval()
730 : @*/
731 3 : PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
732 : {
733 3 : PetscFunctionBegin;
734 3 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
735 9 : PetscValidLogicalCollectiveInt(eps,nev,2);
736 9 : PetscValidLogicalCollectiveInt(eps,ncv,3);
737 9 : PetscValidLogicalCollectiveInt(eps,mpd,4);
738 3 : PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
739 3 : PetscFunctionReturn(PETSC_SUCCESS);
740 : }
741 :
742 6 : static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
743 : {
744 6 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
745 :
746 6 : PetscFunctionBegin;
747 6 : if (nev) *nev = ctx->nev;
748 6 : if (ncv) *ncv = ctx->ncv;
749 6 : if (mpd) *mpd = ctx->mpd;
750 6 : PetscFunctionReturn(PETSC_SUCCESS);
751 : }
752 :
753 : /*@
754 : EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
755 : step in case of doing spectrum slicing for a computational interval.
756 :
757 : Not Collective
758 :
759 : Input Parameter:
760 : . eps - the eigenproblem solver context
761 :
762 : Output Parameters:
763 : + nev - number of eigenvalues to compute
764 : . ncv - the maximum dimension of the subspace to be used by the subsolve
765 : - mpd - the maximum dimension allowed for the projected problem
766 :
767 : Level: advanced
768 :
769 : .seealso: EPSKrylovSchurSetDimensions()
770 : @*/
771 6 : PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
772 : {
773 6 : PetscFunctionBegin;
774 6 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
775 6 : PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
776 6 : PetscFunctionReturn(PETSC_SUCCESS);
777 : }
778 :
779 2 : static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
780 : {
781 2 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
782 2 : PetscInt i;
783 :
784 2 : PetscFunctionBegin;
785 2 : PetscCheck(subint[0]==eps->inta && subint[ctx->npart]==eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"First and last values must match the endpoints of EPSSetInterval()");
786 6 : for (i=0;i<ctx->npart;i++) PetscCheck(subint[i]<=subint[i+1],PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Array must contain values in ascending order");
787 2 : if (ctx->subintervals) PetscCall(PetscFree(ctx->subintervals));
788 2 : PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
789 8 : for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
790 2 : ctx->subintset = PETSC_TRUE;
791 2 : eps->state = EPS_STATE_INITIAL;
792 2 : PetscFunctionReturn(PETSC_SUCCESS);
793 : }
794 :
795 : /*@
796 : EPSKrylovSchurSetSubintervals - Sets the points that delimit the
797 : subintervals to be used in spectrum slicing with several partitions.
798 :
799 : Logically Collective
800 :
801 : Input Parameters:
802 : + eps - the eigenproblem solver context
803 : - subint - array of real values specifying subintervals
804 :
805 : Notes:
806 : This function must be called after EPSKrylovSchurSetPartitions(). For npart
807 : partitions, the argument subint must contain npart+1 real values sorted in
808 : ascending order, subint_0, subint_1, ..., subint_npart, where the first
809 : and last values must coincide with the interval endpoints set with
810 : EPSSetInterval().
811 :
812 : The subintervals are then defined by two consecutive points [subint_0,subint_1],
813 : [subint_1,subint_2], and so on.
814 :
815 : Level: advanced
816 :
817 : .seealso: EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubintervals(), EPSSetInterval()
818 : @*/
819 2 : PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal subint[])
820 : {
821 2 : PetscFunctionBegin;
822 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
823 2 : PetscAssertPointer(subint,2);
824 2 : PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
825 2 : PetscFunctionReturn(PETSC_SUCCESS);
826 : }
827 :
828 2 : static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal **subint)
829 : {
830 2 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
831 2 : PetscInt i;
832 :
833 2 : PetscFunctionBegin;
834 2 : if (!ctx->subintset) {
835 0 : PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
836 0 : PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
837 : }
838 2 : PetscCall(PetscMalloc1(ctx->npart+1,subint));
839 8 : for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
840 2 : PetscFunctionReturn(PETSC_SUCCESS);
841 : }
842 :
843 : /*@C
844 : EPSKrylovSchurGetSubintervals - Returns the points that delimit the
845 : subintervals used in spectrum slicing with several partitions.
846 :
847 : Not Collective
848 :
849 : Input Parameter:
850 : . eps - the eigenproblem solver context
851 :
852 : Output Parameter:
853 : . subint - array of real values specifying subintervals
854 :
855 : Notes:
856 : If the user passed values with EPSKrylovSchurSetSubintervals(), then the
857 : same values are returned. Otherwise, the values computed internally are
858 : obtained.
859 :
860 : This function is only available for spectrum slicing runs.
861 :
862 : The returned array has length npart+1 (see EPSKrylovSchurGetPartitions())
863 : and should be freed by the user.
864 :
865 : Fortran Notes:
866 : The calling sequence from Fortran is
867 : .vb
868 : EPSKrylovSchurGetSubintervals(eps,subint,ierr)
869 : double precision subint(npart+1) output
870 : .ve
871 :
872 : Level: advanced
873 :
874 : .seealso: EPSKrylovSchurSetSubintervals(), EPSKrylovSchurGetPartitions(), EPSSetInterval()
875 : @*/
876 2 : PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal **subint)
877 : {
878 2 : PetscFunctionBegin;
879 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
880 2 : PetscAssertPointer(subint,2);
881 2 : PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
882 2 : PetscFunctionReturn(PETSC_SUCCESS);
883 : }
884 :
885 3 : static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
886 : {
887 3 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
888 3 : PetscInt i,numsh;
889 3 : EPS_SR sr = ctx->sr;
890 :
891 3 : PetscFunctionBegin;
892 3 : PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
893 3 : PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
894 3 : switch (eps->state) {
895 : case EPS_STATE_INITIAL:
896 : break;
897 3 : case EPS_STATE_SETUP:
898 3 : numsh = ctx->npart+1;
899 3 : if (n) *n = numsh;
900 3 : if (shifts) {
901 3 : PetscCall(PetscMalloc1(numsh,shifts));
902 3 : (*shifts)[0] = eps->inta;
903 3 : if (ctx->npart==1) (*shifts)[1] = eps->intb;
904 6 : else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
905 : }
906 3 : if (inertias) {
907 3 : PetscCall(PetscMalloc1(numsh,inertias));
908 3 : (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
909 3 : if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
910 6 : else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
911 : }
912 : break;
913 0 : case EPS_STATE_SOLVED:
914 : case EPS_STATE_EIGENVECTORS:
915 0 : numsh = ctx->nshifts;
916 0 : if (n) *n = numsh;
917 0 : if (shifts) {
918 0 : PetscCall(PetscMalloc1(numsh,shifts));
919 0 : for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
920 : }
921 0 : if (inertias) {
922 0 : PetscCall(PetscMalloc1(numsh,inertias));
923 0 : for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
924 : }
925 : break;
926 : }
927 3 : PetscFunctionReturn(PETSC_SUCCESS);
928 : }
929 :
930 : /*@C
931 : EPSKrylovSchurGetInertias - Gets the values of the shifts and their
932 : corresponding inertias in case of doing spectrum slicing for a
933 : computational interval.
934 :
935 : Not Collective
936 :
937 : Input Parameter:
938 : . eps - the eigenproblem solver context
939 :
940 : Output Parameters:
941 : + n - number of shifts, including the endpoints of the interval
942 : . shifts - the values of the shifts used internally in the solver
943 : - inertias - the values of the inertia in each shift
944 :
945 : Notes:
946 : If called after EPSSolve(), all shifts used internally by the solver are
947 : returned (including both endpoints and any intermediate ones). If called
948 : before EPSSolve() and after EPSSetUp() then only the information of the
949 : endpoints of subintervals is available.
950 :
951 : This function is only available for spectrum slicing runs.
952 :
953 : The returned arrays should be freed by the user. Can pass NULL in any of
954 : the two arrays if not required.
955 :
956 : Fortran Notes:
957 : The calling sequence from Fortran is
958 : .vb
959 : EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
960 : integer n
961 : double precision shifts(*)
962 : integer inertias(*)
963 : .ve
964 : The arrays should be at least of length n. The value of n can be determined
965 : by an initial call
966 : .vb
967 : EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr)
968 : .ve
969 :
970 : Level: advanced
971 :
972 : .seealso: EPSSetInterval(), EPSKrylovSchurSetSubintervals()
973 : @*/
974 3 : PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[])
975 : {
976 3 : PetscFunctionBegin;
977 3 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
978 3 : PetscAssertPointer(n,2);
979 3 : PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
980 3 : PetscFunctionReturn(PETSC_SUCCESS);
981 : }
982 :
983 2 : static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
984 : {
985 2 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
986 2 : EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
987 :
988 2 : PetscFunctionBegin;
989 2 : PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
990 2 : PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
991 2 : if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
992 2 : if (n) *n = sr->numEigs;
993 2 : if (v) PetscCall(BVCreateVec(sr->V,v));
994 2 : PetscFunctionReturn(PETSC_SUCCESS);
995 : }
996 :
997 : /*@
998 : EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
999 : doing spectrum slicing for a computational interval with multiple
1000 : communicators.
1001 :
1002 : Collective on the subcommunicator (if v is given)
1003 :
1004 : Input Parameter:
1005 : . eps - the eigenproblem solver context
1006 :
1007 : Output Parameters:
1008 : + k - index of the subinterval for the calling process
1009 : . n - number of eigenvalues found in the k-th subinterval
1010 : - v - a vector owned by processes in the subcommunicator with dimensions
1011 : compatible for locally computed eigenvectors (or NULL)
1012 :
1013 : Notes:
1014 : This function is only available for spectrum slicing runs.
1015 :
1016 : The returned Vec should be destroyed by the user.
1017 :
1018 : Level: advanced
1019 :
1020 : .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommPairs()
1021 : @*/
1022 2 : PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1023 : {
1024 2 : PetscFunctionBegin;
1025 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1026 2 : PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1027 2 : PetscFunctionReturn(PETSC_SUCCESS);
1028 : }
1029 :
1030 58 : static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1031 : {
1032 58 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1033 58 : EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1034 :
1035 58 : PetscFunctionBegin;
1036 58 : EPSCheckSolved(eps,1);
1037 58 : PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1038 58 : PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1039 58 : if (eig) *eig = sr->eigr[sr->perm[i]];
1040 58 : if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1041 58 : PetscFunctionReturn(PETSC_SUCCESS);
1042 : }
1043 :
1044 : /*@
1045 : EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
1046 : internally in the subcommunicator to which the calling process belongs.
1047 :
1048 : Collective on the subcommunicator (if v is given)
1049 :
1050 : Input Parameters:
1051 : + eps - the eigenproblem solver context
1052 : - i - index of the solution
1053 :
1054 : Output Parameters:
1055 : + eig - the eigenvalue
1056 : - v - the eigenvector
1057 :
1058 : Notes:
1059 : It is allowed to pass NULL for v if the eigenvector is not required.
1060 : Otherwise, the caller must provide a valid Vec objects, i.e.,
1061 : it must be created by the calling program with EPSKrylovSchurGetSubcommInfo().
1062 :
1063 : The index i should be a value between 0 and n-1, where n is the number of
1064 : vectors in the local subinterval, see EPSKrylovSchurGetSubcommInfo().
1065 :
1066 : Level: advanced
1067 :
1068 : .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo(), EPSKrylovSchurGetSubcommMats()
1069 : @*/
1070 58 : PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1071 : {
1072 58 : PetscFunctionBegin;
1073 58 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1074 174 : if (v) PetscValidLogicalCollectiveInt(v,i,2);
1075 58 : PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1076 58 : PetscFunctionReturn(PETSC_SUCCESS);
1077 : }
1078 :
1079 2 : static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1080 : {
1081 2 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1082 :
1083 2 : PetscFunctionBegin;
1084 2 : PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1085 2 : PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1086 2 : PetscCall(EPSGetOperators(ctx->eps,A,B));
1087 2 : PetscFunctionReturn(PETSC_SUCCESS);
1088 : }
1089 :
1090 : /*@
1091 : EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1092 : internally in the subcommunicator to which the calling process belongs.
1093 :
1094 : Collective on the subcommunicator
1095 :
1096 : Input Parameter:
1097 : . eps - the eigenproblem solver context
1098 :
1099 : Output Parameters:
1100 : + A - the matrix associated with the eigensystem
1101 : - B - the second matrix in the case of generalized eigenproblems
1102 :
1103 : Notes:
1104 : This is the analog of EPSGetOperators(), but returns the matrices distributed
1105 : differently (in the subcommunicator rather than in the parent communicator).
1106 :
1107 : These matrices should not be modified by the user.
1108 :
1109 : Level: advanced
1110 :
1111 : .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo()
1112 : @*/
1113 2 : PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1114 : {
1115 2 : PetscFunctionBegin;
1116 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1117 2 : PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1118 2 : PetscFunctionReturn(PETSC_SUCCESS);
1119 : }
1120 :
1121 2 : static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1122 : {
1123 2 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1124 2 : Mat A,B=NULL,Ag,Bg=NULL;
1125 2 : PetscBool reuse=PETSC_TRUE;
1126 :
1127 2 : PetscFunctionBegin;
1128 2 : PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1129 2 : PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1130 2 : PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1131 2 : PetscCall(EPSGetOperators(ctx->eps,&A,&B));
1132 :
1133 2 : PetscCall(MatScale(A,a));
1134 2 : if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1135 2 : if (B) PetscCall(MatScale(B,b));
1136 2 : if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1137 2 : PetscCall(EPSSetOperators(ctx->eps,A,B));
1138 :
1139 : /* Update stored matrix state */
1140 2 : subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1141 2 : PetscCall(MatGetState(A,&subctx->Astate));
1142 2 : if (B) PetscCall(MatGetState(B,&subctx->Bstate));
1143 :
1144 : /* Update matrices in the parent communicator if requested by user */
1145 2 : if (globalup) {
1146 2 : if (ctx->npart>1) {
1147 2 : if (!ctx->isrow) {
1148 2 : PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1149 : reuse = PETSC_FALSE;
1150 : }
1151 2 : if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1152 2 : if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1153 2 : PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1154 2 : PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1155 2 : if (B) {
1156 2 : if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1157 2 : PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1158 2 : PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1159 : }
1160 : }
1161 2 : PetscCall(MatGetState(Ag,&ctx->Astate));
1162 2 : if (Bg) PetscCall(MatGetState(Bg,&ctx->Bstate));
1163 : }
1164 2 : PetscCall(EPSSetOperators(eps,Ag,Bg));
1165 2 : PetscFunctionReturn(PETSC_SUCCESS);
1166 : }
1167 :
1168 : /*@
1169 : EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1170 : internally in the subcommunicator to which the calling process belongs.
1171 :
1172 : Collective
1173 :
1174 : Input Parameters:
1175 : + eps - the eigenproblem solver context
1176 : . s - scalar that multiplies the existing A matrix
1177 : . a - scalar used in the axpy operation on A
1178 : . Au - matrix used in the axpy operation on A
1179 : . t - scalar that multiplies the existing B matrix
1180 : . b - scalar used in the axpy operation on B
1181 : . Bu - matrix used in the axpy operation on B
1182 : . str - structure flag
1183 : - globalup - flag indicating if global matrices must be updated
1184 :
1185 : Notes:
1186 : This function modifies the eigenproblem matrices at the subcommunicator level,
1187 : and optionally updates the global matrices in the parent communicator. The updates
1188 : are expressed as A <-- s*A + a*Au, B <-- t*B + b*Bu.
1189 :
1190 : It is possible to update one of the matrices, or both.
1191 :
1192 : The matrices Au and Bu must be equal in all subcommunicators.
1193 :
1194 : The str flag is passed to the MatAXPY() operations to perform the updates.
1195 :
1196 : If globalup is true, communication is carried out to reconstruct the updated
1197 : matrices in the parent communicator. The user must be warned that if global
1198 : matrices are not in sync with subcommunicator matrices, the errors computed
1199 : by EPSComputeError() will be wrong even if the computed solution is correct
1200 : (the synchronization may be done only once at the end).
1201 :
1202 : Level: advanced
1203 :
1204 : .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommMats()
1205 : @*/
1206 2 : PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1207 : {
1208 2 : PetscFunctionBegin;
1209 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1210 6 : PetscValidLogicalCollectiveScalar(eps,s,2);
1211 6 : PetscValidLogicalCollectiveScalar(eps,a,3);
1212 2 : if (Au) PetscValidHeaderSpecific(Au,MAT_CLASSID,4);
1213 6 : PetscValidLogicalCollectiveScalar(eps,t,5);
1214 6 : PetscValidLogicalCollectiveScalar(eps,b,6);
1215 2 : if (Bu) PetscValidHeaderSpecific(Bu,MAT_CLASSID,7);
1216 6 : PetscValidLogicalCollectiveEnum(eps,str,8);
1217 6 : PetscValidLogicalCollectiveBool(eps,globalup,9);
1218 2 : PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1219 2 : PetscFunctionReturn(PETSC_SUCCESS);
1220 : }
1221 :
1222 38 : PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1223 : {
1224 38 : EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1225 38 : Mat A,B=NULL,Ar=NULL,Br=NULL;
1226 38 : PetscMPIInt rank;
1227 38 : PetscObjectState Astate,Bstate=0;
1228 38 : PetscObjectId Aid,Bid=0;
1229 38 : STType sttype;
1230 38 : PetscInt nmat;
1231 38 : const char *prefix;
1232 38 : MPI_Comm child;
1233 :
1234 38 : PetscFunctionBegin;
1235 38 : PetscCall(EPSGetOperators(eps,&A,&B));
1236 38 : if (ctx->npart==1) {
1237 24 : if (!ctx->eps) PetscCall(EPSCreate(((PetscObject)eps)->comm,&ctx->eps));
1238 24 : PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1239 24 : PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1240 24 : PetscCall(EPSSetOperators(ctx->eps,A,B));
1241 : } else {
1242 14 : PetscCall(MatGetState(A,&Astate));
1243 14 : PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1244 14 : if (B) {
1245 14 : PetscCall(MatGetState(B,&Bstate));
1246 14 : PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1247 : }
1248 14 : if (!ctx->subc) {
1249 : /* Create context for subcommunicators */
1250 5 : PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1251 5 : PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1252 5 : PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1253 5 : PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1254 :
1255 : /* Duplicate matrices */
1256 5 : PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1257 5 : ctx->Astate = Astate;
1258 5 : ctx->Aid = Aid;
1259 5 : PetscCall(MatPropagateSymmetryOptions(A,Ar));
1260 5 : if (B) {
1261 5 : PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1262 5 : ctx->Bstate = Bstate;
1263 5 : ctx->Bid = Bid;
1264 5 : PetscCall(MatPropagateSymmetryOptions(B,Br));
1265 : }
1266 : } else {
1267 9 : PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1268 9 : if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1269 0 : PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1270 0 : if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
1271 0 : PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1272 0 : ctx->Astate = Astate;
1273 0 : ctx->Aid = Aid;
1274 0 : PetscCall(MatPropagateSymmetryOptions(A,Ar));
1275 0 : if (B) {
1276 0 : PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1277 0 : ctx->Bstate = Bstate;
1278 0 : ctx->Bid = Bid;
1279 0 : PetscCall(MatPropagateSymmetryOptions(B,Br));
1280 : }
1281 0 : PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1282 0 : PetscCall(MatDestroy(&Ar));
1283 0 : PetscCall(MatDestroy(&Br));
1284 : }
1285 : }
1286 :
1287 : /* Create auxiliary EPS */
1288 14 : if (!ctx->eps) {
1289 5 : PetscCall(EPSCreate(child,&ctx->eps));
1290 5 : PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1291 5 : PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1292 5 : PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1293 5 : PetscCall(MatDestroy(&Ar));
1294 5 : PetscCall(MatDestroy(&Br));
1295 : }
1296 : /* Create subcommunicator grouping processes with same rank */
1297 14 : if (!ctx->commset) {
1298 5 : PetscCallMPI(MPI_Comm_rank(child,&rank));
1299 5 : PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1300 5 : ctx->commset = PETSC_TRUE;
1301 : }
1302 : }
1303 38 : PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1304 38 : PetscCall(STGetType(eps->st,&sttype));
1305 38 : PetscCall(STSetType(ctx->eps->st,sttype));
1306 :
1307 38 : ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1308 38 : ctx_local->npart = ctx->npart;
1309 38 : ctx_local->global = PETSC_FALSE;
1310 38 : ctx_local->eps = eps;
1311 38 : ctx_local->subc = ctx->subc;
1312 38 : ctx_local->commrank = ctx->commrank;
1313 38 : *childeps = ctx->eps;
1314 38 : PetscFunctionReturn(PETSC_SUCCESS);
1315 : }
1316 :
1317 20 : static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1318 : {
1319 20 : EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1320 20 : ST st;
1321 20 : PetscBool isfilt;
1322 :
1323 20 : PetscFunctionBegin;
1324 20 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1325 20 : PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1326 20 : PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1327 20 : PetscCall(EPSGetST(ctx->eps,&st));
1328 20 : PetscCall(STGetOperator(st,NULL));
1329 20 : PetscCall(STGetKSP(st,ksp));
1330 20 : PetscFunctionReturn(PETSC_SUCCESS);
1331 : }
1332 :
1333 : /*@
1334 : EPSKrylovSchurGetKSP - Retrieve the linear solver object associated with the
1335 : internal EPS object in case of doing spectrum slicing for a computational interval.
1336 :
1337 : Collective
1338 :
1339 : Input Parameter:
1340 : . eps - the eigenproblem solver context
1341 :
1342 : Output Parameter:
1343 : . ksp - the internal KSP object
1344 :
1345 : Notes:
1346 : When invoked to compute all eigenvalues in an interval with spectrum
1347 : slicing, EPSKRYLOVSCHUR creates another EPS object internally that is
1348 : used to compute eigenvalues by chunks near selected shifts. This function
1349 : allows access to the KSP object associated to this internal EPS object.
1350 :
1351 : This function is only available for spectrum slicing runs. In case of
1352 : having more than one partition, the returned KSP will be different
1353 : in MPI processes belonging to different partitions. Hence, if required,
1354 : EPSKrylovSchurSetPartitions() must be called BEFORE this function.
1355 :
1356 : Level: advanced
1357 :
1358 : .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions()
1359 : @*/
1360 6 : PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1361 : {
1362 6 : PetscFunctionBegin;
1363 6 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1364 6 : PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1365 6 : PetscFunctionReturn(PETSC_SUCCESS);
1366 : }
1367 :
1368 20 : static PetscErrorCode EPSKrylovSchurSetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType bse)
1369 : {
1370 20 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1371 :
1372 20 : PetscFunctionBegin;
1373 20 : switch (bse) {
1374 20 : case EPS_KRYLOVSCHUR_BSE_SHAO:
1375 : case EPS_KRYLOVSCHUR_BSE_GRUNING:
1376 : case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
1377 20 : if (ctx->bse != bse) {
1378 13 : ctx->bse = bse;
1379 13 : eps->state = EPS_STATE_INITIAL;
1380 : }
1381 20 : break;
1382 0 : default:
1383 0 : SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid BSE type");
1384 : }
1385 20 : PetscFunctionReturn(PETSC_SUCCESS);
1386 : }
1387 :
1388 : /*@
1389 : EPSKrylovSchurSetBSEType - Sets the method to be used for BSE structured eigenproblems
1390 : in the Krylov-Schur solver.
1391 :
1392 : Logically Collective
1393 :
1394 : Input Parameters:
1395 : + eps - the eigenproblem solver context
1396 : - bse - the BSE method
1397 :
1398 : Options Database Key:
1399 : . -eps_krylovschur_bse_type - Sets the BSE type (either 'shao', 'gruning', or 'projectedbse')
1400 :
1401 : Level: advanced
1402 :
1403 : .seealso: EPSKrylovSchurGetBSEType(), EPSKrylovSchurBSEType, MatCreateBSE()
1404 : @*/
1405 20 : PetscErrorCode EPSKrylovSchurSetBSEType(EPS eps,EPSKrylovSchurBSEType bse)
1406 : {
1407 20 : PetscFunctionBegin;
1408 20 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1409 60 : PetscValidLogicalCollectiveEnum(eps,bse,2);
1410 20 : PetscTryMethod(eps,"EPSKrylovSchurSetBSEType_C",(EPS,EPSKrylovSchurBSEType),(eps,bse));
1411 20 : PetscFunctionReturn(PETSC_SUCCESS);
1412 : }
1413 :
1414 0 : static PetscErrorCode EPSKrylovSchurGetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType *bse)
1415 : {
1416 0 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1417 :
1418 0 : PetscFunctionBegin;
1419 0 : *bse = ctx->bse;
1420 0 : PetscFunctionReturn(PETSC_SUCCESS);
1421 : }
1422 :
1423 : /*@
1424 : EPSKrylovSchurGetBSEType - Gets the method used for BSE structured eigenproblems
1425 : in the Krylov-Schur solver.
1426 :
1427 : Not Collective
1428 :
1429 : Input Parameter:
1430 : . eps - the eigenproblem solver context
1431 :
1432 : Output Parameter:
1433 : . bse - the BSE method
1434 :
1435 : Level: advanced
1436 :
1437 : .seealso: EPSKrylovSchurSetBSEType(), EPSKrylovSchurBSEType, MatCreateBSE()
1438 : @*/
1439 0 : PetscErrorCode EPSKrylovSchurGetBSEType(EPS eps,EPSKrylovSchurBSEType *bse)
1440 : {
1441 0 : PetscFunctionBegin;
1442 0 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
1443 0 : PetscAssertPointer(bse,2);
1444 0 : PetscUseMethod(eps,"EPSKrylovSchurGetBSEType_C",(EPS,EPSKrylovSchurBSEType*),(eps,bse));
1445 0 : PetscFunctionReturn(PETSC_SUCCESS);
1446 : }
1447 :
1448 332 : static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems *PetscOptionsObject)
1449 : {
1450 332 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1451 332 : PetscBool flg,lock,b,f1,f2,f3,isfilt;
1452 332 : PetscReal keep;
1453 332 : PetscInt i,j,k;
1454 332 : KSP ksp;
1455 332 : EPSKrylovSchurBSEType bse;
1456 :
1457 332 : PetscFunctionBegin;
1458 332 : PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");
1459 :
1460 332 : PetscCall(PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg));
1461 332 : if (flg) PetscCall(EPSKrylovSchurSetRestart(eps,keep));
1462 :
1463 332 : PetscCall(PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg));
1464 332 : if (flg) PetscCall(EPSKrylovSchurSetLocking(eps,lock));
1465 :
1466 332 : i = ctx->npart;
1467 332 : PetscCall(PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg));
1468 332 : if (flg) PetscCall(EPSKrylovSchurSetPartitions(eps,i));
1469 :
1470 332 : b = ctx->detect;
1471 332 : PetscCall(PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg));
1472 332 : if (flg) PetscCall(EPSKrylovSchurSetDetectZeros(eps,b));
1473 :
1474 332 : i = 1;
1475 332 : j = k = PETSC_DECIDE;
1476 332 : PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1477 332 : PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1478 332 : PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1479 332 : if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));
1480 :
1481 332 : PetscCall(PetscOptionsEnum("-eps_krylovschur_bse_type","Method for BSE structured eigenproblems","EPSKrylovSchurSetBSEType",EPSKrylovSchurBSETypes,(PetscEnum)ctx->bse,(PetscEnum*)&bse,&flg));
1482 332 : if (flg) PetscCall(EPSKrylovSchurSetBSEType(eps,bse));
1483 :
1484 332 : PetscOptionsHeadEnd();
1485 :
1486 : /* set options of child KSP in spectrum slicing */
1487 332 : if (eps->which==EPS_ALL) {
1488 19 : if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1489 19 : PetscCall(EPSSetDefaultST(eps));
1490 19 : PetscCall(STSetFromOptions(eps->st)); /* need to advance this to check ST type */
1491 19 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1492 19 : if (!isfilt) {
1493 14 : PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1494 14 : PetscCall(KSPSetFromOptions(ksp));
1495 : }
1496 : }
1497 332 : PetscFunctionReturn(PETSC_SUCCESS);
1498 : }
1499 :
1500 4 : static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1501 : {
1502 4 : EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1503 4 : PetscBool isascii,isfilt;
1504 4 : KSP ksp;
1505 4 : PetscViewer sviewer;
1506 :
1507 4 : PetscFunctionBegin;
1508 4 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1509 4 : if (isascii) {
1510 4 : PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1511 4 : PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1512 4 : if (eps->problem_type==EPS_BSE) PetscCall(PetscViewerASCIIPrintf(viewer," BSE method: %s\n",EPSKrylovSchurBSETypes[ctx->bse]));
1513 4 : if (eps->which==EPS_ALL) {
1514 0 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1515 0 : if (isfilt) PetscCall(PetscViewerASCIIPrintf(viewer," using filtering to extract all eigenvalues in an interval\n"));
1516 : else {
1517 0 : PetscCall(PetscViewerASCIIPrintf(viewer," doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd));
1518 0 : if (ctx->npart>1) {
1519 0 : PetscCall(PetscViewerASCIIPrintf(viewer," multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart));
1520 0 : if (ctx->detect) PetscCall(PetscViewerASCIIPrintf(viewer," detecting zeros when factorizing at subinterval boundaries\n"));
1521 : }
1522 : /* view child KSP */
1523 0 : PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1524 0 : PetscCall(PetscViewerASCIIPushTab(viewer));
1525 0 : if (ctx->npart>1 && ctx->subc) {
1526 0 : PetscCall(PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer));
1527 0 : if (!ctx->subc->color) PetscCall(KSPView(ksp,sviewer));
1528 0 : PetscCall(PetscViewerFlush(sviewer));
1529 0 : PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer));
1530 : /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1531 0 : PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1532 0 : } else PetscCall(KSPView(ksp,viewer));
1533 0 : PetscCall(PetscViewerASCIIPopTab(viewer));
1534 : }
1535 : }
1536 : }
1537 4 : PetscFunctionReturn(PETSC_SUCCESS);
1538 : }
1539 :
1540 365 : static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1541 : {
1542 365 : PetscBool isfilt;
1543 :
1544 365 : PetscFunctionBegin;
1545 365 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1546 365 : if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1547 365 : PetscCall(PetscFree(eps->data));
1548 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1549 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1550 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1551 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1552 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1553 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1554 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1555 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1556 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1557 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1558 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1559 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1560 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1561 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1562 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1563 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1564 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1565 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1566 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",NULL));
1567 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",NULL));
1568 365 : PetscFunctionReturn(PETSC_SUCCESS);
1569 : }
1570 :
1571 520 : static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1572 : {
1573 520 : PetscBool isfilt;
1574 :
1575 520 : PetscFunctionBegin;
1576 520 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1577 520 : if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1578 520 : PetscFunctionReturn(PETSC_SUCCESS);
1579 : }
1580 :
1581 877 : static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1582 : {
1583 877 : PetscFunctionBegin;
1584 877 : if (eps->which==EPS_ALL) {
1585 82 : if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1586 : }
1587 877 : PetscFunctionReturn(PETSC_SUCCESS);
1588 : }
1589 :
1590 365 : SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1591 : {
1592 365 : EPS_KRYLOVSCHUR *ctx;
1593 :
1594 365 : PetscFunctionBegin;
1595 365 : PetscCall(PetscNew(&ctx));
1596 365 : eps->data = (void*)ctx;
1597 365 : ctx->lock = PETSC_TRUE;
1598 365 : ctx->nev = 1;
1599 365 : ctx->ncv = PETSC_DETERMINE;
1600 365 : ctx->mpd = PETSC_DETERMINE;
1601 365 : ctx->npart = 1;
1602 365 : ctx->detect = PETSC_FALSE;
1603 365 : ctx->global = PETSC_TRUE;
1604 :
1605 365 : eps->useds = PETSC_TRUE;
1606 :
1607 : /* solve and computevectors determined at setup */
1608 365 : eps->ops->setup = EPSSetUp_KrylovSchur;
1609 365 : eps->ops->setupsort = EPSSetUpSort_KrylovSchur;
1610 365 : eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1611 365 : eps->ops->destroy = EPSDestroy_KrylovSchur;
1612 365 : eps->ops->reset = EPSReset_KrylovSchur;
1613 365 : eps->ops->view = EPSView_KrylovSchur;
1614 365 : eps->ops->backtransform = EPSBackTransform_Default;
1615 365 : eps->ops->setdefaultst = EPSSetDefaultST_KrylovSchur;
1616 :
1617 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1618 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1619 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1620 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1621 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1622 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1623 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1624 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1625 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1626 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1627 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1628 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1629 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1630 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1631 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1632 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1633 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1634 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1635 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",EPSKrylovSchurSetBSEType_KrylovSchur));
1636 365 : PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",EPSKrylovSchurGetBSEType_KrylovSchur));
1637 365 : PetscFunctionReturn(PETSC_SUCCESS);
1638 : }
|
