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