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,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: if (estimaterange) { /* user did not set a range */
79: PetscCall(STGetMatrix(eps->st,0,&A));
80: PetscCall(MatEstimateSpectralRange_EPS(A,&rleft,&rright));
81: PetscCall(PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright));
82: PetscCall(STFilterSetRange(eps->st,rleft,rright));
83: ctx->estimatedrange = PETSC_TRUE;
84: }
85: PetscCall(EPSSetStoppingTest(eps,EPS_STOP_THRESHOLD));
86: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
87: 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");
88: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode EPSSetUp_KrylovSchur(EPS eps)
93: {
94: PetscReal eta;
95: PetscBool isfilt=PETSC_FALSE;
96: BVOrthogType otype;
97: BVOrthogBlockType obtype;
98: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
99: enum { EPS_KS_DEFAULT,EPS_KS_SYMM,EPS_KS_SLICE,EPS_KS_FILTER,EPS_KS_INDEF,EPS_KS_TWOSIDED } variant;
101: PetscFunctionBegin;
102: if (eps->which==EPS_ALL) { /* default values in case of spectrum slicing or polynomial filter */
103: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
104: if (isfilt) PetscCall(EPSSetUp_KrylovSchur_Filter(eps));
105: else PetscCall(EPSSetUp_KrylovSchur_Slice(eps));
106: } else if (eps->isstructured) {
107: if (eps->problem_type==EPS_BSE) PetscCall(EPSSetUp_KrylovSchur_BSE(eps));
108: else if (eps->problem_type==EPS_HAMILT) {
109: PetscCheck(!PetscDefined(USE_COMPLEX),PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The Hamiltonian Krylov-Schur eigensolver is not yet implemented for complex scalars");
110: PetscCall(EPSSetUp_KrylovSchur_Hamilt(eps));
111: } else if (eps->problem_type==EPS_LREP) {
112: PetscCheck(!PetscDefined(USE_COMPLEX),PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The LREP Krylov-Schur eigensolver does not support complex scalars, use BSE instead");
113: PetscCall(EPSSetUp_KrylovSchur_LREP(eps));
114: } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unknown matrix structure");
115: PetscFunctionReturn(PETSC_SUCCESS);
116: } else {
117: PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
118: 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");
119: if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);
120: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
121: }
122: PetscCheck(ctx->lock || eps->mpd>=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
124: EPSCheckDefiniteCondition(eps,eps->arbitrary," with arbitrary selection of eigenpairs");
126: PetscCheck(eps->extraction==EPS_RITZ || eps->extraction==EPS_HARMONIC,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
128: if (!ctx->keep) ctx->keep = 0.5;
130: PetscCall(EPSAllocateSolution(eps,1));
131: PetscCall(EPS_SetInnerProduct(eps));
132: if (eps->arbitrary) PetscCall(EPSSetWorkVecs(eps,2));
133: else if (eps->ishermitian && !eps->ispositive) PetscCall(EPSSetWorkVecs(eps,1));
135: /* dispatch solve method */
136: if (eps->ishermitian) {
137: if (eps->which==EPS_ALL) {
138: EPSCheckDefiniteCondition(eps,eps->which==EPS_ALL," with spectrum slicing");
139: variant = isfilt? EPS_KS_FILTER: EPS_KS_SLICE;
140: } else if (eps->isgeneralized && !eps->ispositive) {
141: variant = EPS_KS_INDEF;
142: } else {
143: switch (eps->extraction) {
144: case EPS_RITZ: variant = EPS_KS_SYMM; break;
145: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
146: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
147: }
148: }
149: } else if (eps->twosided) {
150: variant = EPS_KS_TWOSIDED;
151: } else {
152: switch (eps->extraction) {
153: case EPS_RITZ: variant = EPS_KS_DEFAULT; break;
154: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
155: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
156: }
157: }
158: switch (variant) {
159: case EPS_KS_DEFAULT:
160: eps->ops->solve = EPSSolve_KrylovSchur_Default;
161: eps->ops->computevectors = EPSComputeVectors_Schur;
162: PetscCall(DSSetType(eps->ds,DSNHEP));
163: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
164: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
165: break;
166: case EPS_KS_SYMM:
167: case EPS_KS_FILTER:
168: eps->ops->solve = EPSSolve_KrylovSchur_Default;
169: eps->ops->computevectors = EPSComputeVectors_Hermitian;
170: PetscCall(DSSetType(eps->ds,DSHEP));
171: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
172: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
173: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
174: break;
175: case EPS_KS_SLICE:
176: eps->ops->solve = EPSSolve_KrylovSchur_Slice;
177: eps->ops->computevectors = EPSComputeVectors_Slice;
178: break;
179: case EPS_KS_INDEF:
180: eps->ops->solve = EPSSolve_KrylovSchur_Indefinite;
181: eps->ops->computevectors = EPSComputeVectors_Indefinite;
182: PetscCall(DSSetType(eps->ds,DSGHIEP));
183: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
184: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
185: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
186: /* force reorthogonalization for pseudo-Lanczos */
187: PetscCall(BVGetOrthogonalization(eps->V,&otype,NULL,&eta,&obtype));
188: PetscCall(BVSetOrthogonalization(eps->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
189: break;
190: case EPS_KS_TWOSIDED:
191: eps->ops->solve = EPSSolve_KrylovSchur_TwoSided;
192: eps->ops->computevectors = EPSComputeVectors_Schur;
193: PetscCall(DSSetType(eps->ds,DSNHEPTS));
194: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
195: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
196: break;
197: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
198: }
199: PetscFunctionReturn(PETSC_SUCCESS);
200: }
202: static PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
203: {
204: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
205: SlepcSC sc;
206: PetscBool isfilt;
208: PetscFunctionBegin;
209: PetscCall(EPSSetUpSort_Default(eps));
210: if (eps->which==EPS_ALL) {
211: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
212: if (isfilt) {
213: PetscCall(DSGetSlepcSC(eps->ds,&sc));
214: sc->rg = NULL;
215: sc->comparison = SlepcCompareLargestReal;
216: sc->comparisonctx = NULL;
217: sc->map = NULL;
218: sc->mapobj = NULL;
219: } else {
220: if (!ctx->global && ctx->sr->numEigs>0) {
221: PetscCall(DSGetSlepcSC(eps->ds,&sc));
222: sc->rg = NULL;
223: sc->comparison = SlepcCompareLargestMagnitude;
224: sc->comparisonctx = NULL;
225: sc->map = NULL;
226: sc->mapobj = NULL;
227: }
228: }
229: }
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: PetscErrorCode EPSSolve_KrylovSchur_Default(EPS eps)
234: {
235: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
236: PetscInt i,j,*pj,k,l,nv,ld,nconv;
237: Mat U,Op,H,T;
238: PetscScalar *g;
239: PetscReal beta,gamma=1.0;
240: PetscBool breakdown,harmonic,hermitian;
242: PetscFunctionBegin;
243: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
244: harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
245: hermitian = (eps->ishermitian && !harmonic)?PETSC_TRUE:PETSC_FALSE;
246: if (harmonic) PetscCall(PetscMalloc1(ld,&g));
247: if (eps->arbitrary) pj = &j;
248: else pj = NULL;
250: /* Get the starting Arnoldi vector */
251: PetscCall(EPSGetStartVector(eps,0,NULL));
252: l = 0;
254: /* Restart loop */
255: while (eps->reason == EPS_CONVERGED_ITERATING) {
256: eps->its++;
258: /* Compute an nv-step Arnoldi factorization */
259: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
260: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
261: PetscCall(STGetOperator(eps->st,&Op));
262: if (hermitian) {
263: PetscCall(DSGetMat(eps->ds,DS_MAT_T,&T));
264: PetscCall(BVMatLanczos(eps->V,Op,T,eps->nconv+l,&nv,&beta,&breakdown));
265: PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&T));
266: } else {
267: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
268: PetscCall(BVMatArnoldi(eps->V,Op,H,eps->nconv+l,&nv,&beta,&breakdown));
269: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
270: }
271: PetscCall(STRestoreOperator(eps->st,&Op));
272: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
273: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
274: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
276: /* Compute translation of Krylov decomposition if harmonic extraction used */
277: if (PetscUnlikely(harmonic)) PetscCall(DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,g,&gamma));
279: /* Solve projected problem */
280: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
281: if (PetscUnlikely(eps->arbitrary)) {
282: PetscCall(EPSGetArbitraryValues(eps,eps->rr,eps->ri));
283: j=1;
284: }
285: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj));
286: PetscCall(DSUpdateExtraRow(eps->ds));
287: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
289: /* Check convergence */
290: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
291: EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,eps->errest,k,nv);
292: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
293: nconv = k;
295: /* Update l */
296: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
297: else {
298: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
299: if (!hermitian) PetscCall(DSGetTruncateSize(eps->ds,k,nv,&l));
300: }
301: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
302: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
304: if (eps->reason == EPS_CONVERGED_ITERATING) {
305: if (PetscUnlikely(breakdown || k==nv)) {
306: /* Start a new Arnoldi factorization */
307: PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
308: if (k<eps->nev) {
309: PetscCall(EPSGetStartVector(eps,k,&breakdown));
310: if (breakdown) {
311: eps->reason = EPS_DIVERGED_BREAKDOWN;
312: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
313: }
314: }
315: } else {
316: /* Undo translation of Krylov decomposition */
317: if (PetscUnlikely(harmonic)) {
318: PetscCall(DSSetDimensions(eps->ds,nv,k,l));
319: PetscCall(DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma));
320: /* gamma u^ = u - U*g~ */
321: PetscCall(BVSetActiveColumns(eps->V,0,nv));
322: PetscCall(BVMultColumn(eps->V,-1.0,1.0,nv,g));
323: PetscCall(BVScaleColumn(eps->V,nv,1.0/gamma));
324: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
325: PetscCall(DSSetDimensions(eps->ds,nv,k,nv));
326: }
327: /* Prepare the Rayleigh quotient for restart */
328: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
329: }
330: }
331: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
332: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
333: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
334: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
336: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
337: PetscCall(BVCopyColumn(eps->V,nv,k+l)); /* copy restart vector from the last column */
338: if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) { /* reallocate */
339: eps->ncv = eps->mpd+k;
340: PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
341: for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
342: PetscCall(DSReallocate(eps->ds,eps->ncv+1));
343: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
344: }
345: }
347: eps->nconv = k;
348: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
349: }
351: if (harmonic) PetscCall(PetscFree(g));
352: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
357: {
358: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
360: PetscFunctionBegin;
361: if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
362: else {
363: 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);
364: ctx->keep = keep;
365: }
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: /*@
370: EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
371: method, in particular the proportion of basis vectors that must be kept
372: after restart.
374: Logically Collective
376: Input Parameters:
377: + eps - the linear eigensolver context
378: - keep - the number of vectors to be kept at restart
380: Options Database Key:
381: . -eps_krylovschur_restart keep - sets the restart parameter
383: Notes:
384: Allowed values are in the range [0.1,0.9]. The default is 0.5, which means
385: that at restart the current subspace is compressed into another subspace
386: with a reduction of 50% in size.
388: Implementation details of Krylov-Schur in SLEPc can be found in {cite:p}`Her07b`.
390: Level: advanced
392: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetRestart()`
393: @*/
394: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
395: {
396: PetscFunctionBegin;
399: PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
400: PetscFunctionReturn(PETSC_SUCCESS);
401: }
403: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
404: {
405: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
407: PetscFunctionBegin;
408: *keep = ctx->keep;
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /*@
413: EPSKrylovSchurGetRestart - Gets the restart parameter used in the
414: Krylov-Schur method.
416: Not Collective
418: Input Parameter:
419: . eps - the linear eigensolver context
421: Output Parameter:
422: . keep - the restart parameter
424: Level: advanced
426: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetRestart()`
427: @*/
428: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
429: {
430: PetscFunctionBegin;
432: PetscAssertPointer(keep,2);
433: PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
438: {
439: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
441: PetscFunctionBegin;
442: ctx->lock = lock;
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /*@
447: EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
448: the Krylov-Schur method.
450: Logically Collective
452: Input Parameters:
453: + eps - the linear eigensolver context
454: - lock - true if the locking variant must be selected
456: Options Database Key:
457: . -eps_krylovschur_locking (true|false) - sets the locking flag
459: Notes:
460: The default is to lock converged eigenpairs when the method restarts.
461: This behavior can be changed so that all directions are kept in the
462: working subspace even if already converged to working accuracy (the
463: non-locking variant).
465: Implementation details of Krylov-Schur in SLEPc can be found in {cite:p}`Her07b`.
467: Level: advanced
469: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetLocking()`
470: @*/
471: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
472: {
473: PetscFunctionBegin;
476: PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
481: {
482: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
484: PetscFunctionBegin;
485: *lock = ctx->lock;
486: PetscFunctionReturn(PETSC_SUCCESS);
487: }
489: /*@
490: EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
491: method.
493: Not Collective
495: Input Parameter:
496: . eps - the linear eigensolver context
498: Output Parameter:
499: . lock - the locking flag
501: Level: advanced
503: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetLocking()`
504: @*/
505: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
506: {
507: PetscFunctionBegin;
509: PetscAssertPointer(lock,2);
510: PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
515: {
516: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
517: PetscMPIInt size;
518: PetscInt newnpart;
520: PetscFunctionBegin;
521: if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
522: newnpart = 1;
523: } else {
524: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
525: PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
526: newnpart = npart;
527: }
528: if (ctx->npart!=newnpart) {
529: if (ctx->npart>1) {
530: PetscCall(PetscSubcommDestroy(&ctx->subc));
531: if (ctx->commset) {
532: PetscCallMPI(MPI_Comm_free(&ctx->commrank));
533: ctx->commset = PETSC_FALSE;
534: }
535: }
536: PetscCall(EPSDestroy(&ctx->eps));
537: ctx->npart = newnpart;
538: eps->state = EPS_STATE_INITIAL;
539: }
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: /*@
544: EPSKrylovSchurSetPartitions - Sets the number of partitions for the
545: case of doing spectrum slicing for a computational interval with the
546: communicator split in several sub-communicators.
548: Logically Collective
550: Input Parameters:
551: + eps - the linear eigensolver context
552: - npart - number of partitions
554: Options Database Key:
555: . -eps_krylovschur_partitions npart - sets the number of partitions
557: Notes:
558: This call makes sense only for spectrum slicing runs, that is, when
559: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
560: See more details in section [](#sec:slice).
562: By default, `npart`=1 so all processes in the communicator participate in
563: the processing of the whole interval. If `npart`>1 then the interval is
564: divided into `npart` subintervals, each of them being processed by a
565: subset of processes.
567: The interval is split proportionally unless the separation points are
568: specified with `EPSKrylovSchurSetSubintervals()`.
570: Level: advanced
572: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetSubintervals()`, `EPSSetInterval()`
573: @*/
574: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
575: {
576: PetscFunctionBegin;
579: PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
580: PetscFunctionReturn(PETSC_SUCCESS);
581: }
583: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
584: {
585: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
587: PetscFunctionBegin;
588: *npart = ctx->npart;
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /*@
593: EPSKrylovSchurGetPartitions - Gets the number of partitions of the
594: communicator in case of spectrum slicing.
596: Not Collective
598: Input Parameter:
599: . eps - the linear eigensolver context
601: Output Parameter:
602: . npart - number of partitions
604: Level: advanced
606: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetPartitions()`
607: @*/
608: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
609: {
610: PetscFunctionBegin;
612: PetscAssertPointer(npart,2);
613: PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
618: {
619: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
621: PetscFunctionBegin;
622: ctx->detect = detect;
623: eps->state = EPS_STATE_INITIAL;
624: PetscFunctionReturn(PETSC_SUCCESS);
625: }
627: /*@
628: EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
629: zeros during the factorizations throughout the spectrum slicing computation.
631: Logically Collective
633: Input Parameters:
634: + eps - the linear eigensolver context
635: - detect - check for zeros
637: Options Database Key:
638: . -eps_krylovschur_detect_zeros (true|false) - check for zeros
640: Notes:
641: This flag makes sense only for spectrum slicing runs, that is, when
642: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
643: See more details in section [](#sec:slice).
645: A zero in the factorization indicates that a shift coincides with an eigenvalue.
647: This flag is turned off by default, and may be necessary in some cases,
648: especially when several partitions are being used. This feature currently
649: requires an external package for factorizations with support for zero
650: detection, e.g. MUMPS.
652: Level: advanced
654: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetPartitions()`, `EPSSetInterval()`
655: @*/
656: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
657: {
658: PetscFunctionBegin;
661: PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
666: {
667: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
669: PetscFunctionBegin;
670: *detect = ctx->detect;
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: /*@
675: EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
676: in spectrum slicing.
678: Not Collective
680: Input Parameter:
681: . eps - the linear eigensolver context
683: Output Parameter:
684: . detect - whether zeros detection is enforced during factorizations
686: Level: advanced
688: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetDetectZeros()`
689: @*/
690: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
691: {
692: PetscFunctionBegin;
694: PetscAssertPointer(detect,2);
695: PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
696: PetscFunctionReturn(PETSC_SUCCESS);
697: }
699: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
700: {
701: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
703: PetscFunctionBegin;
704: if (nev != PETSC_CURRENT) {
705: PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
706: ctx->nev = nev;
707: }
708: if (ncv == PETSC_DETERMINE) {
709: ctx->ncv = PETSC_DETERMINE;
710: } else if (ncv != PETSC_CURRENT) {
711: PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
712: ctx->ncv = ncv;
713: }
714: if (mpd == PETSC_DETERMINE) {
715: ctx->mpd = PETSC_DETERMINE;
716: } else if (mpd != PETSC_CURRENT) {
717: PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
718: ctx->mpd = mpd;
719: }
720: eps->state = EPS_STATE_INITIAL;
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: /*@
725: EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
726: step in case of doing spectrum slicing for a computational interval.
728: Logically Collective
730: Input Parameters:
731: + eps - the linear eigensolver context
732: . nev - number of eigenvalues to compute
733: . ncv - the maximum dimension of the subspace to be used by the subsolve
734: - mpd - the maximum dimension allowed for the projected problem
736: Options Database Keys:
737: + -eps_krylovschur_nev nev - sets the number of eigenvalues
738: . -eps_krylovschur_ncv ncv - sets the dimension of the subspace
739: - -eps_krylovschur_mpd mpd - sets the maximum projected dimension
741: Notes:
742: These parameters are relevant only for spectrum slicing runs, that is, when
743: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
744: See more details in section [](#sec:slice).
746: The meaning of the parameters is the same as in `EPSSetDimensions()`, but
747: the ones here apply to every subsolve done by the child `EPS` object.
749: Use `PETSC_DETERMINE` for `ncv` and `mpd` to assign a default value. For any
750: of the arguments, use `PETSC_CURRENT` to preserve the current value.
752: Level: advanced
754: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetDimensions()`, `EPSSetDimensions()`, `EPSSetInterval()`
755: @*/
756: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
757: {
758: PetscFunctionBegin;
763: PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
764: PetscFunctionReturn(PETSC_SUCCESS);
765: }
767: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
768: {
769: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
771: PetscFunctionBegin;
772: if (nev) *nev = ctx->nev;
773: if (ncv) *ncv = ctx->ncv;
774: if (mpd) *mpd = ctx->mpd;
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: /*@
779: EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
780: step in case of doing spectrum slicing for a computational interval.
782: Not Collective
784: Input Parameter:
785: . eps - the linear eigensolver context
787: Output Parameters:
788: + nev - number of eigenvalues to compute
789: . ncv - the maximum dimension of the subspace to be used by the subsolve
790: - mpd - the maximum dimension allowed for the projected problem
792: Level: advanced
794: .seealso: [](ch:eps), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetDimensions()`
795: @*/
796: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
797: {
798: PetscFunctionBegin;
800: PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
805: {
806: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
807: PetscInt i;
809: PetscFunctionBegin;
810: 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()");
811: 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");
812: PetscCall(PetscFree(ctx->subintervals));
813: PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
814: for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
815: ctx->subintset = PETSC_TRUE;
816: eps->state = EPS_STATE_INITIAL;
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
820: /*@
821: EPSKrylovSchurSetSubintervals - Sets the points that delimit the
822: subintervals to be used in spectrum slicing with several partitions.
824: Logically Collective
826: Input Parameters:
827: + eps - the linear eigensolver context
828: - subint - array of real values specifying subintervals
830: Notes:
831: This function is relevant only for spectrum slicing runs, that is, when
832: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
833: See more details in section [](#sec:slice).
835: It must be called after `EPSKrylovSchurSetPartitions()`. For `npart`
836: partitions, the argument `subint` must contain `npart+1` real values sorted in
837: ascending order, `subint_0, subint_1, ..., subint_npart`, where the first
838: and last values must coincide with the interval endpoints set with
839: `EPSSetInterval()`.
841: The subintervals are then defined by two consecutive points `[subint_0,subint_1]`,
842: `[subint_1,subint_2]`, and so on.
844: Level: advanced
846: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubintervals()`, `EPSSetInterval()`
847: @*/
848: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal subint[])
849: {
850: PetscFunctionBegin;
852: PetscAssertPointer(subint,2);
853: PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal *subint[])
858: {
859: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
860: PetscInt i;
862: PetscFunctionBegin;
863: if (!ctx->subintset) {
864: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
865: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
866: }
867: PetscCall(PetscMalloc1(ctx->npart+1,subint));
868: for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
869: PetscFunctionReturn(PETSC_SUCCESS);
870: }
872: /*@C
873: EPSKrylovSchurGetSubintervals - Returns the points that delimit the
874: subintervals used in spectrum slicing with several partitions.
876: Not Collective
878: Input Parameter:
879: . eps - the linear eigensolver context
881: Output Parameter:
882: . subint - array of real values specifying subintervals
884: Notes:
885: If the user passed values with `EPSKrylovSchurSetSubintervals()`, then the
886: same values are returned here. Otherwise, the values computed internally are
887: obtained.
889: This function is only available for spectrum slicing runs, that is, when
890: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
891: See more details in section [](#sec:slice).
893: The returned array has length `npart`+1 (see `EPSKrylovSchurGetPartitions()`)
894: and should be freed by the user.
896: Fortran Notes:
897: The calling sequence from Fortran is
898: .vb
899: EPSKrylovSchurGetSubintervals(eps,subint,ierr)
900: PetscReal subint(npart+1)
901: .ve
903: Level: advanced
905: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetSubintervals()`, `EPSKrylovSchurGetPartitions()`, `EPSSetInterval()`
906: @*/
907: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal *subint[]) PeNS
908: {
909: PetscFunctionBegin;
911: PetscAssertPointer(subint,2);
912: PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
913: PetscFunctionReturn(PETSC_SUCCESS);
914: }
916: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[])
917: {
918: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
919: PetscInt i,numsh;
920: EPS_SR sr = ctx->sr;
922: PetscFunctionBegin;
923: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
924: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
925: switch (eps->state) {
926: case EPS_STATE_INITIAL:
927: break;
928: case EPS_STATE_SETUP:
929: numsh = ctx->npart+1;
930: if (n) *n = numsh;
931: if (shifts) {
932: PetscCall(PetscMalloc1(numsh,shifts));
933: (*shifts)[0] = eps->inta;
934: if (ctx->npart==1) (*shifts)[1] = eps->intb;
935: else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
936: }
937: if (inertias) {
938: PetscCall(PetscMalloc1(numsh,inertias));
939: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
940: if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
941: else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
942: }
943: break;
944: case EPS_STATE_SOLVED:
945: case EPS_STATE_EIGENVECTORS:
946: numsh = ctx->nshifts;
947: if (n) *n = numsh;
948: if (shifts) {
949: PetscCall(PetscMalloc1(numsh,shifts));
950: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
951: }
952: if (inertias) {
953: PetscCall(PetscMalloc1(numsh,inertias));
954: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
955: }
956: break;
957: }
958: PetscFunctionReturn(PETSC_SUCCESS);
959: }
961: /*@C
962: EPSKrylovSchurGetInertias - Gets the values of the shifts and their
963: corresponding inertias in case of doing spectrum slicing for a
964: computational interval.
966: Not Collective
968: Input Parameter:
969: . eps - the linear eigensolver context
971: Output Parameters:
972: + n - number of shifts, including the endpoints of the interval
973: . shifts - the values of the shifts used internally in the solver
974: - inertias - the values of the inertia at each shift
976: Notes:
977: If called after `EPSSolve()`, all shifts used internally by the solver are
978: returned (including both endpoints and any intermediate ones). If called
979: before `EPSSolve()` and after `EPSSetUp()` then only the information of the
980: endpoints of subintervals is available.
982: This function is only available for spectrum slicing runs, that is, when
983: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
984: See more details in section [](#sec:slice).
986: The returned arrays should be freed by the user. Can pass `NULL` in any of
987: the two arrays if not required.
989: Fortran Notes:
990: The calling sequence from Fortran is
991: .vb
992: EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
993: PetscInt n
994: PetscReal shifts(*)
995: PetscInt inertias(*)
996: .ve
997: The arrays should be at least of length `n`. The value of `n` can be determined
998: by an initial call
999: .vb
1000: EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr)
1001: .ve
1003: Level: advanced
1005: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetSubintervals()`
1006: @*/
1007: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[]) PeNS
1008: {
1009: PetscFunctionBegin;
1011: PetscAssertPointer(n,2);
1012: PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1016: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1017: {
1018: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1019: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1021: PetscFunctionBegin;
1022: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1023: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1024: if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
1025: if (n) *n = sr->numEigs;
1026: if (v) PetscCall(BVCreateVec(sr->V,v));
1027: PetscFunctionReturn(PETSC_SUCCESS);
1028: }
1030: /*@
1031: EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
1032: doing spectrum slicing for a computational interval with multiple
1033: communicators.
1035: Collective on the subcommunicator (if `v` is given)
1037: Input Parameter:
1038: . eps - the linear eigensolver context
1040: Output Parameters:
1041: + k - index of the subinterval for the calling process
1042: . n - number of eigenvalues found in the `k`-th subinterval
1043: - v - a vector owned by processes in the subcommunicator with dimensions
1044: compatible for locally computed eigenvectors (or `NULL`)
1046: Notes:
1047: This function is only available for spectrum slicing runs, that is, when
1048: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1049: See more details in section [](#sec:slice).
1051: The returned `Vec` should be destroyed by the user.
1053: Level: advanced
1055: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommPairs()`
1056: @*/
1057: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1058: {
1059: PetscFunctionBegin;
1061: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1062: PetscFunctionReturn(PETSC_SUCCESS);
1063: }
1065: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1066: {
1067: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1068: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1070: PetscFunctionBegin;
1071: EPSCheckSolved(eps,1);
1072: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1073: PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1074: if (eig) *eig = sr->eigr[sr->perm[i]];
1075: if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1076: PetscFunctionReturn(PETSC_SUCCESS);
1077: }
1079: /*@
1080: EPSKrylovSchurGetSubcommPairs - Gets the `i`-th eigenpair stored
1081: internally in the subcommunicator to which the calling process belongs.
1083: Collective on the subcommunicator (if `v` is given)
1085: Input Parameters:
1086: + eps - the linear eigensolver context
1087: - i - index of the solution
1089: Output Parameters:
1090: + eig - the eigenvalue
1091: - v - the eigenvector
1093: Notes:
1094: This function is only available for spectrum slicing runs, that is, when
1095: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1096: And is relevant only when the number of partitions (`EPSKrylovSchurSetPartitions()`)
1097: is larger than one. See more details in section [](#sec:slice).
1099: It is allowed to pass `NULL` for `v` if the eigenvector is not required.
1100: Otherwise, the caller must provide a valid `Vec` object, i.e.,
1101: it must be created by the calling program with `EPSKrylovSchurGetSubcommInfo()`.
1103: The index `i` should be a value between 0 and `n`-1, where `n` is the number of
1104: vectors in the local subinterval, see `EPSKrylovSchurGetSubcommInfo()`.
1106: Level: advanced
1108: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`, `EPSKrylovSchurGetSubcommMats()`
1109: @*/
1110: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1111: {
1112: PetscFunctionBegin;
1115: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1116: PetscFunctionReturn(PETSC_SUCCESS);
1117: }
1119: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1120: {
1121: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1123: PetscFunctionBegin;
1124: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1125: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1126: PetscCall(EPSGetOperators(ctx->eps,A,B));
1127: PetscFunctionReturn(PETSC_SUCCESS);
1128: }
1130: /*@
1131: EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1132: internally in the subcommunicator to which the calling process belongs.
1134: Collective on the subcommunicator
1136: Input Parameter:
1137: . eps - the linear eigensolver context
1139: Output Parameters:
1140: + A - the matrix associated with the eigensystem
1141: - B - the second matrix in the case of generalized eigenproblems
1143: Notes:
1144: This function is only available for spectrum slicing runs, that is, when
1145: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1146: And is relevant only when the number of partitions (`EPSKrylovSchurSetPartitions()`)
1147: is larger than one. See more details in section [](#sec:slice).
1149: This is the analog of `EPSGetOperators()`, but returns the matrices distributed
1150: differently (in the subcommunicator rather than in the parent communicator).
1152: These matrices should not be modified by the user.
1154: Level: advanced
1156: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`
1157: @*/
1158: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1159: {
1160: PetscFunctionBegin;
1162: PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1163: PetscFunctionReturn(PETSC_SUCCESS);
1164: }
1166: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1167: {
1168: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1169: Mat A,B=NULL,Ag,Bg=NULL;
1170: PetscBool reuse=PETSC_TRUE;
1172: PetscFunctionBegin;
1173: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1174: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1175: PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1176: PetscCall(EPSGetOperators(ctx->eps,&A,&B));
1178: PetscCall(MatScale(A,a));
1179: if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1180: if (B) PetscCall(MatScale(B,b));
1181: if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1182: PetscCall(EPSSetOperators(ctx->eps,A,B));
1184: /* Update stored matrix state */
1185: subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1186: PetscCall(MatGetState(A,&subctx->Astate));
1187: if (B) PetscCall(MatGetState(B,&subctx->Bstate));
1189: /* Update matrices in the parent communicator if requested by user */
1190: if (globalup) {
1191: if (ctx->npart>1) {
1192: if (!ctx->isrow) {
1193: PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1194: reuse = PETSC_FALSE;
1195: }
1196: if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1197: if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1198: PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1199: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1200: if (B) {
1201: if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1202: PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1203: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1204: }
1205: }
1206: PetscCall(MatGetState(Ag,&ctx->Astate));
1207: if (Bg) PetscCall(MatGetState(Bg,&ctx->Bstate));
1208: }
1209: PetscCall(EPSSetOperators(eps,Ag,Bg));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: /*@
1214: EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1215: internally in the subcommunicator to which the calling process belongs.
1217: Collective
1219: Input Parameters:
1220: + eps - the linear eigensolver context
1221: . s - scalar that multiplies the existing $A$ matrix
1222: . a - scalar used in the _axpy_ operation on $A$
1223: . Au - matrix used in the _axpy_ operation on $A$
1224: . t - scalar that multiplies the existing $B$ matrix
1225: . b - scalar used in the _axpy_ operation on $B$
1226: . Bu - matrix used in the _axpy_ operation on $B$
1227: . str - structure flag, see `MatStructure`
1228: - globalup - flag indicating if global matrices must be updated
1230: Notes:
1231: This function is only available for spectrum slicing runs, that is, when
1232: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1233: And is relevant only when the number of partitions (`EPSKrylovSchurSetPartitions()`)
1234: is larger than one. See more details in section [](#sec:slice).
1236: This function modifies the eigenproblem matrices at the subcommunicator level,
1237: and optionally updates the global matrices in the parent communicator. The updates
1238: are expressed as $A \leftarrow s A + a A_u$ and $B \leftarrow t B + b B_u$.
1240: It is possible to update one of the matrices, or both.
1242: The matrices `Au` and `Bu` must be equal in all subcommunicators.
1244: The `str` flag is passed to the `MatAXPY()` operations to perform the updates.
1246: If `globalup` is `PETSC_TRUE`, communication is carried out to reconstruct the updated
1247: matrices in the parent communicator. The user must be warned that if global
1248: matrices are not in sync with subcommunicator matrices, the errors computed
1249: by `EPSComputeError()` will be wrong even if the computed solution is correct
1250: (the synchronization may be done only once at the end).
1252: Level: advanced
1254: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommMats()`
1255: @*/
1256: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1257: {
1258: PetscFunctionBegin;
1268: PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1269: PetscFunctionReturn(PETSC_SUCCESS);
1270: }
1272: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1273: {
1274: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1275: Mat A,B=NULL,Ar=NULL,Br=NULL;
1276: PetscMPIInt rank;
1277: PetscObjectState Astate,Bstate=0;
1278: PetscObjectId Aid,Bid=0;
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: PetscCall(MatGetState(A,&Astate));
1293: PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1294: if (B) {
1295: PetscCall(MatGetState(B,&Bstate));
1296: PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1297: }
1298: if (!ctx->subc) {
1299: /* Create context for subcommunicators */
1300: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1301: PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1302: PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1303: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1305: /* Duplicate matrices */
1306: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1307: ctx->Astate = Astate;
1308: ctx->Aid = Aid;
1309: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1310: if (B) {
1311: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1312: ctx->Bstate = Bstate;
1313: ctx->Bid = Bid;
1314: PetscCall(MatPropagateSymmetryOptions(B,Br));
1315: }
1316: } else {
1317: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1318: if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1319: PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1320: if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
1321: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1322: ctx->Astate = Astate;
1323: ctx->Aid = Aid;
1324: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1325: if (B) {
1326: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1327: ctx->Bstate = Bstate;
1328: ctx->Bid = Bid;
1329: PetscCall(MatPropagateSymmetryOptions(B,Br));
1330: }
1331: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1332: PetscCall(MatDestroy(&Ar));
1333: PetscCall(MatDestroy(&Br));
1334: }
1335: }
1337: /* Create auxiliary EPS */
1338: if (!ctx->eps) {
1339: PetscCall(EPSCreate(child,&ctx->eps));
1340: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1341: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1342: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1343: PetscCall(MatDestroy(&Ar));
1344: PetscCall(MatDestroy(&Br));
1345: }
1346: /* Create subcommunicator grouping processes with same rank */
1347: if (!ctx->commset) {
1348: PetscCallMPI(MPI_Comm_rank(child,&rank));
1349: PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1350: ctx->commset = PETSC_TRUE;
1351: }
1352: }
1353: PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1354: PetscCall(STGetType(eps->st,&sttype));
1355: PetscCall(STSetType(ctx->eps->st,sttype));
1357: ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1358: ctx_local->npart = ctx->npart;
1359: ctx_local->global = PETSC_FALSE;
1360: ctx_local->eps = eps;
1361: ctx_local->subc = ctx->subc;
1362: ctx_local->commrank = ctx->commrank;
1363: *childeps = ctx->eps;
1364: PetscFunctionReturn(PETSC_SUCCESS);
1365: }
1367: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1368: {
1369: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1370: ST st;
1371: PetscBool isfilt;
1373: PetscFunctionBegin;
1374: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1375: PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1376: PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1377: PetscCall(EPSGetST(ctx->eps,&st));
1378: PetscCall(STGetOperator(st,NULL));
1379: PetscCall(STGetKSP(st,ksp));
1380: PetscFunctionReturn(PETSC_SUCCESS);
1381: }
1383: /*@
1384: EPSKrylovSchurGetKSP - Retrieve the linear solver object associated with the
1385: internal `EPS` object in case of doing spectrum slicing for a computational interval.
1387: Collective
1389: Input Parameter:
1390: . eps - the linear eigensolver context
1392: Output Parameter:
1393: . ksp - the internal `KSP` object
1395: Notes:
1396: When invoked to compute all eigenvalues in an interval with spectrum
1397: slicing, `EPSKRYLOVSCHUR` creates another `EPS` object internally that is
1398: used to compute eigenvalues by chunks near selected shifts. This function
1399: allows access to the `KSP` object associated to this internal `EPS` object.
1401: This function is only available for spectrum slicing runs, that is, when
1402: an interval has been given with `EPSSetInterval()` and `STSINVERT` is set.
1403: See more details in section [](#sec:slice).
1405: In case of having more than one partition, the returned `KSP` will be different
1406: in MPI processes belonging to different partitions. Hence, if required,
1407: `EPSKrylovSchurSetPartitions()` must be called BEFORE this function.
1409: Level: advanced
1411: .seealso: [](ch:eps), [](#sec:slice), `EPSKRYLOVSCHUR`, `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`
1412: @*/
1413: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1414: {
1415: PetscFunctionBegin;
1417: PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1418: PetscFunctionReturn(PETSC_SUCCESS);
1419: }
1421: static PetscErrorCode EPSKrylovSchurSetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType bse)
1422: {
1423: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1425: PetscFunctionBegin;
1426: switch (bse) {
1427: case EPS_KRYLOVSCHUR_BSE_SHAO:
1428: case EPS_KRYLOVSCHUR_BSE_GRUNING:
1429: case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
1430: if (ctx->bse != bse) {
1431: ctx->bse = bse;
1432: eps->state = EPS_STATE_INITIAL;
1433: }
1434: break;
1435: default:
1436: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid BSE type");
1437: }
1438: PetscFunctionReturn(PETSC_SUCCESS);
1439: }
1441: /*@
1442: EPSKrylovSchurSetBSEType - Sets the method to be used for BSE structured eigenproblems
1443: in the Krylov-Schur solver.
1445: Logically Collective
1447: Input Parameters:
1448: + eps - the linear eigensolver context
1449: - bse - the BSE method, see `EPSKrylovSchurBSEType` for possible values
1451: Options Database Key:
1452: . -eps_krylovschur_bse_type (shao|gruning|projectedbse) - sets the BSE type
1454: Notes:
1455: This function is relevant only for `EPS_BSE` problem types, see section
1456: on [](#sec:structured).
1458: A detailed description of the methods can be found in {cite:p}`Alv25`.
1460: Level: advanced
1462: .seealso: [](ch:eps), [](#sec:structured), `EPS_BSE`, `EPSKRYLOVSCHUR`, `EPSKrylovSchurGetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1463: @*/
1464: PetscErrorCode EPSKrylovSchurSetBSEType(EPS eps,EPSKrylovSchurBSEType bse)
1465: {
1466: PetscFunctionBegin;
1469: PetscTryMethod(eps,"EPSKrylovSchurSetBSEType_C",(EPS,EPSKrylovSchurBSEType),(eps,bse));
1470: PetscFunctionReturn(PETSC_SUCCESS);
1471: }
1473: static PetscErrorCode EPSKrylovSchurGetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType *bse)
1474: {
1475: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1477: PetscFunctionBegin;
1478: *bse = ctx->bse;
1479: PetscFunctionReturn(PETSC_SUCCESS);
1480: }
1482: /*@
1483: EPSKrylovSchurGetBSEType - Gets the method used for BSE structured eigenproblems
1484: in the Krylov-Schur solver.
1486: Not Collective
1488: Input Parameter:
1489: . eps - the linear eigensolver context
1491: Output Parameter:
1492: . bse - the BSE method
1494: Level: advanced
1496: .seealso: [](ch:eps), [](#sec:structured), `EPS_BSE`, `EPSKRYLOVSCHUR`, `EPSKrylovSchurSetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1497: @*/
1498: PetscErrorCode EPSKrylovSchurGetBSEType(EPS eps,EPSKrylovSchurBSEType *bse)
1499: {
1500: PetscFunctionBegin;
1502: PetscAssertPointer(bse,2);
1503: PetscUseMethod(eps,"EPSKrylovSchurGetBSEType_C",(EPS,EPSKrylovSchurBSEType*),(eps,bse));
1504: PetscFunctionReturn(PETSC_SUCCESS);
1505: }
1507: static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems PetscOptionsObject)
1508: {
1509: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1510: PetscBool flg,lock,b,f1,f2,f3,isfilt;
1511: PetscReal keep;
1512: PetscInt i,j,k;
1513: KSP ksp;
1514: EPSKrylovSchurBSEType bse;
1516: PetscFunctionBegin;
1517: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");
1519: PetscCall(PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg));
1520: if (flg) PetscCall(EPSKrylovSchurSetRestart(eps,keep));
1522: PetscCall(PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg));
1523: if (flg) PetscCall(EPSKrylovSchurSetLocking(eps,lock));
1525: i = ctx->npart;
1526: PetscCall(PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg));
1527: if (flg) PetscCall(EPSKrylovSchurSetPartitions(eps,i));
1529: b = ctx->detect;
1530: PetscCall(PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg));
1531: if (flg) PetscCall(EPSKrylovSchurSetDetectZeros(eps,b));
1533: i = 1;
1534: j = k = PETSC_DECIDE;
1535: PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1536: PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1537: PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1538: if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));
1540: PetscCall(PetscOptionsEnum("-eps_krylovschur_bse_type","Method for BSE structured eigenproblems","EPSKrylovSchurSetBSEType",EPSKrylovSchurBSETypes,(PetscEnum)ctx->bse,(PetscEnum*)&bse,&flg));
1541: if (flg) PetscCall(EPSKrylovSchurSetBSEType(eps,bse));
1543: PetscOptionsHeadEnd();
1545: /* set options of child KSP in spectrum slicing */
1546: if (eps->which==EPS_ALL) {
1547: if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1548: PetscCall(EPSSetDefaultST(eps));
1549: PetscCall(STSetFromOptions(eps->st)); /* need to advance this to check ST type */
1550: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1551: if (!isfilt) {
1552: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1553: PetscCall(KSPSetFromOptions(ksp));
1554: }
1555: }
1556: PetscFunctionReturn(PETSC_SUCCESS);
1557: }
1559: static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1560: {
1561: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1562: PetscBool isascii,isfilt;
1563: KSP ksp;
1564: PetscViewer sviewer;
1566: PetscFunctionBegin;
1567: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1568: if (isascii) {
1569: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1570: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1571: if (eps->problem_type==EPS_BSE) PetscCall(PetscViewerASCIIPrintf(viewer," BSE method: %s\n",EPSKrylovSchurBSETypes[ctx->bse]));
1572: if (eps->which==EPS_ALL) {
1573: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1574: if (isfilt) PetscCall(PetscViewerASCIIPrintf(viewer," using polynomial filtering%s\n",eps->nev?"":" to extract all eigenvalues in an interval"));
1575: else {
1576: PetscCall(PetscViewerASCIIPrintf(viewer," doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd));
1577: if (ctx->npart>1) {
1578: PetscCall(PetscViewerASCIIPrintf(viewer," multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart));
1579: if (ctx->detect) PetscCall(PetscViewerASCIIPrintf(viewer," detecting zeros when factorizing at subinterval boundaries\n"));
1580: }
1581: /* view child KSP */
1582: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1583: PetscCall(PetscViewerASCIIPushTab(viewer));
1584: if (ctx->npart>1 && ctx->subc) {
1585: PetscCall(PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer));
1586: if (!ctx->subc->color) PetscCall(KSPView(ksp,sviewer));
1587: PetscCall(PetscViewerFlush(sviewer));
1588: PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer));
1589: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1590: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1591: } else PetscCall(KSPView(ksp,viewer));
1592: PetscCall(PetscViewerASCIIPopTab(viewer));
1593: }
1594: }
1595: }
1596: PetscFunctionReturn(PETSC_SUCCESS);
1597: }
1599: static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1600: {
1601: PetscBool isfilt;
1603: PetscFunctionBegin;
1604: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1605: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1606: PetscCall(PetscFree(eps->data));
1607: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1608: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1609: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1610: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1611: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1612: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1613: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1614: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1615: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1616: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1617: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1618: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1619: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1620: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1621: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1622: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1623: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1624: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1625: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",NULL));
1626: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",NULL));
1627: PetscFunctionReturn(PETSC_SUCCESS);
1628: }
1630: static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1631: {
1632: PetscBool isfilt;
1634: PetscFunctionBegin;
1635: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1636: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1637: PetscFunctionReturn(PETSC_SUCCESS);
1638: }
1640: static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1641: {
1642: PetscFunctionBegin;
1643: if (eps->which==EPS_ALL) {
1644: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1645: }
1646: PetscFunctionReturn(PETSC_SUCCESS);
1647: }
1649: /*MC
1650: EPSKRYLOVSCHUR - EPSKRYLOVSCHUR = "krylovschur" - Krylov-Schur method.
1652: Notes:
1653: This is the default solver, and is recommended in most situations.
1655: The implemented algorithm is the single-vector Krylov-Schur method proposed
1656: by {cite:t}`Ste01b` for non-Hermitian eigenproblems. When the problem is
1657: Hermitian, the solver will behave as a thick-restart Lanczos method
1658: {cite:p}`Wu00` with full reorthogonalization.
1660: The solver includes support for many features\:
1661: - Harmonic extraction, see `EPSSetExtraction()`.
1662: - Inertia-based spectrum slicing to compute all eigenvalues in an interval,
1663: see {cite:p}`Cam12`.
1664: - Polynomial filter to compute eigenvalues in an interval via `STFILTER`.
1665: - Indefinite Lanczos to solve `EPS_GHIEP` problems.
1666: - Structured variants for problem types such as `EPS_BSE`.
1667: - A two-sided variant that also computes left eigenvectors.
1668: - Arbitrary selection of eigenvalues, see `EPSSetArbitrarySelection()`.
1670: Developer Note:
1671: In the future, we would like to have a block version of Krylov-Schur.
1673: Level: beginner
1675: .seealso: [](ch:eps), `EPS`, `EPSType`, `EPSSetType()`, `EPSSetProblemType()`, `EPSSetExtraction()`, `EPSSetArbitrarySelection()`, `EPSSetTwoSided()`
1676: M*/
1677: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1678: {
1679: EPS_KRYLOVSCHUR *ctx;
1681: PetscFunctionBegin;
1682: PetscCall(PetscNew(&ctx));
1683: eps->data = (void*)ctx;
1684: ctx->lock = PETSC_TRUE;
1685: ctx->nev = 1;
1686: ctx->ncv = PETSC_DETERMINE;
1687: ctx->mpd = PETSC_DETERMINE;
1688: ctx->npart = 1;
1689: ctx->detect = PETSC_FALSE;
1690: ctx->global = PETSC_TRUE;
1692: eps->useds = PETSC_TRUE;
1694: /* solve and computevectors determined at setup */
1695: eps->ops->setup = EPSSetUp_KrylovSchur;
1696: eps->ops->setupsort = EPSSetUpSort_KrylovSchur;
1697: eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1698: eps->ops->destroy = EPSDestroy_KrylovSchur;
1699: eps->ops->reset = EPSReset_KrylovSchur;
1700: eps->ops->view = EPSView_KrylovSchur;
1701: eps->ops->backtransform = EPSBackTransform_Default;
1702: eps->ops->setdefaultst = EPSSetDefaultST_KrylovSchur;
1704: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1705: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1706: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1707: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1708: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1709: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1710: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1711: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1712: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1713: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1714: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1715: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1716: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1717: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1718: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1719: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1720: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1721: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1722: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",EPSKrylovSchurSetBSEType_KrylovSchur));
1723: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",EPSKrylovSchurGetBSEType_KrylovSchur));
1724: PetscFunctionReturn(PETSC_SUCCESS);
1725: }