Actual source code: krylovschur.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
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"
13: Method: Krylov-Schur
15: Algorithm:
17: Single-vector Krylov-Schur method for non-symmetric problems,
18: including harmonic extraction.
20: References:
22: [1] "Krylov-Schur Methods in SLEPc", SLEPc Technical Report STR-7,
23: available at https://slepc.upv.es.
25: [2] G.W. Stewart, "A Krylov-Schur Algorithm for Large Eigenproblems",
26: SIAM J. Matrix Anal. App. 23(3):601-614, 2001.
28: [3] "Practical Implementation of Harmonic Krylov-Schur", SLEPc Technical
29: Report STR-9, available at https://slepc.upv.es.
30: */
32: #include <slepc/private/epsimpl.h>
33: #include "krylovschur.h"
35: PetscErrorCode EPSGetArbitraryValues(EPS eps,PetscScalar *rr,PetscScalar *ri)
36: {
37: PetscInt i,newi,ld,n,l;
38: Vec xr=eps->work[0],xi=eps->work[1];
39: PetscScalar re,im,*Zr,*Zi,*X;
41: PetscFunctionBegin;
42: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
43: PetscCall(DSGetDimensions(eps->ds,&n,&l,NULL,NULL));
44: for (i=l;i<n;i++) {
45: re = eps->eigr[i];
46: im = eps->eigi[i];
47: PetscCall(STBackTransform(eps->st,1,&re,&im));
48: newi = i;
49: PetscCall(DSVectors(eps->ds,DS_MAT_X,&newi,NULL));
50: PetscCall(DSGetArray(eps->ds,DS_MAT_X,&X));
51: Zr = X+i*ld;
52: if (newi==i+1) Zi = X+newi*ld;
53: else Zi = NULL;
54: PetscCall(EPSComputeRitzVector(eps,Zr,Zi,eps->V,xr,xi));
55: PetscCall(DSRestoreArray(eps->ds,DS_MAT_X,&X));
56: PetscCall((*eps->arbitrary)(re,im,xr,xi,rr+i,ri+i,eps->arbitraryctx));
57: }
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode EPSSetUp_KrylovSchur_Filter(EPS eps)
62: {
63: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
64: PetscBool estimaterange=PETSC_TRUE,same;
65: PetscReal rleft,rright;
66: Mat A;
68: PetscFunctionBegin;
69: EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
70: EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
71: 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: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
73: PetscCall(STFilterSetInterval(eps->st,eps->inta,eps->intb));
74: if (!ctx->estimatedrange) {
75: PetscCall(STFilterGetRange(eps->st,&rleft,&rright));
76: estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
77: }
78: PetscCall(STGetMatrix(eps->st,0,&A));
79: PetscCall(MatStateCompareUpdate(A,&ctx->fstate,&same));
80: if (estimaterange && !same) { /* user did not set a range, and is not available from previous run */
81: PetscCall(MatEstimateSpectralRange_EPS(A,&rleft,&rright));
82: PetscCall(PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright));
83: PetscCall(STFilterSetRange(eps->st,rleft,rright));
84: ctx->estimatedrange = PETSC_TRUE;
85: }
86: PetscCall(EPSSetStoppingTest(eps,EPS_STOP_THRESHOLD));
87: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
88: PetscCheck(eps->nev==0 || 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: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: static PetscErrorCode EPSSetUp_KrylovSchur(EPS eps)
94: {
95: PetscReal eta;
96: PetscBool isfilt=PETSC_FALSE;
97: BVOrthogType otype;
98: BVOrthogBlockType obtype;
99: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
100: enum { EPS_KS_DEFAULT,EPS_KS_SYMM,EPS_KS_SLICE,EPS_KS_FILTER,EPS_KS_INDEF,EPS_KS_TWOSIDED } variant;
102: PetscFunctionBegin;
103: if (eps->which==EPS_ALL) { /* default values in case of spectrum slicing or polynomial filter */
104: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
105: if (isfilt) PetscCall(EPSSetUp_KrylovSchur_Filter(eps));
106: else PetscCall(EPSSetUp_KrylovSchur_Slice(eps));
107: } else if (eps->isstructured) {
108: if (eps->problem_type==EPS_BSE) PetscCall(EPSSetUp_KrylovSchur_BSE(eps));
109: else if (eps->problem_type==EPS_HAMILT) {
110: PetscCheck(!PetscDefined(USE_COMPLEX),PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The Hamiltonian Krylov-Schur eigensolver is not yet implemented for complex scalars");
111: PetscCall(EPSSetUp_KrylovSchur_Hamilt(eps));
112: } else if (eps->problem_type==EPS_LREP) {
113: PetscCheck(!PetscDefined(USE_COMPLEX),PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The LREP Krylov-Schur eigensolver does not support complex scalars, use BSE instead");
114: PetscCall(EPSSetUp_KrylovSchur_LREP(eps));
115: } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unknown matrix structure");
116: PetscFunctionReturn(PETSC_SUCCESS);
117: } else {
118: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
119: 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");
120: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);
121: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
122: }
123: PetscCheck(ctx->lock || eps->mpd>=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
125: EPSCheckDefiniteCondition(eps,eps->arbitrary," with arbitrary selection of eigenpairs");
127: PetscCheck(eps->extraction==EPS_RITZ || eps->extraction==EPS_HARMONIC,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
129: if (!ctx->keep) ctx->keep = 0.5;
131: PetscCall(EPSAllocateSolution(eps,1));
132: PetscCall(EPS_SetInnerProduct(eps));
133: if (eps->arbitrary) PetscCall(EPSSetWorkVecs(eps,2));
134: else if (eps->ishermitian && !eps->ispositive) PetscCall(EPSSetWorkVecs(eps,1));
136: /* dispatch solve method */
137: if (eps->ishermitian) {
138: if (eps->which==EPS_ALL) {
139: EPSCheckDefiniteCondition(eps,eps->which==EPS_ALL," with spectrum slicing");
140: variant = isfilt? EPS_KS_FILTER: EPS_KS_SLICE;
141: } else if (eps->isgeneralized && !eps->ispositive) {
142: variant = EPS_KS_INDEF;
143: } else {
144: switch (eps->extraction) {
145: case EPS_RITZ: variant = EPS_KS_SYMM; break;
146: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
147: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
148: }
149: }
150: } else if (eps->twosided) {
151: variant = EPS_KS_TWOSIDED;
152: } else {
153: switch (eps->extraction) {
154: case EPS_RITZ: variant = EPS_KS_DEFAULT; break;
155: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
156: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
157: }
158: }
159: switch (variant) {
160: case EPS_KS_DEFAULT:
161: eps->ops->solve = EPSSolve_KrylovSchur_Default;
162: eps->ops->computevectors = EPSComputeVectors_Schur;
163: PetscCall(DSSetType(eps->ds,DSNHEP));
164: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
165: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
166: break;
167: case EPS_KS_SYMM:
168: case EPS_KS_FILTER:
169: eps->ops->solve = EPSSolve_KrylovSchur_Default;
170: eps->ops->computevectors = EPSComputeVectors_Hermitian;
171: PetscCall(DSSetType(eps->ds,DSHEP));
172: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
173: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
174: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
175: break;
176: case EPS_KS_SLICE:
177: eps->ops->solve = EPSSolve_KrylovSchur_Slice;
178: eps->ops->computevectors = EPSComputeVectors_Slice;
179: break;
180: case EPS_KS_INDEF:
181: eps->ops->solve = EPSSolve_KrylovSchur_Indefinite;
182: eps->ops->computevectors = EPSComputeVectors_Indefinite;
183: PetscCall(DSSetType(eps->ds,DSGHIEP));
184: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
185: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
186: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
187: /* force reorthogonalization for pseudo-Lanczos */
188: PetscCall(BVGetOrthogonalization(eps->V,&otype,NULL,&eta,&obtype));
189: PetscCall(BVSetOrthogonalization(eps->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
190: break;
191: case EPS_KS_TWOSIDED:
192: eps->ops->solve = EPSSolve_KrylovSchur_TwoSided;
193: eps->ops->computevectors = EPSComputeVectors_Schur;
194: PetscCall(DSSetType(eps->ds,DSNHEPTS));
195: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
196: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
197: break;
198: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
199: }
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: static PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
204: {
205: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
206: SlepcSC sc;
207: PetscBool isfilt;
209: PetscFunctionBegin;
210: PetscCall(EPSSetUpSort_Default(eps));
211: if (eps->which==EPS_ALL) {
212: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
213: if (isfilt) {
214: PetscCall(DSGetSlepcSC(eps->ds,&sc));
215: sc->rg = NULL;
216: sc->comparison = SlepcCompareLargestReal;
217: sc->comparisonctx = NULL;
218: sc->map = NULL;
219: sc->mapobj = NULL;
220: } else {
221: if (!ctx->global && ctx->sr->numEigs>0) {
222: PetscCall(DSGetSlepcSC(eps->ds,&sc));
223: sc->rg = NULL;
224: sc->comparison = SlepcCompareLargestMagnitude;
225: sc->comparisonctx = NULL;
226: sc->map = NULL;
227: sc->mapobj = NULL;
228: }
229: }
230: }
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: PetscErrorCode EPSSolve_KrylovSchur_Default(EPS eps)
235: {
236: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
237: PetscInt i,j,*pj,k,l,nv,ld,nconv;
238: Mat U,Op,H,T;
239: PetscScalar *g;
240: PetscReal beta,gamma=1.0;
241: PetscBool breakdown,harmonic,hermitian;
243: PetscFunctionBegin;
244: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
245: harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
246: hermitian = (eps->ishermitian && !harmonic)?PETSC_TRUE:PETSC_FALSE;
247: if (harmonic) PetscCall(PetscMalloc1(ld,&g));
248: if (eps->arbitrary) pj = &j;
249: else pj = NULL;
251: /* Get the starting Arnoldi vector */
252: PetscCall(EPSGetStartVector(eps,0,NULL));
253: l = 0;
255: /* Restart loop */
256: while (eps->reason == EPS_CONVERGED_ITERATING) {
257: eps->its++;
259: /* Compute an nv-step Arnoldi factorization */
260: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
261: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
262: PetscCall(STGetOperator(eps->st,&Op));
263: if (hermitian) {
264: PetscCall(DSGetMat(eps->ds,DS_MAT_T,&T));
265: PetscCall(BVMatLanczos(eps->V,Op,T,eps->nconv+l,&nv,&beta,&breakdown));
266: PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&T));
267: } else {
268: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
269: PetscCall(BVMatArnoldi(eps->V,Op,H,eps->nconv+l,&nv,&beta,&breakdown));
270: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
271: }
272: PetscCall(STRestoreOperator(eps->st,&Op));
273: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
274: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
275: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
277: /* Compute translation of Krylov decomposition if harmonic extraction used */
278: if (PetscUnlikely(harmonic)) PetscCall(DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,g,&gamma));
280: /* Solve projected problem */
281: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
282: if (PetscUnlikely(eps->arbitrary)) {
283: PetscCall(EPSGetArbitraryValues(eps,eps->rr,eps->ri));
284: j=1;
285: }
286: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj));
287: PetscCall(DSUpdateExtraRow(eps->ds));
288: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
290: /* Check convergence */
291: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
292: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,k,nv);
293: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
294: nconv = k;
296: /* Update l */
297: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
298: else {
299: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
300: if (!hermitian) PetscCall(DSGetTruncateSize(eps->ds,k,nv,&l));
301: }
302: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
303: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
305: if (eps->reason == EPS_CONVERGED_ITERATING) {
306: if (PetscUnlikely(breakdown || k==nv)) {
307: /* Start a new Arnoldi factorization */
308: PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
309: if (k<eps->nev) {
310: PetscCall(EPSGetStartVector(eps,k,&breakdown));
311: if (breakdown) {
312: eps->reason = EPS_DIVERGED_BREAKDOWN;
313: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
314: }
315: }
316: } else {
317: /* Undo translation of Krylov decomposition */
318: if (PetscUnlikely(harmonic)) {
319: PetscCall(DSSetDimensions(eps->ds,nv,k,l));
320: PetscCall(DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma));
321: /* gamma u^ = u - U*g~ */
322: PetscCall(BVSetActiveColumns(eps->V,0,nv));
323: PetscCall(BVMultColumn(eps->V,-1.0,1.0,nv,g));
324: PetscCall(BVScaleColumn(eps->V,nv,1.0/gamma));
325: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
326: PetscCall(DSSetDimensions(eps->ds,nv,k,nv));
327: }
328: /* Prepare the Rayleigh quotient for restart */
329: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
330: }
331: }
332: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
333: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
334: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
335: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
337: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
338: PetscCall(BVCopyColumn(eps->V,nv,k+l)); /* copy restart vector from the last column */
339: if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
340: eps->ncv = eps->mpd+k;
341: PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
342: for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
343: PetscCall(DSReallocate(eps->ds,eps->ncv+1));
344: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
345: }
346: }
348: eps->nconv = k;
349: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
350: }
352: if (harmonic) PetscCall(PetscFree(g));
353: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
358: {
359: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
361: PetscFunctionBegin;
362: if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
363: else {
364: 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);
365: ctx->keep = keep;
366: }
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: /*@
371: EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
372: method, in particular the proportion of basis vectors that must be kept
373: after restart.
375: Logically Collective
377: Input Parameters:
378: + eps - the linear eigensolver context
379: - keep - the number of vectors to be kept at restart
381: Options Database Key:
382: . -eps_krylovschur_restart keep - sets the restart parameter
384: Notes:
385: Allowed values are in the range [0.1,0.9]. The default is 0.5, which means
386: that at restart the current subspace is compressed into another subspace
387: with a reduction of 50% in size.
389: Implementation details of Krylov-Schur in SLEPc can be found in {cite:p}`Her07b`.
391: Level: advanced
393: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetRestart()`
394: @*/
395: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
396: {
397: PetscFunctionBegin;
400: PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
405: {
406: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
408: PetscFunctionBegin;
409: *keep = ctx->keep;
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /*@
414: EPSKrylovSchurGetRestart - Gets the restart parameter used in the
415: Krylov-Schur method.
417: Not Collective
419: Input Parameter:
420: . eps - the linear eigensolver context
422: Output Parameter:
423: . keep - the restart parameter
425: Level: advanced
427: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetRestart()`
428: @*/
429: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
430: {
431: PetscFunctionBegin;
433: PetscAssertPointer(keep,2);
434: PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
439: {
440: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
442: PetscFunctionBegin;
443: ctx->lock = lock;
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: /*@
448: EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
449: the Krylov-Schur method.
451: Logically Collective
453: Input Parameters:
454: + eps - the linear eigensolver context
455: - lock - true if the locking variant must be selected
457: Options Database Key:
458: . -eps_krylovschur_locking (true|false) - sets the locking flag
460: Notes:
461: The default is to lock converged eigenpairs when the method restarts.
462: This behavior can be changed so that all directions are kept in the
463: working subspace even if already converged to working accuracy (the
464: non-locking variant).
466: Implementation details of Krylov-Schur in SLEPc can be found in {cite:p}`Her07b`.
468: Level: advanced
470: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetLocking()`
471: @*/
472: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
473: {
474: PetscFunctionBegin;
477: PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
482: {
483: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
485: PetscFunctionBegin;
486: *lock = ctx->lock;
487: PetscFunctionReturn(PETSC_SUCCESS);
488: }
490: /*@
491: EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
492: method.
494: Not Collective
496: Input Parameter:
497: . eps - the linear eigensolver context
499: Output Parameter:
500: . lock - the locking flag
502: Level: advanced
504: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetLocking()`
505: @*/
506: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
507: {
508: PetscFunctionBegin;
510: PetscAssertPointer(lock,2);
511: PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
512: PetscFunctionReturn(PETSC_SUCCESS);
513: }
515: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
516: {
517: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
518: PetscMPIInt size;
519: PetscInt newnpart;
521: PetscFunctionBegin;
522: if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
523: newnpart = 1;
524: } else {
525: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
526: PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
527: newnpart = npart;
528: }
529: if (ctx->npart!=newnpart) {
530: if (ctx->npart>1) {
531: PetscCall(PetscSubcommDestroy(&ctx->subc));
532: if (ctx->commset) {
533: PetscCallMPI(MPI_Comm_free(&ctx->commrank));
534: ctx->commset = PETSC_FALSE;
535: }
536: }
537: PetscCall(EPSDestroy(&ctx->eps));
538: ctx->npart = newnpart;
539: eps->state = EPS_STATE_INITIAL;
540: }
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /*@
545: EPSKrylovSchurSetPartitions - Sets the number of partitions for the
546: case of doing spectrum slicing for a computational interval with the
547: communicator split in several sub-communicators.
549: Logically Collective
551: Input Parameters:
552: + eps - the linear eigensolver context
553: - npart - number of partitions
555: Options Database Key:
556: . -eps_krylovschur_partitions npart - sets the number of partitions
558: Notes:
559: This call makes sense only for spectrum slicing runs, that is, when
560: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
561: See more details in section [](#sec:slice).
563: By default, `npart`=1 so all processes in the communicator participate in
564: the processing of the whole interval. If `npart`>1 then the interval is
565: divided into `npart` subintervals, each of them being processed by a
566: subset of processes.
568: The interval is split proportionally unless the separation points are
569: specified with `EPSKrylovSchurSetSubintervals()`.
571: Level: advanced
573: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetSubintervals()`, `EPSSetInterval()`
574: @*/
575: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
576: {
577: PetscFunctionBegin;
580: PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
585: {
586: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
588: PetscFunctionBegin;
589: *npart = ctx->npart;
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: /*@
594: EPSKrylovSchurGetPartitions - Gets the number of partitions of the
595: communicator in case of spectrum slicing.
597: Not Collective
599: Input Parameter:
600: . eps - the linear eigensolver context
602: Output Parameter:
603: . npart - number of partitions
605: Level: advanced
607: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetPartitions()`
608: @*/
609: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
610: {
611: PetscFunctionBegin;
613: PetscAssertPointer(npart,2);
614: PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
615: PetscFunctionReturn(PETSC_SUCCESS);
616: }
618: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
619: {
620: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
622: PetscFunctionBegin;
623: ctx->detect = detect;
624: eps->state = EPS_STATE_INITIAL;
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }
628: /*@
629: EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
630: zeros during the factorizations throughout the spectrum slicing computation.
632: Logically Collective
634: Input Parameters:
635: + eps - the linear eigensolver context
636: - detect - check for zeros
638: Options Database Key:
639: . -eps_krylovschur_detect_zeros (true|false) - check for zeros
641: Notes:
642: This flag makes sense only for spectrum slicing runs, that is, when
643: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
644: See more details in section [](#sec:slice).
646: A zero in the factorization indicates that a shift coincides with an eigenvalue.
648: This flag is turned off by default, and may be necessary in some cases,
649: especially when several partitions are being used. This feature currently
650: requires an external package for factorizations with support for zero
651: detection, e.g. MUMPS.
653: Level: advanced
655: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetPartitions()`, `EPSSetInterval()`
656: @*/
657: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
658: {
659: PetscFunctionBegin;
662: PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
667: {
668: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
670: PetscFunctionBegin;
671: *detect = ctx->detect;
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: /*@
676: EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
677: in spectrum slicing.
679: Not Collective
681: Input Parameter:
682: . eps - the linear eigensolver context
684: Output Parameter:
685: . detect - whether zeros detection is enforced during factorizations
687: Level: advanced
689: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetDetectZeros()`
690: @*/
691: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
692: {
693: PetscFunctionBegin;
695: PetscAssertPointer(detect,2);
696: PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
697: PetscFunctionReturn(PETSC_SUCCESS);
698: }
700: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
701: {
702: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
704: PetscFunctionBegin;
705: if (nev != PETSC_CURRENT) {
706: PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
707: ctx->nev = nev;
708: }
709: if (ncv == PETSC_DETERMINE) {
710: ctx->ncv = PETSC_DETERMINE;
711: } else if (ncv != PETSC_CURRENT) {
712: PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
713: ctx->ncv = ncv;
714: }
715: if (mpd == PETSC_DETERMINE) {
716: ctx->mpd = PETSC_DETERMINE;
717: } else if (mpd != PETSC_CURRENT) {
718: PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
719: ctx->mpd = mpd;
720: }
721: eps->state = EPS_STATE_INITIAL;
722: PetscFunctionReturn(PETSC_SUCCESS);
723: }
725: /*@
726: EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
727: step in case of doing spectrum slicing for a computational interval.
729: Logically Collective
731: Input Parameters:
732: + eps - the linear eigensolver context
733: . nev - number of eigenvalues to compute
734: . ncv - the maximum dimension of the subspace to be used by the subsolve
735: - mpd - the maximum dimension allowed for the projected problem
737: Options Database Keys:
738: + -eps_krylovschur_nev nev - sets the number of eigenvalues
739: . -eps_krylovschur_ncv ncv - sets the dimension of the subspace
740: - -eps_krylovschur_mpd mpd - sets the maximum projected dimension
742: Notes:
743: These parameters are relevant only for spectrum slicing runs, that is, when
744: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
745: See more details in section [](#sec:slice).
747: The meaning of the parameters is the same as in `EPSSetDimensions()`, but
748: the ones here apply to every subsolve done by the child `EPS` object.
750: Use `PETSC_DETERMINE` for `ncv` and `mpd` to assign a default value. For any
751: of the arguments, use `PETSC_CURRENT` to preserve the current value.
753: Level: advanced
755: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetDimensions()`, `EPSSetDimensions()`, `EPSSetInterval()`
756: @*/
757: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
758: {
759: PetscFunctionBegin;
764: PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
769: {
770: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
772: PetscFunctionBegin;
773: if (nev) *nev = ctx->nev;
774: if (ncv) *ncv = ctx->ncv;
775: if (mpd) *mpd = ctx->mpd;
776: PetscFunctionReturn(PETSC_SUCCESS);
777: }
779: /*@
780: EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
781: step in case of doing spectrum slicing for a computational interval.
783: Not Collective
785: Input Parameter:
786: . eps - the linear eigensolver context
788: Output Parameters:
789: + nev - number of eigenvalues to compute
790: . ncv - the maximum dimension of the subspace to be used by the subsolve
791: - mpd - the maximum dimension allowed for the projected problem
793: Level: advanced
795: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetDimensions()`
796: @*/
797: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
798: {
799: PetscFunctionBegin;
801: PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
802: PetscFunctionReturn(PETSC_SUCCESS);
803: }
805: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
806: {
807: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
808: PetscInt i;
810: PetscFunctionBegin;
811: 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()");
812: 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");
813: PetscCall(PetscFree(ctx->subintervals));
814: PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
815: for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
816: ctx->subintset = PETSC_TRUE;
817: eps->state = EPS_STATE_INITIAL;
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: /*@
822: EPSKrylovSchurSetSubintervals - Sets the points that delimit the
823: subintervals to be used in spectrum slicing with several partitions.
825: Logically Collective
827: Input Parameters:
828: + eps - the linear eigensolver context
829: - subint - array of real values specifying subintervals
831: Notes:
832: This function is relevant only for spectrum slicing runs, that is, when
833: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
834: See more details in section [](#sec:slice).
836: It must be called after `EPSKrylovSchurSetPartitions()`. For `npart`
837: partitions, the argument `subint` must contain `npart+1` real values sorted in
838: ascending order, `subint_0, subint_1, ..., subint_npart`, where the first
839: and last values must coincide with the interval endpoints set with
840: `EPSSetInterval()`.
842: The subintervals are then defined by two consecutive points `[subint_0,subint_1]`,
843: `[subint_1,subint_2]`, and so on.
845: Level: advanced
847: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubintervals()`, `EPSSetInterval()`
848: @*/
849: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal subint[])
850: {
851: PetscFunctionBegin;
853: PetscAssertPointer(subint,2);
854: PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
855: PetscFunctionReturn(PETSC_SUCCESS);
856: }
858: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal *subint[])
859: {
860: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
861: PetscInt i;
863: PetscFunctionBegin;
864: if (!ctx->subintset) {
865: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
866: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
867: }
868: PetscCall(PetscMalloc1(ctx->npart+1,subint));
869: for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: /*@
874: EPSKrylovSchurGetSubintervals - Returns the points that delimit the
875: subintervals used in spectrum slicing with several partitions.
877: Not Collective
879: Input Parameter:
880: . eps - the linear eigensolver context
882: Output Parameter:
883: . subint - array of real values specifying subintervals
885: Notes:
886: If the user passed values with `EPSKrylovSchurSetSubintervals()`, then the
887: same values are returned here. Otherwise, the values computed internally are
888: obtained.
890: This function is only available for spectrum slicing runs, that is, when
891: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
892: See more details in section [](#sec:slice).
894: The returned array has length `npart`+1 (see `EPSKrylovSchurGetPartitions()`)
895: and should be freed by the user.
897: Fortran Notes:
898: The calling sequence from Fortran is
899: .vb
900: EPSKrylovSchurGetSubintervals(eps,subint,ierr)
901: PetscReal subint(npart+1)
902: .ve
904: Level: advanced
906: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetSubintervals()`, `EPSKrylovSchurGetPartitions()`, `EPSSetInterval()`
907: @*/
908: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal *subint[]) PeNS
909: {
910: PetscFunctionBegin;
912: PetscAssertPointer(subint,2);
913: PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[])
918: {
919: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
920: PetscInt i,numsh;
921: EPS_SR sr = ctx->sr;
923: PetscFunctionBegin;
924: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
925: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
926: switch (eps->state) {
927: case EPS_STATE_INITIAL:
928: break;
929: case EPS_STATE_SETUP:
930: numsh = ctx->npart+1;
931: if (n) *n = numsh;
932: if (shifts) {
933: PetscCall(PetscMalloc1(numsh,shifts));
934: (*shifts)[0] = eps->inta;
935: if (ctx->npart==1) (*shifts)[1] = eps->intb;
936: else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
937: }
938: if (inertias) {
939: PetscCall(PetscMalloc1(numsh,inertias));
940: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
941: if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
942: else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
943: }
944: break;
945: case EPS_STATE_SOLVED:
946: case EPS_STATE_EIGENVECTORS:
947: numsh = ctx->nshifts;
948: if (n) *n = numsh;
949: if (shifts) {
950: PetscCall(PetscMalloc1(numsh,shifts));
951: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
952: }
953: if (inertias) {
954: PetscCall(PetscMalloc1(numsh,inertias));
955: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
956: }
957: break;
958: }
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }
962: /*@
963: EPSKrylovSchurGetInertias - Gets the values of the shifts and their
964: corresponding inertias in case of doing spectrum slicing for a
965: computational interval.
967: Not Collective
969: Input Parameter:
970: . eps - the linear eigensolver context
972: Output Parameters:
973: + n - number of shifts, including the endpoints of the interval
974: . shifts - the values of the shifts used internally in the solver
975: - inertias - the values of the inertia at each shift
977: Notes:
978: If called after `EPSSolve()`, all shifts used internally by the solver are
979: returned (including both endpoints and any intermediate ones). If called
980: before `EPSSolve()` and after `EPSSetUp()` then only the information of the
981: endpoints of subintervals is available.
983: This function is only available for spectrum slicing runs, that is, when
984: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
985: See more details in section [](#sec:slice).
987: The returned arrays should be freed by the user. Can pass `NULL` in any of
988: the two arrays if not required.
990: Fortran Notes:
991: The calling sequence from Fortran is
992: .vb
993: EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
994: PetscInt n
995: PetscReal shifts(*)
996: PetscInt inertias(*)
997: .ve
998: The arrays should be at least of length `n`. The value of `n` can be determined
999: by an initial call
1000: .vb
1001: EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr)
1002: .ve
1004: Level: advanced
1006: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetSubintervals()`
1007: @*/
1008: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[]) PeNS
1009: {
1010: PetscFunctionBegin;
1012: PetscAssertPointer(n,2);
1013: PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
1014: PetscFunctionReturn(PETSC_SUCCESS);
1015: }
1017: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1018: {
1019: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1020: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1022: PetscFunctionBegin;
1023: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1024: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1025: if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
1026: if (n) *n = sr->numEigs;
1027: if (v) PetscCall(BVCreateVec(sr->V,v));
1028: PetscFunctionReturn(PETSC_SUCCESS);
1029: }
1031: /*@
1032: EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
1033: doing spectrum slicing for a computational interval with multiple
1034: communicators.
1036: Collective on the subcommunicator (if `v` is given)
1038: Input Parameter:
1039: . eps - the linear eigensolver context
1041: Output Parameters:
1042: + k - index of the subinterval for the calling process
1043: . n - number of eigenvalues found in the `k`-th subinterval
1044: - v - a vector owned by processes in the subcommunicator with dimensions
1045: compatible for locally computed eigenvectors (or `NULL`)
1047: Notes:
1048: This function is only available for spectrum slicing runs, that is, when
1049: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1050: See more details in section [](#sec:slice).
1052: The returned `Vec` should be destroyed by the user.
1054: Level: advanced
1056: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommPairs()`
1057: @*/
1058: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1059: {
1060: PetscFunctionBegin;
1062: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1063: PetscFunctionReturn(PETSC_SUCCESS);
1064: }
1066: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1067: {
1068: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1069: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1071: PetscFunctionBegin;
1072: EPSCheckSolved(eps,1);
1073: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1074: PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1075: if (eig) *eig = sr->eigr[sr->perm[i]];
1076: if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1077: PetscFunctionReturn(PETSC_SUCCESS);
1078: }
1080: /*@
1081: EPSKrylovSchurGetSubcommPairs - Gets the `i`-th eigenpair stored
1082: internally in the subcommunicator to which the calling process belongs.
1084: Collective on the subcommunicator (if `v` is given)
1086: Input Parameters:
1087: + eps - the linear eigensolver context
1088: - i - index of the solution
1090: Output Parameters:
1091: + eig - the eigenvalue
1092: - v - the eigenvector
1094: Notes:
1095: This function is only available for spectrum slicing runs, that is, when
1096: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1097: And is relevant only when the number of partitions (`EPSKrylovSchurSetPartitions()`)
1098: is larger than one. See more details in section [](#sec:slice).
1100: It is allowed to pass `NULL` for `v` if the eigenvector is not required.
1101: Otherwise, the caller must provide a valid `Vec` object, i.e.,
1102: it must be created by the calling program with `EPSKrylovSchurGetSubcommInfo()`.
1104: The index `i` should be a value between 0 and `n`-1, where `n` is the number of
1105: vectors in the local subinterval, see `EPSKrylovSchurGetSubcommInfo()`.
1107: Level: advanced
1109: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`, `EPSKrylovSchurGetSubcommMats()`
1110: @*/
1111: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1112: {
1113: PetscFunctionBegin;
1116: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1117: PetscFunctionReturn(PETSC_SUCCESS);
1118: }
1120: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1121: {
1122: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1124: PetscFunctionBegin;
1125: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1126: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1127: PetscCall(EPSGetOperators(ctx->eps,A,B));
1128: PetscFunctionReturn(PETSC_SUCCESS);
1129: }
1131: /*@
1132: EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1133: internally in the subcommunicator to which the calling process belongs.
1135: Collective on the subcommunicator
1137: Input Parameter:
1138: . eps - the linear eigensolver context
1140: Output Parameters:
1141: + A - the matrix associated with the eigensystem
1142: - B - the second matrix in the case of generalized eigenproblems
1144: Notes:
1145: This function is only available for spectrum slicing runs, that is, when
1146: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1147: And is relevant only when the number of partitions (`EPSKrylovSchurSetPartitions()`)
1148: is larger than one. See more details in section [](#sec:slice).
1150: This is the analog of `EPSGetOperators()`, but returns the matrices distributed
1151: differently (in the subcommunicator rather than in the parent communicator).
1153: These matrices should not be modified by the user.
1155: Level: advanced
1157: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`
1158: @*/
1159: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1160: {
1161: PetscFunctionBegin;
1163: PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1164: PetscFunctionReturn(PETSC_SUCCESS);
1165: }
1167: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1168: {
1169: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1170: Mat A,B=NULL,Ag,Bg=NULL;
1171: PetscBool reuse=PETSC_TRUE;
1173: PetscFunctionBegin;
1174: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1175: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1176: PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1177: PetscCall(EPSGetOperators(ctx->eps,&A,&B));
1179: PetscCall(MatScale(A,a));
1180: if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1181: if (B) PetscCall(MatScale(B,b));
1182: if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1183: PetscCall(EPSSetOperators(ctx->eps,A,B));
1185: /* Update stored matrix state */
1186: subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1187: PetscCall(MatGetState(A,&subctx->Astate));
1188: if (B) PetscCall(MatGetState(B,&subctx->Bstate));
1190: /* Update matrices in the parent communicator if requested by user */
1191: if (globalup) {
1192: if (ctx->npart>1) {
1193: if (!ctx->isrow) {
1194: PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1195: reuse = PETSC_FALSE;
1196: }
1197: if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1198: if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1199: PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1200: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1201: if (B) {
1202: if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1203: PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1204: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1205: }
1206: }
1207: PetscCall(MatGetState(Ag,&ctx->Astate));
1208: if (Bg) PetscCall(MatGetState(Bg,&ctx->Bstate));
1209: }
1210: PetscCall(EPSSetOperators(eps,Ag,Bg));
1211: PetscFunctionReturn(PETSC_SUCCESS);
1212: }
1214: /*@
1215: EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1216: internally in the subcommunicator to which the calling process belongs.
1218: Collective
1220: Input Parameters:
1221: + eps - the linear eigensolver context
1222: . s - scalar that multiplies the existing $A$ matrix
1223: . a - scalar used in the _axpy_ operation on $A$
1224: . Au - matrix used in the _axpy_ operation on $A$
1225: . t - scalar that multiplies the existing $B$ matrix
1226: . b - scalar used in the _axpy_ operation on $B$
1227: . Bu - matrix used in the _axpy_ operation on $B$
1228: . str - structure flag, see `MatStructure`
1229: - globalup - flag indicating if global matrices must be updated
1231: Notes:
1232: This function is only available for spectrum slicing runs, that is, when
1233: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1234: And is relevant only when the number of partitions (`EPSKrylovSchurSetPartitions()`)
1235: is larger than one. See more details in section [](#sec:slice).
1237: This function modifies the eigenproblem matrices at the subcommunicator level,
1238: and optionally updates the global matrices in the parent communicator. The updates
1239: are expressed as $A \leftarrow s A + a A_u$ and $B \leftarrow t B + b B_u$.
1241: It is possible to update one of the matrices, or both.
1243: The matrices `Au` and `Bu` must be equal in all subcommunicators.
1245: The `str` flag is passed to the `MatAXPY()` operations to perform the updates.
1247: If `globalup` is `PETSC_TRUE`, communication is carried out to reconstruct the updated
1248: matrices in the parent communicator. The user must be warned that if global
1249: matrices are not in sync with subcommunicator matrices, the errors computed
1250: by `EPSComputeError()` will be wrong even if the computed solution is correct
1251: (the synchronization may be done only once at the end).
1253: Level: advanced
1255: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommMats()`
1256: @*/
1257: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1258: {
1259: PetscFunctionBegin;
1269: PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1270: PetscFunctionReturn(PETSC_SUCCESS);
1271: }
1273: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1274: {
1275: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1276: Mat A,B=NULL,Ar=NULL,Br=NULL;
1277: PetscMPIInt rank;
1278: PetscBool Asame,Bsame;
1279: STType sttype;
1280: PetscInt nmat;
1281: const char *prefix;
1282: MPI_Comm child;
1284: PetscFunctionBegin;
1285: PetscCall(EPSGetOperators(eps,&A,&B));
1286: if (ctx->npart==1) {
1287: if (!ctx->eps) PetscCall(EPSCreate(((PetscObject)eps)->comm,&ctx->eps));
1288: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1289: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1290: PetscCall(EPSSetOperators(ctx->eps,A,B));
1291: } else {
1292: if (!ctx->subc) {
1293: /* Create context for subcommunicators */
1294: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1295: PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1296: PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1297: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1299: /* Duplicate matrices */
1300: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1301: PetscCall(MatGetState(A,&ctx->Astate));
1302: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1303: if (B) {
1304: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1305: PetscCall(MatGetState(B,&ctx->Bstate));
1306: PetscCall(MatPropagateSymmetryOptions(B,Br));
1307: }
1308: } else {
1309: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1310: PetscCall(MatStateCompareUpdate(A,&ctx->Astate,&Asame));
1311: if (B) PetscCall(MatStateCompareUpdate(B,&ctx->Bstate,&Bsame));
1312: else Bsame = PETSC_TRUE;
1313: if (!Asame || !Bsame) {
1314: PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1315: if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
1316: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1317: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1318: if (B) {
1319: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1320: PetscCall(MatPropagateSymmetryOptions(B,Br));
1321: }
1322: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1323: PetscCall(MatDestroy(&Ar));
1324: PetscCall(MatDestroy(&Br));
1325: }
1326: }
1328: /* Create auxiliary EPS */
1329: if (!ctx->eps) {
1330: PetscCall(EPSCreate(child,&ctx->eps));
1331: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1332: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1333: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1334: PetscCall(MatDestroy(&Ar));
1335: PetscCall(MatDestroy(&Br));
1336: }
1337: /* Create subcommunicator grouping processes with same rank */
1338: if (!ctx->commset) {
1339: PetscCallMPI(MPI_Comm_rank(child,&rank));
1340: PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1341: ctx->commset = PETSC_TRUE;
1342: }
1343: }
1344: PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1345: PetscCall(STGetType(eps->st,&sttype));
1346: PetscCall(STSetType(ctx->eps->st,sttype));
1348: ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1349: ctx_local->npart = ctx->npart;
1350: ctx_local->global = PETSC_FALSE;
1351: ctx_local->eps = eps;
1352: ctx_local->subc = ctx->subc;
1353: ctx_local->commrank = ctx->commrank;
1354: *childeps = ctx->eps;
1355: PetscFunctionReturn(PETSC_SUCCESS);
1356: }
1358: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1359: {
1360: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1361: ST st;
1362: PetscBool isfilt;
1364: PetscFunctionBegin;
1365: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1366: PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1367: PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1368: PetscCall(EPSGetST(ctx->eps,&st));
1369: PetscCall(STGetOperator(st,NULL));
1370: PetscCall(STGetKSP(st,ksp));
1371: PetscFunctionReturn(PETSC_SUCCESS);
1372: }
1374: /*@
1375: EPSKrylovSchurGetKSP - Retrieve the linear solver object associated with the
1376: internal `EPS` object in case of doing spectrum slicing for a computational interval.
1378: Collective
1380: Input Parameter:
1381: . eps - the linear eigensolver context
1383: Output Parameter:
1384: . ksp - the internal `KSP` object
1386: Notes:
1387: When invoked to compute all eigenvalues in an interval with spectrum
1388: slicing, `EPSKRYLOVSCHUR` creates another `EPS` object internally that is
1389: used to compute eigenvalues by chunks near selected shifts. This function
1390: allows access to the `KSP` object associated to this internal `EPS` object.
1392: This function is only available for spectrum slicing runs, that is, when
1393: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1394: See more details in section [](#sec:slice).
1396: In case of having more than one partition, the returned `KSP` will be different
1397: in MPI processes belonging to different partitions. Hence, if required,
1398: `EPSKrylovSchurSetPartitions()` must be called BEFORE this function.
1400: Level: advanced
1402: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`
1403: @*/
1404: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1405: {
1406: PetscFunctionBegin;
1408: PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1409: PetscFunctionReturn(PETSC_SUCCESS);
1410: }
1412: static PetscErrorCode EPSKrylovSchurSetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType bse)
1413: {
1414: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1416: PetscFunctionBegin;
1417: switch (bse) {
1418: case EPS_KRYLOVSCHUR_BSE_SHAO:
1419: case EPS_KRYLOVSCHUR_BSE_GRUNING:
1420: case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
1421: if (ctx->bse != bse) {
1422: ctx->bse = bse;
1423: eps->state = EPS_STATE_INITIAL;
1424: }
1425: break;
1426: default:
1427: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid BSE type");
1428: }
1429: PetscFunctionReturn(PETSC_SUCCESS);
1430: }
1432: /*@
1433: EPSKrylovSchurSetBSEType - Sets the method to be used for BSE structured eigenproblems
1434: in the Krylov-Schur solver.
1436: Logically Collective
1438: Input Parameters:
1439: + eps - the linear eigensolver context
1440: - bse - the BSE method, see `EPSKrylovSchurBSEType` for possible values
1442: Options Database Key:
1443: . -eps_krylovschur_bse_type (shao|gruning|projectedbse) - sets the BSE type
1445: Notes:
1446: This function is relevant only for `EPS_BSE` problem types, see section
1447: on [](#sec:structured).
1449: A detailed description of the methods can be found in {cite:p}`Alv25`.
1451: Level: advanced
1453: .seealso: [](ch:eps), [](#sec:structured), `EPS_BSE`, `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1454: @*/
1455: PetscErrorCode EPSKrylovSchurSetBSEType(EPS eps,EPSKrylovSchurBSEType bse)
1456: {
1457: PetscFunctionBegin;
1460: PetscTryMethod(eps,"EPSKrylovSchurSetBSEType_C",(EPS,EPSKrylovSchurBSEType),(eps,bse));
1461: PetscFunctionReturn(PETSC_SUCCESS);
1462: }
1464: static PetscErrorCode EPSKrylovSchurGetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType *bse)
1465: {
1466: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1468: PetscFunctionBegin;
1469: *bse = ctx->bse;
1470: PetscFunctionReturn(PETSC_SUCCESS);
1471: }
1473: /*@
1474: EPSKrylovSchurGetBSEType - Gets the method used for BSE structured eigenproblems
1475: in the Krylov-Schur solver.
1477: Not Collective
1479: Input Parameter:
1480: . eps - the linear eigensolver context
1482: Output Parameter:
1483: . bse - the BSE method
1485: Level: advanced
1487: .seealso: [](ch:eps), [](#sec:structured), `EPS_BSE`, `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1488: @*/
1489: PetscErrorCode EPSKrylovSchurGetBSEType(EPS eps,EPSKrylovSchurBSEType *bse)
1490: {
1491: PetscFunctionBegin;
1493: PetscAssertPointer(bse,2);
1494: PetscUseMethod(eps,"EPSKrylovSchurGetBSEType_C",(EPS,EPSKrylovSchurBSEType*),(eps,bse));
1495: PetscFunctionReturn(PETSC_SUCCESS);
1496: }
1498: static PetscErrorCode EPSKrylovSchurSetLREPType_KrylovSchur(EPS eps,EPSKrylovSchurLREPType lrep)
1499: {
1500: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1502: PetscFunctionBegin;
1503: switch (lrep) {
1504: case EPS_KRYLOVSCHUR_LREP_TENG:
1505: case EPS_KRYLOVSCHUR_LREP_ZHONG:
1506: if (ctx->lrep != lrep) {
1507: ctx->lrep = lrep;
1508: eps->state = EPS_STATE_INITIAL;
1509: }
1510: break;
1511: default:
1512: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid LREP type");
1513: }
1514: PetscFunctionReturn(PETSC_SUCCESS);
1515: }
1517: /*@
1518: EPSKrylovSchurSetLREPType - Sets the method to be used for LREP structured eigenproblems
1519: in the Krylov-Schur solver.
1521: Logically Collective
1523: Input Parameters:
1524: + eps - the linear eigensolver context
1525: - lrep - the LREP method, see `EPSKrylovSchurLREPType` for possible values
1527: Options Database Key:
1528: . -eps_krylovschur_lrep_type (teng|zhong) - sets the LREP type
1530: Notes:
1531: This function is relevant only for `EPS_LREP` problem types, see section
1532: on [](#sec:structured).
1534: Level: advanced
1536: .seealso: [](ch:eps), [](#sec:structured), `EPS_LREP`, `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetLREPType()`, `EPSKrylovSchurLREPType`, `MatCreateLREP()`
1537: @*/
1538: PetscErrorCode EPSKrylovSchurSetLREPType(EPS eps,EPSKrylovSchurLREPType lrep)
1539: {
1540: PetscFunctionBegin;
1543: PetscTryMethod(eps,"EPSKrylovSchurSetLREPType_C",(EPS,EPSKrylovSchurLREPType),(eps,lrep));
1544: PetscFunctionReturn(PETSC_SUCCESS);
1545: }
1547: static PetscErrorCode EPSKrylovSchurGetLREPType_KrylovSchur(EPS eps,EPSKrylovSchurLREPType *lrep)
1548: {
1549: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1551: PetscFunctionBegin;
1552: *lrep = ctx->lrep;
1553: PetscFunctionReturn(PETSC_SUCCESS);
1554: }
1556: /*@
1557: EPSKrylovSchurGetLREPType - Gets the method used for LREP structured eigenproblems
1558: in the Krylov-Schur solver.
1560: Not Collective
1562: Input Parameter:
1563: . eps - the linear eigensolver context
1565: Output Parameter:
1566: . lrep - the LREP method
1568: Level: advanced
1570: .seealso: [](ch:eps), [](#sec:structured), `EPS_LREP`, `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetLREPType()`, `EPSKrylovSchurLREPType`, `MatCreateLREP()`
1571: @*/
1572: PetscErrorCode EPSKrylovSchurGetLREPType(EPS eps,EPSKrylovSchurLREPType *lrep)
1573: {
1574: PetscFunctionBegin;
1576: PetscAssertPointer(lrep,2);
1577: PetscUseMethod(eps,"EPSKrylovSchurGetLREPType_C",(EPS,EPSKrylovSchurLREPType*),(eps,lrep));
1578: PetscFunctionReturn(PETSC_SUCCESS);
1579: }
1581: static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems PetscOptionsObject)
1582: {
1583: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1584: PetscBool flg,lock,b,f1,f2,f3,isfilt;
1585: PetscReal keep;
1586: PetscInt i,j,k;
1587: KSP ksp;
1588: EPSKrylovSchurBSEType bse;
1589: EPSKrylovSchurLREPType lrep;
1591: PetscFunctionBegin;
1592: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");
1594: PetscCall(PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg));
1595: if (flg) PetscCall(EPSKrylovSchurSetRestart(eps,keep));
1597: PetscCall(PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg));
1598: if (flg) PetscCall(EPSKrylovSchurSetLocking(eps,lock));
1600: i = ctx->npart;
1601: PetscCall(PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg));
1602: if (flg) PetscCall(EPSKrylovSchurSetPartitions(eps,i));
1604: b = ctx->detect;
1605: PetscCall(PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg));
1606: if (flg) PetscCall(EPSKrylovSchurSetDetectZeros(eps,b));
1608: i = 1;
1609: j = k = PETSC_DECIDE;
1610: PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1611: PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1612: PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1613: if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));
1615: PetscCall(PetscOptionsEnum("-eps_krylovschur_bse_type","Method for BSE structured eigenproblems","EPSKrylovSchurSetBSEType",EPSKrylovSchurBSETypes,(PetscEnum)ctx->bse,(PetscEnum*)&bse,&flg));
1616: if (flg) PetscCall(EPSKrylovSchurSetBSEType(eps,bse));
1618: PetscCall(PetscOptionsEnum("-eps_krylovschur_lrep_type","Method for LREP structured eigenproblems","EPSKrylovSchurSetLREPType",EPSKrylovSchurLREPTypes,(PetscEnum)ctx->lrep,(PetscEnum*)&lrep,&flg));
1619: if (flg) PetscCall(EPSKrylovSchurSetLREPType(eps,lrep));
1621: PetscOptionsHeadEnd();
1623: /* set options of child KSP in spectrum slicing */
1624: if (eps->which==EPS_ALL) {
1625: if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1626: PetscCall(EPSSetDefaultST(eps));
1627: PetscCall(STSetFromOptions(eps->st)); /* need to advance this to check ST type */
1628: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1629: if (!isfilt) {
1630: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1631: PetscCall(KSPSetFromOptions(ksp));
1632: }
1633: }
1634: PetscFunctionReturn(PETSC_SUCCESS);
1635: }
1637: static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1638: {
1639: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1640: PetscBool isascii,isfilt;
1641: KSP ksp;
1642: PetscViewer sviewer;
1644: PetscFunctionBegin;
1645: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1646: if (isascii) {
1647: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1648: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1649: if (eps->problem_type==EPS_BSE) PetscCall(PetscViewerASCIIPrintf(viewer," BSE method: %s\n",EPSKrylovSchurBSETypes[ctx->bse]));
1650: else if (eps->problem_type==EPS_LREP) PetscCall(PetscViewerASCIIPrintf(viewer," LREP method: %s\n",EPSKrylovSchurLREPTypes[ctx->lrep]));
1651: if (eps->which==EPS_ALL) {
1652: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1653: if (isfilt) PetscCall(PetscViewerASCIIPrintf(viewer," using polynomial filtering%s\n",eps->nev?"":" to extract all eigenvalues in an interval"));
1654: else {
1655: PetscCall(PetscViewerASCIIPrintf(viewer," doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd));
1656: if (ctx->npart>1) {
1657: PetscCall(PetscViewerASCIIPrintf(viewer," multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart));
1658: if (ctx->detect) PetscCall(PetscViewerASCIIPrintf(viewer," detecting zeros when factorizing at subinterval boundaries\n"));
1659: }
1660: /* view child KSP */
1661: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1662: PetscCall(PetscViewerASCIIPushTab(viewer));
1663: if (ctx->npart>1 && ctx->subc) {
1664: PetscCall(PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer));
1665: if (!ctx->subc->color) PetscCall(KSPView(ksp,sviewer));
1666: PetscCall(PetscViewerFlush(sviewer));
1667: PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer));
1668: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1669: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1670: } else PetscCall(KSPView(ksp,viewer));
1671: PetscCall(PetscViewerASCIIPopTab(viewer));
1672: }
1673: }
1674: }
1675: PetscFunctionReturn(PETSC_SUCCESS);
1676: }
1678: static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1679: {
1680: PetscBool isfilt;
1682: PetscFunctionBegin;
1683: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1684: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1685: PetscCall(PetscFree(eps->data));
1686: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1687: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1688: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1689: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1690: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1691: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1692: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1693: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1694: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1695: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1696: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1697: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1698: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1699: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1700: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1701: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1702: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1703: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1704: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",NULL));
1705: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",NULL));
1706: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLREPType_C",NULL));
1707: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLREPType_C",NULL));
1708: PetscFunctionReturn(PETSC_SUCCESS);
1709: }
1711: static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1712: {
1713: PetscBool isfilt;
1715: PetscFunctionBegin;
1716: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1717: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1718: PetscFunctionReturn(PETSC_SUCCESS);
1719: }
1721: static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1722: {
1723: PetscFunctionBegin;
1724: if (eps->which==EPS_ALL) {
1725: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1726: }
1727: PetscFunctionReturn(PETSC_SUCCESS);
1728: }
1730: /*MC
1731: EPSKRYLOVSCHUR - EPSKRYLOVSCHUR = "krylovschur" - Krylov-Schur method.
1733: Notes:
1734: This is the default solver, and is recommended in most situations.
1736: The implemented algorithm is the single-vector Krylov-Schur method proposed
1737: by {cite:t}`Ste01b` for non-Hermitian eigenproblems. When the problem is
1738: Hermitian, the solver will behave as a thick-restart Lanczos method
1739: {cite:p}`Wu00` with full reorthogonalization.
1741: The solver includes support for many features\:
1742: - Harmonic extraction, see `EPSSetExtraction()`.
1743: - Inertia-based spectrum slicing to compute all eigenvalues in an interval,
1744: see {cite:p}`Cam12`.
1745: - Polynomial filter to compute eigenvalues in an interval via `STFILTER`.
1746: - Indefinite Lanczos to solve `EPS_GHIEP` problems.
1747: - Structured variants for problem types such as `EPS_BSE`.
1748: - A two-sided variant that also computes left eigenvectors.
1749: - Arbitrary selection of eigenvalues, see `EPSSetArbitrarySelection()`.
1751: Developer Note:
1752: In the future, we would like to have a block version of Krylov-Schur.
1754: Level: beginner
1756: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`, `EPSSetProblemType()`, `EPSSetExtraction()`, `EPSSetArbitrarySelection()`, `EPSSetTwoSided()`
1757: M*/
1758: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1759: {
1760: EPS_KRYLOVSCHUR *ctx;
1762: PetscFunctionBegin;
1763: PetscCall(PetscNew(&ctx));
1764: eps->data = (void*)ctx;
1765: ctx->lock = PETSC_TRUE;
1766: ctx->nev = 1;
1767: ctx->ncv = PETSC_DETERMINE;
1768: ctx->mpd = PETSC_DETERMINE;
1769: ctx->npart = 1;
1770: ctx->detect = PETSC_FALSE;
1771: ctx->global = PETSC_TRUE;
1773: eps->useds = PETSC_TRUE;
1775: /* solve and computevectors determined at setup */
1776: eps->ops->setup = EPSSetUp_KrylovSchur;
1777: eps->ops->setupsort = EPSSetUpSort_KrylovSchur;
1778: eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1779: eps->ops->destroy = EPSDestroy_KrylovSchur;
1780: eps->ops->reset = EPSReset_KrylovSchur;
1781: eps->ops->view = EPSView_KrylovSchur;
1782: eps->ops->backtransform = EPSBackTransform_Default;
1783: eps->ops->setdefaultst = EPSSetDefaultST_KrylovSchur;
1785: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1786: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1787: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1788: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1789: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1790: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1791: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1792: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1793: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1794: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1795: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1796: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1797: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1798: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1799: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1800: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1801: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1802: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1803: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",EPSKrylovSchurSetBSEType_KrylovSchur));
1804: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",EPSKrylovSchurGetBSEType_KrylovSchur));
1805: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLREPType_C",EPSKrylovSchurSetLREPType_KrylovSchur));
1806: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLREPType_C",EPSKrylovSchurGetLREPType_KrylovSchur));
1807: PetscFunctionReturn(PETSC_SUCCESS);
1808: }