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;
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 | EPS_FEATURE_THRESHOLD,PETSC_TRUE," with polynomial filter");
73: if (eps->tol==(PetscReal)PETSC_DETERMINE) eps->tol = SLEPC_DEFAULT_TOL*1e-2; /* use tighter tolerance */
74: PetscCall(STFilterSetInterval(eps->st,eps->inta,eps->intb));
75: if (!ctx->estimatedrange) {
76: PetscCall(STFilterGetRange(eps->st,&rleft,&rright));
77: estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
78: }
79: if (estimaterange) { /* user did not set a range */
80: PetscCall(STGetMatrix(eps->st,0,&A));
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: if (eps->ncv==PETSC_DETERMINE && eps->nev==0) eps->nev = 40; /* user did not provide nev estimation */
87: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
88: 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: 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 SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unknown matrix structure");
113: PetscFunctionReturn(PETSC_SUCCESS);
114: } else {
115: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
116: 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");
117: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);
118: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
119: }
120: PetscCheck(ctx->lock || eps->mpd>=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
122: EPSCheckDefiniteCondition(eps,eps->arbitrary," with arbitrary selection of eigenpairs");
124: PetscCheck(eps->extraction==EPS_RITZ || eps->extraction==EPS_HARMONIC,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
126: if (!ctx->keep) ctx->keep = 0.5;
128: PetscCall(EPSAllocateSolution(eps,1));
129: PetscCall(EPS_SetInnerProduct(eps));
130: if (eps->arbitrary) PetscCall(EPSSetWorkVecs(eps,2));
131: else if (eps->ishermitian && !eps->ispositive) PetscCall(EPSSetWorkVecs(eps,1));
133: /* dispatch solve method */
134: if (eps->ishermitian) {
135: if (eps->which==EPS_ALL) {
136: EPSCheckDefiniteCondition(eps,eps->which==EPS_ALL," with spectrum slicing");
137: variant = isfilt? EPS_KS_FILTER: EPS_KS_SLICE;
138: } else if (eps->isgeneralized && !eps->ispositive) {
139: variant = EPS_KS_INDEF;
140: } else {
141: switch (eps->extraction) {
142: case EPS_RITZ: variant = EPS_KS_SYMM; break;
143: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
144: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
145: }
146: }
147: } else if (eps->twosided) {
148: variant = EPS_KS_TWOSIDED;
149: } else {
150: switch (eps->extraction) {
151: case EPS_RITZ: variant = EPS_KS_DEFAULT; break;
152: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
153: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
154: }
155: }
156: switch (variant) {
157: case EPS_KS_DEFAULT:
158: eps->ops->solve = EPSSolve_KrylovSchur_Default;
159: eps->ops->computevectors = EPSComputeVectors_Schur;
160: PetscCall(DSSetType(eps->ds,DSNHEP));
161: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
162: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
163: break;
164: case EPS_KS_SYMM:
165: case EPS_KS_FILTER:
166: eps->ops->solve = EPSSolve_KrylovSchur_Default;
167: eps->ops->computevectors = EPSComputeVectors_Hermitian;
168: PetscCall(DSSetType(eps->ds,DSHEP));
169: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
170: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
171: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
172: break;
173: case EPS_KS_SLICE:
174: eps->ops->solve = EPSSolve_KrylovSchur_Slice;
175: eps->ops->computevectors = EPSComputeVectors_Slice;
176: break;
177: case EPS_KS_INDEF:
178: eps->ops->solve = EPSSolve_KrylovSchur_Indefinite;
179: eps->ops->computevectors = EPSComputeVectors_Indefinite;
180: PetscCall(DSSetType(eps->ds,DSGHIEP));
181: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
182: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
183: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
184: /* force reorthogonalization for pseudo-Lanczos */
185: PetscCall(BVGetOrthogonalization(eps->V,&otype,NULL,&eta,&obtype));
186: PetscCall(BVSetOrthogonalization(eps->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
187: break;
188: case EPS_KS_TWOSIDED:
189: eps->ops->solve = EPSSolve_KrylovSchur_TwoSided;
190: eps->ops->computevectors = EPSComputeVectors_Schur;
191: PetscCall(DSSetType(eps->ds,DSNHEPTS));
192: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
193: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
194: break;
195: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
196: }
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: static PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
201: {
202: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
203: SlepcSC sc;
204: PetscBool isfilt;
206: PetscFunctionBegin;
207: PetscCall(EPSSetUpSort_Default(eps));
208: if (eps->which==EPS_ALL) {
209: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
210: if (isfilt) {
211: PetscCall(DSGetSlepcSC(eps->ds,&sc));
212: sc->rg = NULL;
213: sc->comparison = SlepcCompareLargestReal;
214: sc->comparisonctx = NULL;
215: sc->map = NULL;
216: sc->mapobj = NULL;
217: } else {
218: if (!ctx->global && ctx->sr->numEigs>0) {
219: PetscCall(DSGetSlepcSC(eps->ds,&sc));
220: sc->rg = NULL;
221: sc->comparison = SlepcCompareLargestMagnitude;
222: sc->comparisonctx = NULL;
223: sc->map = NULL;
224: sc->mapobj = NULL;
225: }
226: }
227: }
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: PetscErrorCode EPSSolve_KrylovSchur_Default(EPS eps)
232: {
233: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
234: PetscInt i,j,*pj,k,l,nv,ld,nconv;
235: Mat U,Op,H,T;
236: PetscScalar *g;
237: PetscReal beta,gamma=1.0;
238: PetscBool breakdown,harmonic,hermitian;
240: PetscFunctionBegin;
241: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
242: harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
243: hermitian = (eps->ishermitian && !harmonic)?PETSC_TRUE:PETSC_FALSE;
244: if (harmonic) PetscCall(PetscMalloc1(ld,&g));
245: if (eps->arbitrary) pj = &j;
246: else pj = NULL;
248: /* Get the starting Arnoldi vector */
249: PetscCall(EPSGetStartVector(eps,0,NULL));
250: l = 0;
252: /* Restart loop */
253: while (eps->reason == EPS_CONVERGED_ITERATING) {
254: eps->its++;
256: /* Compute an nv-step Arnoldi factorization */
257: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
258: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
259: PetscCall(STGetOperator(eps->st,&Op));
260: if (hermitian) {
261: PetscCall(DSGetMat(eps->ds,DS_MAT_T,&T));
262: PetscCall(BVMatLanczos(eps->V,Op,T,eps->nconv+l,&nv,&beta,&breakdown));
263: PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&T));
264: } else {
265: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
266: PetscCall(BVMatArnoldi(eps->V,Op,H,eps->nconv+l,&nv,&beta,&breakdown));
267: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
268: }
269: PetscCall(STRestoreOperator(eps->st,&Op));
270: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
271: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
272: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
274: /* Compute translation of Krylov decomposition if harmonic extraction used */
275: if (PetscUnlikely(harmonic)) PetscCall(DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,g,&gamma));
277: /* Solve projected problem */
278: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
279: if (PetscUnlikely(eps->arbitrary)) {
280: PetscCall(EPSGetArbitraryValues(eps,eps->rr,eps->ri));
281: j=1;
282: }
283: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj));
284: PetscCall(DSUpdateExtraRow(eps->ds));
285: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
287: /* Check convergence */
288: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
289: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
290: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
291: nconv = k;
293: /* Update l */
294: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
295: else {
296: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
297: if (!hermitian) PetscCall(DSGetTruncateSize(eps->ds,k,nv,&l));
298: }
299: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
300: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
302: if (eps->reason == EPS_CONVERGED_ITERATING) {
303: if (PetscUnlikely(breakdown || k==nv)) {
304: /* Start a new Arnoldi factorization */
305: PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
306: if (k<eps->nev) {
307: PetscCall(EPSGetStartVector(eps,k,&breakdown));
308: if (breakdown) {
309: eps->reason = EPS_DIVERGED_BREAKDOWN;
310: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
311: }
312: }
313: } else {
314: /* Undo translation of Krylov decomposition */
315: if (PetscUnlikely(harmonic)) {
316: PetscCall(DSSetDimensions(eps->ds,nv,k,l));
317: PetscCall(DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma));
318: /* gamma u^ = u - U*g~ */
319: PetscCall(BVSetActiveColumns(eps->V,0,nv));
320: PetscCall(BVMultColumn(eps->V,-1.0,1.0,nv,g));
321: PetscCall(BVScaleColumn(eps->V,nv,1.0/gamma));
322: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
323: PetscCall(DSSetDimensions(eps->ds,nv,k,nv));
324: }
325: /* Prepare the Rayleigh quotient for restart */
326: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
327: }
328: }
329: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
330: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
331: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
332: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
334: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
335: PetscCall(BVCopyColumn(eps->V,nv,k+l)); /* copy restart vector from the last column */
336: if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
337: eps->ncv = eps->mpd+k;
338: PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
339: for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
340: PetscCall(DSReallocate(eps->ds,eps->ncv+1));
341: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
342: }
343: }
345: eps->nconv = k;
346: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
347: }
349: if (harmonic) PetscCall(PetscFree(g));
350: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
355: {
356: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
358: PetscFunctionBegin;
359: if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
360: else {
361: 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);
362: ctx->keep = keep;
363: }
364: PetscFunctionReturn(PETSC_SUCCESS);
365: }
367: /*@
368: EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
369: method, in particular the proportion of basis vectors that must be kept
370: after restart.
372: Logically Collective
374: Input Parameters:
375: + eps - the linear eigensolver context
376: - keep - the number of vectors to be kept at restart
378: Options Database Key:
379: . -eps_krylovschur_restart \<keep\> - sets the restart parameter
381: Notes:
382: Allowed values are in the range [0.1,0.9]. The default is 0.5, which means
383: that at restart the current subspace is compressed into another subspace
384: with a reduction of 50% in size.
386: Implementation details of Krylov-Schur in SLEPc can be found in {cite:p}`Her07b`.
388: Level: advanced
390: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetRestart()`
391: @*/
392: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
393: {
394: PetscFunctionBegin;
397: PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
402: {
403: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
405: PetscFunctionBegin;
406: *keep = ctx->keep;
407: PetscFunctionReturn(PETSC_SUCCESS);
408: }
410: /*@
411: EPSKrylovSchurGetRestart - Gets the restart parameter used in the
412: Krylov-Schur method.
414: Not Collective
416: Input Parameter:
417: . eps - the linear eigensolver context
419: Output Parameter:
420: . keep - the restart parameter
422: Level: advanced
424: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetRestart()`
425: @*/
426: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
427: {
428: PetscFunctionBegin;
430: PetscAssertPointer(keep,2);
431: PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
436: {
437: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
439: PetscFunctionBegin;
440: ctx->lock = lock;
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: /*@
445: EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
446: the Krylov-Schur method.
448: Logically Collective
450: Input Parameters:
451: + eps - the linear eigensolver context
452: - lock - true if the locking variant must be selected
454: Options Database Key:
455: . -eps_krylovschur_locking - sets the locking flag
457: Notes:
458: The default is to lock converged eigenpairs when the method restarts.
459: This behavior can be changed so that all directions are kept in the
460: working subspace even if already converged to working accuracy (the
461: non-locking variant).
463: Implementation details of Krylov-Schur in SLEPc can be found in {cite:p}`Her07b`.
465: Level: advanced
467: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetLocking()`
468: @*/
469: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
470: {
471: PetscFunctionBegin;
474: PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
479: {
480: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
482: PetscFunctionBegin;
483: *lock = ctx->lock;
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: /*@
488: EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
489: method.
491: Not Collective
493: Input Parameter:
494: . eps - the linear eigensolver context
496: Output Parameter:
497: . lock - the locking flag
499: Level: advanced
501: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetLocking()`
502: @*/
503: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
504: {
505: PetscFunctionBegin;
507: PetscAssertPointer(lock,2);
508: PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
513: {
514: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
515: PetscMPIInt size;
516: PetscInt newnpart;
518: PetscFunctionBegin;
519: if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
520: newnpart = 1;
521: } else {
522: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
523: PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
524: newnpart = npart;
525: }
526: if (ctx->npart!=newnpart) {
527: if (ctx->npart>1) {
528: PetscCall(PetscSubcommDestroy(&ctx->subc));
529: if (ctx->commset) {
530: PetscCallMPI(MPI_Comm_free(&ctx->commrank));
531: ctx->commset = PETSC_FALSE;
532: }
533: }
534: PetscCall(EPSDestroy(&ctx->eps));
535: ctx->npart = newnpart;
536: eps->state = EPS_STATE_INITIAL;
537: }
538: PetscFunctionReturn(PETSC_SUCCESS);
539: }
541: /*@
542: EPSKrylovSchurSetPartitions - Sets the number of partitions for the
543: case of doing spectrum slicing for a computational interval with the
544: communicator split in several sub-communicators.
546: Logically Collective
548: Input Parameters:
549: + eps - the linear eigensolver context
550: - npart - number of partitions
552: Options Database Key:
553: . -eps_krylovschur_partitions \<npart\> - sets the number of partitions
555: Notes:
556: This call makes sense only for spectrum slicing runs, that is, when
557: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
558: See more details in section [](#sec:slice).
560: By default, `npart`=1 so all processes in the communicator participate in
561: the processing of the whole interval. If `npart`>1 then the interval is
562: divided into `npart` subintervals, each of them being processed by a
563: subset of processes.
565: The interval is split proportionally unless the separation points are
566: specified with `EPSKrylovSchurSetSubintervals()`.
568: Level: advanced
570: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetSubintervals()`, `EPSSetInterval()`
571: @*/
572: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
573: {
574: PetscFunctionBegin;
577: PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
582: {
583: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
585: PetscFunctionBegin;
586: *npart = ctx->npart;
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: /*@
591: EPSKrylovSchurGetPartitions - Gets the number of partitions of the
592: communicator in case of spectrum slicing.
594: Not Collective
596: Input Parameter:
597: . eps - the linear eigensolver context
599: Output Parameter:
600: . npart - number of partitions
602: Level: advanced
604: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetPartitions()`
605: @*/
606: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
607: {
608: PetscFunctionBegin;
610: PetscAssertPointer(npart,2);
611: PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
616: {
617: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
619: PetscFunctionBegin;
620: ctx->detect = detect;
621: eps->state = EPS_STATE_INITIAL;
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: /*@
626: EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
627: zeros during the factorizations throughout the spectrum slicing computation.
629: Logically Collective
631: Input Parameters:
632: + eps - the linear eigensolver context
633: - detect - check for zeros
635: Options Database Key:
636: . -eps_krylovschur_detect_zeros - check for zeros
638: Notes:
639: This flag makes sense only for spectrum slicing runs, that is, when
640: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
641: See more details in section [](#sec:slice).
643: A zero in the factorization indicates that a shift coincides with an eigenvalue.
645: This flag is turned off by default, and may be necessary in some cases,
646: especially when several partitions are being used. This feature currently
647: requires an external package for factorizations with support for zero
648: detection, e.g. MUMPS.
650: Level: advanced
652: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetPartitions()`, `EPSSetInterval()`
653: @*/
654: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
655: {
656: PetscFunctionBegin;
659: PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
660: PetscFunctionReturn(PETSC_SUCCESS);
661: }
663: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
664: {
665: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
667: PetscFunctionBegin;
668: *detect = ctx->detect;
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: /*@
673: EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
674: in spectrum slicing.
676: Not Collective
678: Input Parameter:
679: . eps - the linear eigensolver context
681: Output Parameter:
682: . detect - whether zeros detection is enforced during factorizations
684: Level: advanced
686: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetDetectZeros()`
687: @*/
688: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
689: {
690: PetscFunctionBegin;
692: PetscAssertPointer(detect,2);
693: PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
694: PetscFunctionReturn(PETSC_SUCCESS);
695: }
697: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
698: {
699: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
701: PetscFunctionBegin;
702: if (nev != PETSC_CURRENT) {
703: PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
704: ctx->nev = nev;
705: }
706: if (ncv == PETSC_DETERMINE) {
707: ctx->ncv = PETSC_DETERMINE;
708: } else if (ncv != PETSC_CURRENT) {
709: PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
710: ctx->ncv = ncv;
711: }
712: if (mpd == PETSC_DETERMINE) {
713: ctx->mpd = PETSC_DETERMINE;
714: } else if (mpd != PETSC_CURRENT) {
715: PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
716: ctx->mpd = mpd;
717: }
718: eps->state = EPS_STATE_INITIAL;
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }
722: /*@
723: EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
724: step in case of doing spectrum slicing for a computational interval.
726: Logically Collective
728: Input Parameters:
729: + eps - the linear eigensolver context
730: . nev - number of eigenvalues to compute
731: . ncv - the maximum dimension of the subspace to be used by the subsolve
732: - mpd - the maximum dimension allowed for the projected problem
734: Options Database Keys:
735: + -eps_krylovschur_nev \<nev\> - sets the number of eigenvalues
736: . -eps_krylovschur_ncv \<ncv\> - sets the dimension of the subspace
737: - -eps_krylovschur_mpd \<mpd\> - sets the maximum projected dimension
739: Notes:
740: These parameters are relevant only for spectrum slicing runs, that is, when
741: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
742: See more details in section [](#sec:slice).
744: The meaning of the parameters is the same as in `EPSSetDimensions()`, but
745: the ones here apply to every subsolve done by the child `EPS` object.
747: Use `PETSC_DETERMINE` for `ncv` and `mpd` to assign a default value. For any
748: of the arguments, use `PETSC_CURRENT` to preserve the current value.
750: Level: advanced
752: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetDimensions()`, `EPSSetDimensions()`, `EPSSetInterval()`
753: @*/
754: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
755: {
756: PetscFunctionBegin;
761: PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
762: PetscFunctionReturn(PETSC_SUCCESS);
763: }
765: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
766: {
767: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
769: PetscFunctionBegin;
770: if (nev) *nev = ctx->nev;
771: if (ncv) *ncv = ctx->ncv;
772: if (mpd) *mpd = ctx->mpd;
773: PetscFunctionReturn(PETSC_SUCCESS);
774: }
776: /*@
777: EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
778: step in case of doing spectrum slicing for a computational interval.
780: Not Collective
782: Input Parameter:
783: . eps - the linear eigensolver context
785: Output Parameters:
786: + nev - number of eigenvalues to compute
787: . ncv - the maximum dimension of the subspace to be used by the subsolve
788: - mpd - the maximum dimension allowed for the projected problem
790: Level: advanced
792: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetDimensions()`
793: @*/
794: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
795: {
796: PetscFunctionBegin;
798: PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
802: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
803: {
804: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
805: PetscInt i;
807: PetscFunctionBegin;
808: 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()");
809: 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");
810: if (ctx->subintervals) PetscCall(PetscFree(ctx->subintervals));
811: PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
812: for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
813: ctx->subintset = PETSC_TRUE;
814: eps->state = EPS_STATE_INITIAL;
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }
818: /*@
819: EPSKrylovSchurSetSubintervals - Sets the points that delimit the
820: subintervals to be used in spectrum slicing with several partitions.
822: Logically Collective
824: Input Parameters:
825: + eps - the linear eigensolver context
826: - subint - array of real values specifying subintervals
828: Notes:
829: This function is relevant only for spectrum slicing runs, that is, when
830: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
831: See more details in section [](#sec:slice).
833: It must be called after `EPSKrylovSchurSetPartitions()`. For `npart`
834: partitions, the argument `subint` must contain `npart+1` real values sorted in
835: ascending order, `subint_0, subint_1, ..., subint_npart`, where the first
836: and last values must coincide with the interval endpoints set with
837: `EPSSetInterval()`.
839: The subintervals are then defined by two consecutive points `[subint_0,subint_1]`,
840: `[subint_1,subint_2]`, and so on.
842: Level: advanced
844: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubintervals()`, `EPSSetInterval()`
845: @*/
846: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal subint[])
847: {
848: PetscFunctionBegin;
850: PetscAssertPointer(subint,2);
851: PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal *subint[])
856: {
857: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
858: PetscInt i;
860: PetscFunctionBegin;
861: if (!ctx->subintset) {
862: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
863: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
864: }
865: PetscCall(PetscMalloc1(ctx->npart+1,subint));
866: for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
867: PetscFunctionReturn(PETSC_SUCCESS);
868: }
870: /*@C
871: EPSKrylovSchurGetSubintervals - Returns the points that delimit the
872: subintervals used in spectrum slicing with several partitions.
874: Not Collective
876: Input Parameter:
877: . eps - the linear eigensolver context
879: Output Parameter:
880: . subint - array of real values specifying subintervals
882: Notes:
883: If the user passed values with `EPSKrylovSchurSetSubintervals()`, then the
884: same values are returned here. Otherwise, the values computed internally are
885: obtained.
887: This function is only available for spectrum slicing runs, that is, when
888: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
889: See more details in section [](#sec:slice).
891: The returned array has length `npart`+1 (see `EPSKrylovSchurGetPartitions()`)
892: and should be freed by the user.
894: Fortran Notes:
895: The calling sequence from Fortran is
896: .vb
897: EPSKrylovSchurGetSubintervals(eps,subint,ierr)
898: PetscReal subint(npart+1)
899: .ve
901: Level: advanced
903: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetSubintervals()`, `EPSKrylovSchurGetPartitions()`, `EPSSetInterval()`
904: @*/
905: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal *subint[]) PeNS
906: {
907: PetscFunctionBegin;
909: PetscAssertPointer(subint,2);
910: PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[])
915: {
916: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
917: PetscInt i,numsh;
918: EPS_SR sr = ctx->sr;
920: PetscFunctionBegin;
921: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
922: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
923: switch (eps->state) {
924: case EPS_STATE_INITIAL:
925: break;
926: case EPS_STATE_SETUP:
927: numsh = ctx->npart+1;
928: if (n) *n = numsh;
929: if (shifts) {
930: PetscCall(PetscMalloc1(numsh,shifts));
931: (*shifts)[0] = eps->inta;
932: if (ctx->npart==1) (*shifts)[1] = eps->intb;
933: else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
934: }
935: if (inertias) {
936: PetscCall(PetscMalloc1(numsh,inertias));
937: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
938: if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
939: else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
940: }
941: break;
942: case EPS_STATE_SOLVED:
943: case EPS_STATE_EIGENVECTORS:
944: numsh = ctx->nshifts;
945: if (n) *n = numsh;
946: if (shifts) {
947: PetscCall(PetscMalloc1(numsh,shifts));
948: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
949: }
950: if (inertias) {
951: PetscCall(PetscMalloc1(numsh,inertias));
952: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
953: }
954: break;
955: }
956: PetscFunctionReturn(PETSC_SUCCESS);
957: }
959: /*@C
960: EPSKrylovSchurGetInertias - Gets the values of the shifts and their
961: corresponding inertias in case of doing spectrum slicing for a
962: computational interval.
964: Not Collective
966: Input Parameter:
967: . eps - the linear eigensolver context
969: Output Parameters:
970: + n - number of shifts, including the endpoints of the interval
971: . shifts - the values of the shifts used internally in the solver
972: - inertias - the values of the inertia at each shift
974: Notes:
975: If called after `EPSSolve()`, all shifts used internally by the solver are
976: returned (including both endpoints and any intermediate ones). If called
977: before `EPSSolve()` and after `EPSSetUp()` then only the information of the
978: endpoints of subintervals is available.
980: This function is only available for spectrum slicing runs, that is, when
981: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
982: See more details in section [](#sec:slice).
984: The returned arrays should be freed by the user. Can pass `NULL` in any of
985: the two arrays if not required.
987: Fortran Notes:
988: The calling sequence from Fortran is
989: .vb
990: EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
991: PetscInt n
992: PetscReal shifts(*)
993: PetscInt inertias(*)
994: .ve
995: The arrays should be at least of length `n`. The value of `n` can be determined
996: by an initial call
997: .vb
998: EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr)
999: .ve
1001: Level: advanced
1003: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetSubintervals()`
1004: @*/
1005: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[]) PeNS
1006: {
1007: PetscFunctionBegin;
1009: PetscAssertPointer(n,2);
1010: PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
1011: PetscFunctionReturn(PETSC_SUCCESS);
1012: }
1014: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1015: {
1016: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1017: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1019: PetscFunctionBegin;
1020: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1021: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1022: if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
1023: if (n) *n = sr->numEigs;
1024: if (v) PetscCall(BVCreateVec(sr->V,v));
1025: PetscFunctionReturn(PETSC_SUCCESS);
1026: }
1028: /*@
1029: EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
1030: doing spectrum slicing for a computational interval with multiple
1031: communicators.
1033: Collective on the subcommunicator (if `v` is given)
1035: Input Parameter:
1036: . eps - the linear eigensolver context
1038: Output Parameters:
1039: + k - index of the subinterval for the calling process
1040: . n - number of eigenvalues found in the `k`-th subinterval
1041: - v - a vector owned by processes in the subcommunicator with dimensions
1042: compatible for locally computed eigenvectors (or `NULL`)
1044: Notes:
1045: This function is only available for spectrum slicing runs, that is, when
1046: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1047: See more details in section [](#sec:slice).
1049: The returned `Vec` should be destroyed by the user.
1051: Level: advanced
1053: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommPairs()`
1054: @*/
1055: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1056: {
1057: PetscFunctionBegin;
1059: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1060: PetscFunctionReturn(PETSC_SUCCESS);
1061: }
1063: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1064: {
1065: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1066: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1068: PetscFunctionBegin;
1069: EPSCheckSolved(eps,1);
1070: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1071: PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1072: if (eig) *eig = sr->eigr[sr->perm[i]];
1073: if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: /*@
1078: EPSKrylovSchurGetSubcommPairs - Gets the `i`-th eigenpair stored
1079: internally in the subcommunicator to which the calling process belongs.
1081: Collective on the subcommunicator (if `v` is given)
1083: Input Parameters:
1084: + eps - the linear eigensolver context
1085: - i - index of the solution
1087: Output Parameters:
1088: + eig - the eigenvalue
1089: - v - the eigenvector
1091: Notes:
1092: This function is only available for spectrum slicing runs, that is, when
1093: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1094: And is relevant only when the number of partitions (`EPSKrylovSchurSetPartitions()`)
1095: is larger than one. See more details in section [](#sec:slice).
1097: It is allowed to pass `NULL` for `v` if the eigenvector is not required.
1098: Otherwise, the caller must provide a valid `Vec` object, i.e.,
1099: it must be created by the calling program with `EPSKrylovSchurGetSubcommInfo()`.
1101: The index `i` should be a value between 0 and `n`-1, where `n` is the number of
1102: vectors in the local subinterval, see `EPSKrylovSchurGetSubcommInfo()`.
1104: Level: advanced
1106: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`, `EPSKrylovSchurGetSubcommMats()`
1107: @*/
1108: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1109: {
1110: PetscFunctionBegin;
1113: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1114: PetscFunctionReturn(PETSC_SUCCESS);
1115: }
1117: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1118: {
1119: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1121: PetscFunctionBegin;
1122: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1123: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1124: PetscCall(EPSGetOperators(ctx->eps,A,B));
1125: PetscFunctionReturn(PETSC_SUCCESS);
1126: }
1128: /*@
1129: EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1130: internally in the subcommunicator to which the calling process belongs.
1132: Collective on the subcommunicator
1134: Input Parameter:
1135: . eps - the linear eigensolver context
1137: Output Parameters:
1138: + A - the matrix associated with the eigensystem
1139: - B - the second matrix in the case of generalized eigenproblems
1141: Notes:
1142: This function is only available for spectrum slicing runs, that is, when
1143: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1144: And is relevant only when the number of partitions (`EPSKrylovSchurSetPartitions()`)
1145: is larger than one. See more details in section [](#sec:slice).
1147: This is the analog of `EPSGetOperators()`, but returns the matrices distributed
1148: differently (in the subcommunicator rather than in the parent communicator).
1150: These matrices should not be modified by the user.
1152: Level: advanced
1154: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`
1155: @*/
1156: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1157: {
1158: PetscFunctionBegin;
1160: PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1161: PetscFunctionReturn(PETSC_SUCCESS);
1162: }
1164: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1165: {
1166: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1167: Mat A,B=NULL,Ag,Bg=NULL;
1168: PetscBool reuse=PETSC_TRUE;
1170: PetscFunctionBegin;
1171: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1172: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1173: PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1174: PetscCall(EPSGetOperators(ctx->eps,&A,&B));
1176: PetscCall(MatScale(A,a));
1177: if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1178: if (B) PetscCall(MatScale(B,b));
1179: if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1180: PetscCall(EPSSetOperators(ctx->eps,A,B));
1182: /* Update stored matrix state */
1183: subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1184: PetscCall(MatGetState(A,&subctx->Astate));
1185: if (B) PetscCall(MatGetState(B,&subctx->Bstate));
1187: /* Update matrices in the parent communicator if requested by user */
1188: if (globalup) {
1189: if (ctx->npart>1) {
1190: if (!ctx->isrow) {
1191: PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1192: reuse = PETSC_FALSE;
1193: }
1194: if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1195: if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1196: PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1197: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1198: if (B) {
1199: if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1200: PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1201: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1202: }
1203: }
1204: PetscCall(MatGetState(Ag,&ctx->Astate));
1205: if (Bg) PetscCall(MatGetState(Bg,&ctx->Bstate));
1206: }
1207: PetscCall(EPSSetOperators(eps,Ag,Bg));
1208: PetscFunctionReturn(PETSC_SUCCESS);
1209: }
1211: /*@
1212: EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1213: internally in the subcommunicator to which the calling process belongs.
1215: Collective
1217: Input Parameters:
1218: + eps - the linear eigensolver context
1219: . s - scalar that multiplies the existing $A$ matrix
1220: . a - scalar used in the _axpy_ operation on $A$
1221: . Au - matrix used in the _axpy_ operation on $A$
1222: . t - scalar that multiplies the existing $B$ matrix
1223: . b - scalar used in the _axpy_ operation on $B$
1224: . Bu - matrix used in the _axpy_ operation on $B$
1225: . str - structure flag, see `MatStructure`
1226: - globalup - flag indicating if global matrices must be updated
1228: Notes:
1229: This function is only available for spectrum slicing runs, that is, when
1230: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1231: And is relevant only when the number of partitions (`EPSKrylovSchurSetPartitions()`)
1232: is larger than one. See more details in section [](#sec:slice).
1234: This function modifies the eigenproblem matrices at the subcommunicator level,
1235: and optionally updates the global matrices in the parent communicator. The updates
1236: are expressed as $A \leftarrow s A + a A_u$ and $B \leftarrow t B + b B_u$.
1238: It is possible to update one of the matrices, or both.
1240: The matrices `Au` and `Bu` must be equal in all subcommunicators.
1242: The `str` flag is passed to the `MatAXPY()` operations to perform the updates.
1244: If `globalup` is `PETSC_TRUE`, communication is carried out to reconstruct the updated
1245: matrices in the parent communicator. The user must be warned that if global
1246: matrices are not in sync with subcommunicator matrices, the errors computed
1247: by `EPSComputeError()` will be wrong even if the computed solution is correct
1248: (the synchronization may be done only once at the end).
1250: Level: advanced
1252: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommMats()`
1253: @*/
1254: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1255: {
1256: PetscFunctionBegin;
1266: PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1271: {
1272: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1273: Mat A,B=NULL,Ar=NULL,Br=NULL;
1274: PetscMPIInt rank;
1275: PetscObjectState Astate,Bstate=0;
1276: PetscObjectId Aid,Bid=0;
1277: STType sttype;
1278: PetscInt nmat;
1279: const char *prefix;
1280: MPI_Comm child;
1282: PetscFunctionBegin;
1283: PetscCall(EPSGetOperators(eps,&A,&B));
1284: if (ctx->npart==1) {
1285: if (!ctx->eps) PetscCall(EPSCreate(((PetscObject)eps)->comm,&ctx->eps));
1286: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1287: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1288: PetscCall(EPSSetOperators(ctx->eps,A,B));
1289: } else {
1290: PetscCall(MatGetState(A,&Astate));
1291: PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1292: if (B) {
1293: PetscCall(MatGetState(B,&Bstate));
1294: PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1295: }
1296: if (!ctx->subc) {
1297: /* Create context for subcommunicators */
1298: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1299: PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1300: PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1301: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1303: /* Duplicate matrices */
1304: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1305: ctx->Astate = Astate;
1306: ctx->Aid = Aid;
1307: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1308: if (B) {
1309: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1310: ctx->Bstate = Bstate;
1311: ctx->Bid = Bid;
1312: PetscCall(MatPropagateSymmetryOptions(B,Br));
1313: }
1314: } else {
1315: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1316: if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1317: PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1318: if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
1319: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1320: ctx->Astate = Astate;
1321: ctx->Aid = Aid;
1322: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1323: if (B) {
1324: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1325: ctx->Bstate = Bstate;
1326: ctx->Bid = Bid;
1327: PetscCall(MatPropagateSymmetryOptions(B,Br));
1328: }
1329: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1330: PetscCall(MatDestroy(&Ar));
1331: PetscCall(MatDestroy(&Br));
1332: }
1333: }
1335: /* Create auxiliary EPS */
1336: if (!ctx->eps) {
1337: PetscCall(EPSCreate(child,&ctx->eps));
1338: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1339: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1340: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1341: PetscCall(MatDestroy(&Ar));
1342: PetscCall(MatDestroy(&Br));
1343: }
1344: /* Create subcommunicator grouping processes with same rank */
1345: if (!ctx->commset) {
1346: PetscCallMPI(MPI_Comm_rank(child,&rank));
1347: PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1348: ctx->commset = PETSC_TRUE;
1349: }
1350: }
1351: PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1352: PetscCall(STGetType(eps->st,&sttype));
1353: PetscCall(STSetType(ctx->eps->st,sttype));
1355: ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1356: ctx_local->npart = ctx->npart;
1357: ctx_local->global = PETSC_FALSE;
1358: ctx_local->eps = eps;
1359: ctx_local->subc = ctx->subc;
1360: ctx_local->commrank = ctx->commrank;
1361: *childeps = ctx->eps;
1362: PetscFunctionReturn(PETSC_SUCCESS);
1363: }
1365: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1366: {
1367: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1368: ST st;
1369: PetscBool isfilt;
1371: PetscFunctionBegin;
1372: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1373: PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1374: PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1375: PetscCall(EPSGetST(ctx->eps,&st));
1376: PetscCall(STGetOperator(st,NULL));
1377: PetscCall(STGetKSP(st,ksp));
1378: PetscFunctionReturn(PETSC_SUCCESS);
1379: }
1381: /*@
1382: EPSKrylovSchurGetKSP - Retrieve the linear solver object associated with the
1383: internal `EPS` object in case of doing spectrum slicing for a computational interval.
1385: Collective
1387: Input Parameter:
1388: . eps - the linear eigensolver context
1390: Output Parameter:
1391: . ksp - the internal `KSP` object
1393: Notes:
1394: When invoked to compute all eigenvalues in an interval with spectrum
1395: slicing, `EPSKRYLOVSCHUR` creates another `EPS` object internally that is
1396: used to compute eigenvalues by chunks near selected shifts. This function
1397: allows access to the `KSP` object associated to this internal `EPS` object.
1399: This function is only available for spectrum slicing runs, that is, when
1400: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1401: See more details in section [](#sec:slice).
1403: In case of having more than one partition, the returned `KSP` will be different
1404: in MPI processes belonging to different partitions. Hence, if required,
1405: `EPSKrylovSchurSetPartitions()` must be called BEFORE this function.
1407: Level: advanced
1409: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`
1410: @*/
1411: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1412: {
1413: PetscFunctionBegin;
1415: PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1416: PetscFunctionReturn(PETSC_SUCCESS);
1417: }
1419: static PetscErrorCode EPSKrylovSchurSetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType bse)
1420: {
1421: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1423: PetscFunctionBegin;
1424: switch (bse) {
1425: case EPS_KRYLOVSCHUR_BSE_SHAO:
1426: case EPS_KRYLOVSCHUR_BSE_GRUNING:
1427: case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
1428: if (ctx->bse != bse) {
1429: ctx->bse = bse;
1430: eps->state = EPS_STATE_INITIAL;
1431: }
1432: break;
1433: default:
1434: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid BSE type");
1435: }
1436: PetscFunctionReturn(PETSC_SUCCESS);
1437: }
1439: /*@
1440: EPSKrylovSchurSetBSEType - Sets the method to be used for BSE structured eigenproblems
1441: in the Krylov-Schur solver.
1443: Logically Collective
1445: Input Parameters:
1446: + eps - the linear eigensolver context, see `EPSKrylovSchurBSEType` for possible values
1447: - bse - the BSE method
1449: Options Database Key:
1450: . -eps_krylovschur_bse_type \<bse\> - sets the BSE type, either `shao`, `gruning`, or `projectedbse`
1452: Notes:
1453: This function is relevant only for `EPS_BSE` problem types, see section
1454: on [](#sec:structured).
1456: A detailed description of the methods can be found in {cite:p}`Alv25`.
1458: Level: advanced
1460: .seealso: [](ch:eps), [](#sec:structured), `EPS_BSE`, `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1461: @*/
1462: PetscErrorCode EPSKrylovSchurSetBSEType(EPS eps,EPSKrylovSchurBSEType bse)
1463: {
1464: PetscFunctionBegin;
1467: PetscTryMethod(eps,"EPSKrylovSchurSetBSEType_C",(EPS,EPSKrylovSchurBSEType),(eps,bse));
1468: PetscFunctionReturn(PETSC_SUCCESS);
1469: }
1471: static PetscErrorCode EPSKrylovSchurGetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType *bse)
1472: {
1473: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1475: PetscFunctionBegin;
1476: *bse = ctx->bse;
1477: PetscFunctionReturn(PETSC_SUCCESS);
1478: }
1480: /*@
1481: EPSKrylovSchurGetBSEType - Gets the method used for BSE structured eigenproblems
1482: in the Krylov-Schur solver.
1484: Not Collective
1486: Input Parameter:
1487: . eps - the linear eigensolver context
1489: Output Parameter:
1490: . bse - the BSE method
1492: Level: advanced
1494: .seealso: [](ch:eps), [](#sec:structured), `EPS_BSE`, `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1495: @*/
1496: PetscErrorCode EPSKrylovSchurGetBSEType(EPS eps,EPSKrylovSchurBSEType *bse)
1497: {
1498: PetscFunctionBegin;
1500: PetscAssertPointer(bse,2);
1501: PetscUseMethod(eps,"EPSKrylovSchurGetBSEType_C",(EPS,EPSKrylovSchurBSEType*),(eps,bse));
1502: PetscFunctionReturn(PETSC_SUCCESS);
1503: }
1505: static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems PetscOptionsObject)
1506: {
1507: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1508: PetscBool flg,lock,b,f1,f2,f3,isfilt;
1509: PetscReal keep;
1510: PetscInt i,j,k;
1511: KSP ksp;
1512: EPSKrylovSchurBSEType bse;
1514: PetscFunctionBegin;
1515: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");
1517: PetscCall(PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg));
1518: if (flg) PetscCall(EPSKrylovSchurSetRestart(eps,keep));
1520: PetscCall(PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg));
1521: if (flg) PetscCall(EPSKrylovSchurSetLocking(eps,lock));
1523: i = ctx->npart;
1524: PetscCall(PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg));
1525: if (flg) PetscCall(EPSKrylovSchurSetPartitions(eps,i));
1527: b = ctx->detect;
1528: PetscCall(PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg));
1529: if (flg) PetscCall(EPSKrylovSchurSetDetectZeros(eps,b));
1531: i = 1;
1532: j = k = PETSC_DECIDE;
1533: PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1534: PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1535: PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1536: if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));
1538: PetscCall(PetscOptionsEnum("-eps_krylovschur_bse_type","Method for BSE structured eigenproblems","EPSKrylovSchurSetBSEType",EPSKrylovSchurBSETypes,(PetscEnum)ctx->bse,(PetscEnum*)&bse,&flg));
1539: if (flg) PetscCall(EPSKrylovSchurSetBSEType(eps,bse));
1541: PetscOptionsHeadEnd();
1543: /* set options of child KSP in spectrum slicing */
1544: if (eps->which==EPS_ALL) {
1545: if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1546: PetscCall(EPSSetDefaultST(eps));
1547: PetscCall(STSetFromOptions(eps->st)); /* need to advance this to check ST type */
1548: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1549: if (!isfilt) {
1550: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1551: PetscCall(KSPSetFromOptions(ksp));
1552: }
1553: }
1554: PetscFunctionReturn(PETSC_SUCCESS);
1555: }
1557: static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1558: {
1559: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1560: PetscBool isascii,isfilt;
1561: KSP ksp;
1562: PetscViewer sviewer;
1564: PetscFunctionBegin;
1565: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1566: if (isascii) {
1567: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1568: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1569: if (eps->problem_type==EPS_BSE) PetscCall(PetscViewerASCIIPrintf(viewer," BSE method: %s\n",EPSKrylovSchurBSETypes[ctx->bse]));
1570: if (eps->which==EPS_ALL) {
1571: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1572: if (isfilt) PetscCall(PetscViewerASCIIPrintf(viewer," using filtering to extract all eigenvalues in an interval\n"));
1573: else {
1574: PetscCall(PetscViewerASCIIPrintf(viewer," doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd));
1575: if (ctx->npart>1) {
1576: PetscCall(PetscViewerASCIIPrintf(viewer," multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart));
1577: if (ctx->detect) PetscCall(PetscViewerASCIIPrintf(viewer," detecting zeros when factorizing at subinterval boundaries\n"));
1578: }
1579: /* view child KSP */
1580: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1581: PetscCall(PetscViewerASCIIPushTab(viewer));
1582: if (ctx->npart>1 && ctx->subc) {
1583: PetscCall(PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer));
1584: if (!ctx->subc->color) PetscCall(KSPView(ksp,sviewer));
1585: PetscCall(PetscViewerFlush(sviewer));
1586: PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer));
1587: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1588: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1589: } else PetscCall(KSPView(ksp,viewer));
1590: PetscCall(PetscViewerASCIIPopTab(viewer));
1591: }
1592: }
1593: }
1594: PetscFunctionReturn(PETSC_SUCCESS);
1595: }
1597: static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1598: {
1599: PetscBool isfilt;
1601: PetscFunctionBegin;
1602: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1603: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1604: PetscCall(PetscFree(eps->data));
1605: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1606: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1607: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1608: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1609: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1610: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1611: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1612: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1613: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1614: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1615: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1616: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1617: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1618: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1619: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1620: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1621: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1622: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1623: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",NULL));
1624: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",NULL));
1625: PetscFunctionReturn(PETSC_SUCCESS);
1626: }
1628: static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1629: {
1630: PetscBool isfilt;
1632: PetscFunctionBegin;
1633: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1634: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1635: PetscFunctionReturn(PETSC_SUCCESS);
1636: }
1638: static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1639: {
1640: PetscFunctionBegin;
1641: if (eps->which==EPS_ALL) {
1642: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1643: }
1644: PetscFunctionReturn(PETSC_SUCCESS);
1645: }
1647: /*MC
1648: EPSKRYLOVSCHUR - EPSKRYLOVSCHUR = "krylovschur" - Krylov-Schur method.
1650: Notes:
1651: This is the default solver, and is recommended in most situations.
1653: The implemented algorithm is the single-vector Krylov-Schur method proposed
1654: by {cite:t}`Ste01b` for non-Hermitian eigenproblems. When the problem is
1655: Hermitian, the solver will behave as a thick-restart Lanczos method
1656: {cite:p}`Wu00` with full reorthogonalization.
1658: The solver includes support for many features\:
1659: - Harmonic extraction, see `EPSSetExtraction()`.
1660: - Inertia-based spectrum slicing to compute all eigenvalues in an interval,
1661: see {cite:p}`Cam12`.
1662: - Polynomial filter to compute eigenvalues in an interval via `STFILTER`.
1663: - Indefinite Lanczos to solve `EPSGHIEP` problems.
1664: - Structured variants for problem types such as `EPSBSE`.
1665: - A two-sided variant that also computes left eigenvectors.
1666: - Arbitrary selection of eigenvalues, see `EPSSetArbitrarySelection()`.
1668: Developer Note:
1669: In the future, we would like to have a block version of Krylov-Schur.
1671: Level: beginner
1673: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`, `EPSSetProblemType()`, `EPSSetExtraction()`, `EPSSetArbitrarySelection()`, `EPSSetTwoSided()`
1674: M*/
1675: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1676: {
1677: EPS_KRYLOVSCHUR *ctx;
1679: PetscFunctionBegin;
1680: PetscCall(PetscNew(&ctx));
1681: eps->data = (void*)ctx;
1682: ctx->lock = PETSC_TRUE;
1683: ctx->nev = 1;
1684: ctx->ncv = PETSC_DETERMINE;
1685: ctx->mpd = PETSC_DETERMINE;
1686: ctx->npart = 1;
1687: ctx->detect = PETSC_FALSE;
1688: ctx->global = PETSC_TRUE;
1690: eps->useds = PETSC_TRUE;
1692: /* solve and computevectors determined at setup */
1693: eps->ops->setup = EPSSetUp_KrylovSchur;
1694: eps->ops->setupsort = EPSSetUpSort_KrylovSchur;
1695: eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1696: eps->ops->destroy = EPSDestroy_KrylovSchur;
1697: eps->ops->reset = EPSReset_KrylovSchur;
1698: eps->ops->view = EPSView_KrylovSchur;
1699: eps->ops->backtransform = EPSBackTransform_Default;
1700: eps->ops->setdefaultst = EPSSetDefaultST_KrylovSchur;
1702: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1703: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1704: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1705: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1706: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1707: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1708: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1709: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1710: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1711: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1712: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1713: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1714: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1715: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1716: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1717: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1718: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1719: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1720: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",EPSKrylovSchurSetBSEType_KrylovSchur));
1721: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",EPSKrylovSchurGetBSEType_KrylovSchur));
1722: PetscFunctionReturn(PETSC_SUCCESS);
1723: }