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