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 - Sets the restart parameter
381: Notes:
382: Allowed values are in the range [0.1,0.9]. The default is 0.5.
384: Level: advanced
386: .seealso: [](ch:eps), `EPSKrylovSchurGetRestart()`
387: @*/
388: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
389: {
390: PetscFunctionBegin;
393: PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
398: {
399: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
401: PetscFunctionBegin;
402: *keep = ctx->keep;
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: /*@
407: EPSKrylovSchurGetRestart - Gets the restart parameter used in the
408: Krylov-Schur method.
410: Not Collective
412: Input Parameter:
413: . eps - the linear eigensolver context
415: Output Parameter:
416: . keep - the restart parameter
418: Level: advanced
420: .seealso: [](ch:eps), `EPSKrylovSchurSetRestart()`
421: @*/
422: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
423: {
424: PetscFunctionBegin;
426: PetscAssertPointer(keep,2);
427: PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
428: PetscFunctionReturn(PETSC_SUCCESS);
429: }
431: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
432: {
433: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
435: PetscFunctionBegin;
436: ctx->lock = lock;
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
440: /*@
441: EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
442: the Krylov-Schur method.
444: Logically Collective
446: Input Parameters:
447: + eps - the linear eigensolver context
448: - lock - true if the locking variant must be selected
450: Options Database Key:
451: . -eps_krylovschur_locking - Sets the locking flag
453: Notes:
454: The default is to lock converged eigenpairs when the method restarts.
455: This behavior can be changed so that all directions are kept in the
456: working subspace even if already converged to working accuracy (the
457: non-locking variant).
459: Level: advanced
461: .seealso: [](ch:eps), `EPSKrylovSchurGetLocking()`
462: @*/
463: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
464: {
465: PetscFunctionBegin;
468: PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
473: {
474: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
476: PetscFunctionBegin;
477: *lock = ctx->lock;
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@
482: EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
483: method.
485: Not Collective
487: Input Parameter:
488: . eps - the linear eigensolver context
490: Output Parameter:
491: . lock - the locking flag
493: Level: advanced
495: .seealso: [](ch:eps), `EPSKrylovSchurSetLocking()`
496: @*/
497: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
498: {
499: PetscFunctionBegin;
501: PetscAssertPointer(lock,2);
502: PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
507: {
508: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
509: PetscMPIInt size;
510: PetscInt newnpart;
512: PetscFunctionBegin;
513: if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
514: newnpart = 1;
515: } else {
516: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
517: PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
518: newnpart = npart;
519: }
520: if (ctx->npart!=newnpart) {
521: if (ctx->npart>1) {
522: PetscCall(PetscSubcommDestroy(&ctx->subc));
523: if (ctx->commset) {
524: PetscCallMPI(MPI_Comm_free(&ctx->commrank));
525: ctx->commset = PETSC_FALSE;
526: }
527: }
528: PetscCall(EPSDestroy(&ctx->eps));
529: ctx->npart = newnpart;
530: eps->state = EPS_STATE_INITIAL;
531: }
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: /*@
536: EPSKrylovSchurSetPartitions - Sets the number of partitions for the
537: case of doing spectrum slicing for a computational interval with the
538: communicator split in several sub-communicators.
540: Logically Collective
542: Input Parameters:
543: + eps - the linear eigensolver context
544: - npart - number of partitions
546: Options Database Key:
547: . -eps_krylovschur_partitions <npart> - Sets the number of partitions
549: Notes:
550: By default, npart=1 so all processes in the communicator participate in
551: the processing of the whole interval. If npart>1 then the interval is
552: divided into npart subintervals, each of them being processed by a
553: subset of processes.
555: The interval is split proportionally unless the separation points are
556: specified with EPSKrylovSchurSetSubintervals().
558: Level: advanced
560: .seealso: [](ch:eps), `EPSKrylovSchurSetSubintervals()`, `EPSSetInterval()`
561: @*/
562: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
563: {
564: PetscFunctionBegin;
567: PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
572: {
573: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
575: PetscFunctionBegin;
576: *npart = ctx->npart;
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: /*@
581: EPSKrylovSchurGetPartitions - Gets the number of partitions of the
582: communicator in case of spectrum slicing.
584: Not Collective
586: Input Parameter:
587: . eps - the linear eigensolver context
589: Output Parameter:
590: . npart - number of partitions
592: Level: advanced
594: .seealso: [](ch:eps), `EPSKrylovSchurSetPartitions()`
595: @*/
596: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
597: {
598: PetscFunctionBegin;
600: PetscAssertPointer(npart,2);
601: PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
606: {
607: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
609: PetscFunctionBegin;
610: ctx->detect = detect;
611: eps->state = EPS_STATE_INITIAL;
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /*@
616: EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
617: zeros during the factorizations throughout the spectrum slicing computation.
619: Logically Collective
621: Input Parameters:
622: + eps - the linear eigensolver context
623: - detect - check for zeros
625: Options Database Key:
626: . -eps_krylovschur_detect_zeros - Check for zeros; this takes an optional
627: bool value (0/1/no/yes/true/false)
629: Notes:
630: A zero in the factorization indicates that a shift coincides with an eigenvalue.
632: This flag is turned off by default, and may be necessary in some cases,
633: especially when several partitions are being used. This feature currently
634: requires an external package for factorizations with support for zero
635: detection, e.g. MUMPS.
637: Level: advanced
639: .seealso: [](ch:eps), `EPSKrylovSchurSetPartitions()`, `EPSSetInterval()`
640: @*/
641: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
642: {
643: PetscFunctionBegin;
646: PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
651: {
652: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
654: PetscFunctionBegin;
655: *detect = ctx->detect;
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: /*@
660: EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
661: in spectrum slicing.
663: Not Collective
665: Input Parameter:
666: . eps - the linear eigensolver context
668: Output Parameter:
669: . detect - whether zeros detection is enforced during factorizations
671: Level: advanced
673: .seealso: [](ch:eps), `EPSKrylovSchurSetDetectZeros()`
674: @*/
675: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
676: {
677: PetscFunctionBegin;
679: PetscAssertPointer(detect,2);
680: PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }
684: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
685: {
686: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
688: PetscFunctionBegin;
689: if (nev != PETSC_CURRENT) {
690: PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
691: ctx->nev = nev;
692: }
693: if (ncv == PETSC_DETERMINE) {
694: ctx->ncv = PETSC_DETERMINE;
695: } else if (ncv != PETSC_CURRENT) {
696: PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
697: ctx->ncv = ncv;
698: }
699: if (mpd == PETSC_DETERMINE) {
700: ctx->mpd = PETSC_DETERMINE;
701: } else if (mpd != PETSC_CURRENT) {
702: PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
703: ctx->mpd = mpd;
704: }
705: eps->state = EPS_STATE_INITIAL;
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: /*@
710: EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
711: step in case of doing spectrum slicing for a computational interval.
712: The meaning of the parameters is the same as in EPSSetDimensions().
714: Logically Collective
716: Input Parameters:
717: + eps - the linear eigensolver context
718: . nev - number of eigenvalues to compute
719: . ncv - the maximum dimension of the subspace to be used by the subsolve
720: - mpd - the maximum dimension allowed for the projected problem
722: Options Database Keys:
723: + -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
724: . -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
725: - -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension
727: Note:
728: Use PETSC_DETERMINE for ncv and mpd to assign a default value. For any
729: of the arguments, use PETSC_CURRENT to preserve the current value.
731: Level: advanced
733: .seealso: [](ch:eps), `EPSKrylovSchurGetDimensions()`, `EPSSetDimensions()`, `EPSSetInterval()`
734: @*/
735: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
736: {
737: PetscFunctionBegin;
742: PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
747: {
748: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
750: PetscFunctionBegin;
751: if (nev) *nev = ctx->nev;
752: if (ncv) *ncv = ctx->ncv;
753: if (mpd) *mpd = ctx->mpd;
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: /*@
758: EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
759: step in case of doing spectrum slicing for a computational interval.
761: Not Collective
763: Input Parameter:
764: . eps - the linear eigensolver context
766: Output Parameters:
767: + nev - number of eigenvalues to compute
768: . ncv - the maximum dimension of the subspace to be used by the subsolve
769: - mpd - the maximum dimension allowed for the projected problem
771: Level: advanced
773: .seealso: [](ch:eps), `EPSKrylovSchurSetDimensions()`
774: @*/
775: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
776: {
777: PetscFunctionBegin;
779: PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
780: PetscFunctionReturn(PETSC_SUCCESS);
781: }
783: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
784: {
785: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
786: PetscInt i;
788: PetscFunctionBegin;
789: 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()");
790: 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");
791: if (ctx->subintervals) PetscCall(PetscFree(ctx->subintervals));
792: PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
793: for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
794: ctx->subintset = PETSC_TRUE;
795: eps->state = EPS_STATE_INITIAL;
796: PetscFunctionReturn(PETSC_SUCCESS);
797: }
799: /*@
800: EPSKrylovSchurSetSubintervals - Sets the points that delimit the
801: subintervals to be used in spectrum slicing with several partitions.
803: Logically Collective
805: Input Parameters:
806: + eps - the linear eigensolver context
807: - subint - array of real values specifying subintervals
809: Notes:
810: This function must be called after EPSKrylovSchurSetPartitions(). For npart
811: partitions, the argument subint must contain npart+1 real values sorted in
812: ascending order, subint_0, subint_1, ..., subint_npart, where the first
813: and last values must coincide with the interval endpoints set with
814: EPSSetInterval().
816: The subintervals are then defined by two consecutive points [subint_0,subint_1],
817: [subint_1,subint_2], and so on.
819: Level: advanced
821: .seealso: [](ch:eps), `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubintervals()`, `EPSSetInterval()`
822: @*/
823: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal subint[])
824: {
825: PetscFunctionBegin;
827: PetscAssertPointer(subint,2);
828: PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal *subint[])
833: {
834: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
835: PetscInt i;
837: PetscFunctionBegin;
838: if (!ctx->subintset) {
839: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
840: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
841: }
842: PetscCall(PetscMalloc1(ctx->npart+1,subint));
843: for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: /*@C
848: EPSKrylovSchurGetSubintervals - Returns the points that delimit the
849: subintervals used in spectrum slicing with several partitions.
851: Not Collective
853: Input Parameter:
854: . eps - the linear eigensolver context
856: Output Parameter:
857: . subint - array of real values specifying subintervals
859: Notes:
860: If the user passed values with EPSKrylovSchurSetSubintervals(), then the
861: same values are returned. Otherwise, the values computed internally are
862: obtained.
864: This function is only available for spectrum slicing runs.
866: The returned array has length npart+1 (see EPSKrylovSchurGetPartitions())
867: and should be freed by the user.
869: Fortran Notes:
870: The calling sequence from Fortran is
871: .vb
872: EPSKrylovSchurGetSubintervals(eps,subint,ierr)
873: double precision subint(npart+1) output
874: .ve
876: Level: advanced
878: .seealso: [](ch:eps), `EPSKrylovSchurSetSubintervals()`, `EPSKrylovSchurGetPartitions()`, `EPSSetInterval()`
879: @*/
880: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal *subint[]) PeNS
881: {
882: PetscFunctionBegin;
884: PetscAssertPointer(subint,2);
885: PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
886: PetscFunctionReturn(PETSC_SUCCESS);
887: }
889: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[])
890: {
891: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
892: PetscInt i,numsh;
893: EPS_SR sr = ctx->sr;
895: PetscFunctionBegin;
896: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
897: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
898: switch (eps->state) {
899: case EPS_STATE_INITIAL:
900: break;
901: case EPS_STATE_SETUP:
902: numsh = ctx->npart+1;
903: if (n) *n = numsh;
904: if (shifts) {
905: PetscCall(PetscMalloc1(numsh,shifts));
906: (*shifts)[0] = eps->inta;
907: if (ctx->npart==1) (*shifts)[1] = eps->intb;
908: else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
909: }
910: if (inertias) {
911: PetscCall(PetscMalloc1(numsh,inertias));
912: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
913: if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
914: else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
915: }
916: break;
917: case EPS_STATE_SOLVED:
918: case EPS_STATE_EIGENVECTORS:
919: numsh = ctx->nshifts;
920: if (n) *n = numsh;
921: if (shifts) {
922: PetscCall(PetscMalloc1(numsh,shifts));
923: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
924: }
925: if (inertias) {
926: PetscCall(PetscMalloc1(numsh,inertias));
927: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
928: }
929: break;
930: }
931: PetscFunctionReturn(PETSC_SUCCESS);
932: }
934: /*@C
935: EPSKrylovSchurGetInertias - Gets the values of the shifts and their
936: corresponding inertias in case of doing spectrum slicing for a
937: computational interval.
939: Not Collective
941: Input Parameter:
942: . eps - the linear eigensolver context
944: Output Parameters:
945: + n - number of shifts, including the endpoints of the interval
946: . shifts - the values of the shifts used internally in the solver
947: - inertias - the values of the inertia in each shift
949: Notes:
950: If called after EPSSolve(), all shifts used internally by the solver are
951: returned (including both endpoints and any intermediate ones). If called
952: before EPSSolve() and after EPSSetUp() then only the information of the
953: endpoints of subintervals is available.
955: This function is only available for spectrum slicing runs.
957: The returned arrays should be freed by the user. Can pass NULL in any of
958: the two arrays if not required.
960: Fortran Notes:
961: The calling sequence from Fortran is
962: .vb
963: EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
964: integer n
965: double precision shifts(*)
966: integer inertias(*)
967: .ve
968: The arrays should be at least of length n. The value of n can be determined
969: by an initial call
970: .vb
971: EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr)
972: .ve
974: Level: advanced
976: .seealso: [](ch:eps), `EPSSetInterval()`, `EPSKrylovSchurSetSubintervals()`
977: @*/
978: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[]) PeNS
979: {
980: PetscFunctionBegin;
982: PetscAssertPointer(n,2);
983: PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
984: PetscFunctionReturn(PETSC_SUCCESS);
985: }
987: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
988: {
989: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
990: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
992: PetscFunctionBegin;
993: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
994: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
995: if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
996: if (n) *n = sr->numEigs;
997: if (v) PetscCall(BVCreateVec(sr->V,v));
998: PetscFunctionReturn(PETSC_SUCCESS);
999: }
1001: /*@
1002: EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
1003: doing spectrum slicing for a computational interval with multiple
1004: communicators.
1006: Collective on the subcommunicator (if v is given)
1008: Input Parameter:
1009: . eps - the linear eigensolver context
1011: Output Parameters:
1012: + k - index of the subinterval for the calling process
1013: . n - number of eigenvalues found in the k-th subinterval
1014: - v - a vector owned by processes in the subcommunicator with dimensions
1015: compatible for locally computed eigenvectors (or NULL)
1017: Notes:
1018: This function is only available for spectrum slicing runs.
1020: The returned Vec should be destroyed by the user.
1022: Level: advanced
1024: .seealso: [](ch:eps), `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommPairs()`
1025: @*/
1026: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1027: {
1028: PetscFunctionBegin;
1030: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1031: PetscFunctionReturn(PETSC_SUCCESS);
1032: }
1034: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1035: {
1036: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1037: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1039: PetscFunctionBegin;
1040: EPSCheckSolved(eps,1);
1041: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1042: PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1043: if (eig) *eig = sr->eigr[sr->perm[i]];
1044: if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1045: PetscFunctionReturn(PETSC_SUCCESS);
1046: }
1048: /*@
1049: EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
1050: internally in the subcommunicator to which the calling process belongs.
1052: Collective on the subcommunicator (if v is given)
1054: Input Parameters:
1055: + eps - the linear eigensolver context
1056: - i - index of the solution
1058: Output Parameters:
1059: + eig - the eigenvalue
1060: - v - the eigenvector
1062: Notes:
1063: It is allowed to pass NULL for v if the eigenvector is not required.
1064: Otherwise, the caller must provide a valid Vec objects, i.e.,
1065: it must be created by the calling program with EPSKrylovSchurGetSubcommInfo().
1067: The index i should be a value between 0 and n-1, where n is the number of
1068: vectors in the local subinterval, see EPSKrylovSchurGetSubcommInfo().
1070: Level: advanced
1072: .seealso: [](ch:eps), `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`, `EPSKrylovSchurGetSubcommMats()`
1073: @*/
1074: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1075: {
1076: PetscFunctionBegin;
1079: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1084: {
1085: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1087: PetscFunctionBegin;
1088: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1089: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1090: PetscCall(EPSGetOperators(ctx->eps,A,B));
1091: PetscFunctionReturn(PETSC_SUCCESS);
1092: }
1094: /*@
1095: EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1096: internally in the subcommunicator to which the calling process belongs.
1098: Collective on the subcommunicator
1100: Input Parameter:
1101: . eps - the linear eigensolver context
1103: Output Parameters:
1104: + A - the matrix associated with the eigensystem
1105: - B - the second matrix in the case of generalized eigenproblems
1107: Notes:
1108: This is the analog of EPSGetOperators(), but returns the matrices distributed
1109: differently (in the subcommunicator rather than in the parent communicator).
1111: These matrices should not be modified by the user.
1113: Level: advanced
1115: .seealso: [](ch:eps), `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`
1116: @*/
1117: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1118: {
1119: PetscFunctionBegin;
1121: PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1122: PetscFunctionReturn(PETSC_SUCCESS);
1123: }
1125: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1126: {
1127: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1128: Mat A,B=NULL,Ag,Bg=NULL;
1129: PetscBool reuse=PETSC_TRUE;
1131: PetscFunctionBegin;
1132: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1133: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1134: PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1135: PetscCall(EPSGetOperators(ctx->eps,&A,&B));
1137: PetscCall(MatScale(A,a));
1138: if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1139: if (B) PetscCall(MatScale(B,b));
1140: if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1141: PetscCall(EPSSetOperators(ctx->eps,A,B));
1143: /* Update stored matrix state */
1144: subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1145: PetscCall(MatGetState(A,&subctx->Astate));
1146: if (B) PetscCall(MatGetState(B,&subctx->Bstate));
1148: /* Update matrices in the parent communicator if requested by user */
1149: if (globalup) {
1150: if (ctx->npart>1) {
1151: if (!ctx->isrow) {
1152: PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1153: reuse = PETSC_FALSE;
1154: }
1155: if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1156: if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1157: PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1158: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1159: if (B) {
1160: if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1161: PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1162: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1163: }
1164: }
1165: PetscCall(MatGetState(Ag,&ctx->Astate));
1166: if (Bg) PetscCall(MatGetState(Bg,&ctx->Bstate));
1167: }
1168: PetscCall(EPSSetOperators(eps,Ag,Bg));
1169: PetscFunctionReturn(PETSC_SUCCESS);
1170: }
1172: /*@
1173: EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1174: internally in the subcommunicator to which the calling process belongs.
1176: Collective
1178: Input Parameters:
1179: + eps - the linear eigensolver context
1180: . s - scalar that multiplies the existing A matrix
1181: . a - scalar used in the axpy operation on A
1182: . Au - matrix used in the axpy operation on A
1183: . t - scalar that multiplies the existing B matrix
1184: . b - scalar used in the axpy operation on B
1185: . Bu - matrix used in the axpy operation on B
1186: . str - structure flag
1187: - globalup - flag indicating if global matrices must be updated
1189: Notes:
1190: This function modifies the eigenproblem matrices at the subcommunicator level,
1191: and optionally updates the global matrices in the parent communicator. The updates
1192: are expressed as A <-- s*A + a*Au, B <-- t*B + b*Bu.
1194: It is possible to update one of the matrices, or both.
1196: The matrices Au and Bu must be equal in all subcommunicators.
1198: The str flag is passed to the MatAXPY() operations to perform the updates.
1200: If globalup is true, communication is carried out to reconstruct the updated
1201: matrices in the parent communicator. The user must be warned that if global
1202: matrices are not in sync with subcommunicator matrices, the errors computed
1203: by EPSComputeError() will be wrong even if the computed solution is correct
1204: (the synchronization may be done only once at the end).
1206: Level: advanced
1208: .seealso: [](ch:eps), `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommMats()`
1209: @*/
1210: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1211: {
1212: PetscFunctionBegin;
1222: PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1223: PetscFunctionReturn(PETSC_SUCCESS);
1224: }
1226: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1227: {
1228: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1229: Mat A,B=NULL,Ar=NULL,Br=NULL;
1230: PetscMPIInt rank;
1231: PetscObjectState Astate,Bstate=0;
1232: PetscObjectId Aid,Bid=0;
1233: STType sttype;
1234: PetscInt nmat;
1235: const char *prefix;
1236: MPI_Comm child;
1238: PetscFunctionBegin;
1239: PetscCall(EPSGetOperators(eps,&A,&B));
1240: if (ctx->npart==1) {
1241: if (!ctx->eps) PetscCall(EPSCreate(((PetscObject)eps)->comm,&ctx->eps));
1242: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1243: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1244: PetscCall(EPSSetOperators(ctx->eps,A,B));
1245: } else {
1246: PetscCall(MatGetState(A,&Astate));
1247: PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1248: if (B) {
1249: PetscCall(MatGetState(B,&Bstate));
1250: PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1251: }
1252: if (!ctx->subc) {
1253: /* Create context for subcommunicators */
1254: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1255: PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1256: PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1257: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1259: /* Duplicate matrices */
1260: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1261: ctx->Astate = Astate;
1262: ctx->Aid = Aid;
1263: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1264: if (B) {
1265: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1266: ctx->Bstate = Bstate;
1267: ctx->Bid = Bid;
1268: PetscCall(MatPropagateSymmetryOptions(B,Br));
1269: }
1270: } else {
1271: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1272: if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1273: PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1274: if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
1275: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1276: ctx->Astate = Astate;
1277: ctx->Aid = Aid;
1278: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1279: if (B) {
1280: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1281: ctx->Bstate = Bstate;
1282: ctx->Bid = Bid;
1283: PetscCall(MatPropagateSymmetryOptions(B,Br));
1284: }
1285: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1286: PetscCall(MatDestroy(&Ar));
1287: PetscCall(MatDestroy(&Br));
1288: }
1289: }
1291: /* Create auxiliary EPS */
1292: if (!ctx->eps) {
1293: PetscCall(EPSCreate(child,&ctx->eps));
1294: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1295: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1296: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1297: PetscCall(MatDestroy(&Ar));
1298: PetscCall(MatDestroy(&Br));
1299: }
1300: /* Create subcommunicator grouping processes with same rank */
1301: if (!ctx->commset) {
1302: PetscCallMPI(MPI_Comm_rank(child,&rank));
1303: PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1304: ctx->commset = PETSC_TRUE;
1305: }
1306: }
1307: PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1308: PetscCall(STGetType(eps->st,&sttype));
1309: PetscCall(STSetType(ctx->eps->st,sttype));
1311: ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1312: ctx_local->npart = ctx->npart;
1313: ctx_local->global = PETSC_FALSE;
1314: ctx_local->eps = eps;
1315: ctx_local->subc = ctx->subc;
1316: ctx_local->commrank = ctx->commrank;
1317: *childeps = ctx->eps;
1318: PetscFunctionReturn(PETSC_SUCCESS);
1319: }
1321: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1322: {
1323: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1324: ST st;
1325: PetscBool isfilt;
1327: PetscFunctionBegin;
1328: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1329: PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1330: PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1331: PetscCall(EPSGetST(ctx->eps,&st));
1332: PetscCall(STGetOperator(st,NULL));
1333: PetscCall(STGetKSP(st,ksp));
1334: PetscFunctionReturn(PETSC_SUCCESS);
1335: }
1337: /*@
1338: EPSKrylovSchurGetKSP - Retrieve the linear solver object associated with the
1339: internal EPS object in case of doing spectrum slicing for a computational interval.
1341: Collective
1343: Input Parameter:
1344: . eps - the linear eigensolver context
1346: Output Parameter:
1347: . ksp - the internal KSP object
1349: Notes:
1350: When invoked to compute all eigenvalues in an interval with spectrum
1351: slicing, EPSKRYLOVSCHUR creates another EPS object internally that is
1352: used to compute eigenvalues by chunks near selected shifts. This function
1353: allows access to the KSP object associated to this internal EPS object.
1355: This function is only available for spectrum slicing runs. In case of
1356: having more than one partition, the returned KSP will be different
1357: in MPI processes belonging to different partitions. Hence, if required,
1358: EPSKrylovSchurSetPartitions() must be called BEFORE this function.
1360: Level: advanced
1362: .seealso: [](ch:eps), `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`
1363: @*/
1364: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1365: {
1366: PetscFunctionBegin;
1368: PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1369: PetscFunctionReturn(PETSC_SUCCESS);
1370: }
1372: static PetscErrorCode EPSKrylovSchurSetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType bse)
1373: {
1374: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1376: PetscFunctionBegin;
1377: switch (bse) {
1378: case EPS_KRYLOVSCHUR_BSE_SHAO:
1379: case EPS_KRYLOVSCHUR_BSE_GRUNING:
1380: case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
1381: if (ctx->bse != bse) {
1382: ctx->bse = bse;
1383: eps->state = EPS_STATE_INITIAL;
1384: }
1385: break;
1386: default:
1387: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid BSE type");
1388: }
1389: PetscFunctionReturn(PETSC_SUCCESS);
1390: }
1392: /*@
1393: EPSKrylovSchurSetBSEType - Sets the method to be used for BSE structured eigenproblems
1394: in the Krylov-Schur solver.
1396: Logically Collective
1398: Input Parameters:
1399: + eps - the linear eigensolver context
1400: - bse - the BSE method
1402: Options Database Key:
1403: . -eps_krylovschur_bse_type - Sets the BSE type (either 'shao', 'gruning', or 'projectedbse')
1405: Level: advanced
1407: .seealso: [](ch:eps), `EPSKrylovSchurGetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1408: @*/
1409: PetscErrorCode EPSKrylovSchurSetBSEType(EPS eps,EPSKrylovSchurBSEType bse)
1410: {
1411: PetscFunctionBegin;
1414: PetscTryMethod(eps,"EPSKrylovSchurSetBSEType_C",(EPS,EPSKrylovSchurBSEType),(eps,bse));
1415: PetscFunctionReturn(PETSC_SUCCESS);
1416: }
1418: static PetscErrorCode EPSKrylovSchurGetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType *bse)
1419: {
1420: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1422: PetscFunctionBegin;
1423: *bse = ctx->bse;
1424: PetscFunctionReturn(PETSC_SUCCESS);
1425: }
1427: /*@
1428: EPSKrylovSchurGetBSEType - Gets the method used for BSE structured eigenproblems
1429: in the Krylov-Schur solver.
1431: Not Collective
1433: Input Parameter:
1434: . eps - the linear eigensolver context
1436: Output Parameter:
1437: . bse - the BSE method
1439: Level: advanced
1441: .seealso: [](ch:eps), `EPSKrylovSchurSetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1442: @*/
1443: PetscErrorCode EPSKrylovSchurGetBSEType(EPS eps,EPSKrylovSchurBSEType *bse)
1444: {
1445: PetscFunctionBegin;
1447: PetscAssertPointer(bse,2);
1448: PetscUseMethod(eps,"EPSKrylovSchurGetBSEType_C",(EPS,EPSKrylovSchurBSEType*),(eps,bse));
1449: PetscFunctionReturn(PETSC_SUCCESS);
1450: }
1452: static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems PetscOptionsObject)
1453: {
1454: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1455: PetscBool flg,lock,b,f1,f2,f3,isfilt;
1456: PetscReal keep;
1457: PetscInt i,j,k;
1458: KSP ksp;
1459: EPSKrylovSchurBSEType bse;
1461: PetscFunctionBegin;
1462: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");
1464: PetscCall(PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg));
1465: if (flg) PetscCall(EPSKrylovSchurSetRestart(eps,keep));
1467: PetscCall(PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg));
1468: if (flg) PetscCall(EPSKrylovSchurSetLocking(eps,lock));
1470: i = ctx->npart;
1471: PetscCall(PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg));
1472: if (flg) PetscCall(EPSKrylovSchurSetPartitions(eps,i));
1474: b = ctx->detect;
1475: PetscCall(PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg));
1476: if (flg) PetscCall(EPSKrylovSchurSetDetectZeros(eps,b));
1478: i = 1;
1479: j = k = PETSC_DECIDE;
1480: PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1481: PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1482: PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1483: if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));
1485: PetscCall(PetscOptionsEnum("-eps_krylovschur_bse_type","Method for BSE structured eigenproblems","EPSKrylovSchurSetBSEType",EPSKrylovSchurBSETypes,(PetscEnum)ctx->bse,(PetscEnum*)&bse,&flg));
1486: if (flg) PetscCall(EPSKrylovSchurSetBSEType(eps,bse));
1488: PetscOptionsHeadEnd();
1490: /* set options of child KSP in spectrum slicing */
1491: if (eps->which==EPS_ALL) {
1492: if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1493: PetscCall(EPSSetDefaultST(eps));
1494: PetscCall(STSetFromOptions(eps->st)); /* need to advance this to check ST type */
1495: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1496: if (!isfilt) {
1497: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1498: PetscCall(KSPSetFromOptions(ksp));
1499: }
1500: }
1501: PetscFunctionReturn(PETSC_SUCCESS);
1502: }
1504: static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1505: {
1506: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1507: PetscBool isascii,isfilt;
1508: KSP ksp;
1509: PetscViewer sviewer;
1511: PetscFunctionBegin;
1512: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1513: if (isascii) {
1514: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1515: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1516: if (eps->problem_type==EPS_BSE) PetscCall(PetscViewerASCIIPrintf(viewer," BSE method: %s\n",EPSKrylovSchurBSETypes[ctx->bse]));
1517: if (eps->which==EPS_ALL) {
1518: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1519: if (isfilt) PetscCall(PetscViewerASCIIPrintf(viewer," using filtering to extract all eigenvalues in an interval\n"));
1520: else {
1521: PetscCall(PetscViewerASCIIPrintf(viewer," doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd));
1522: if (ctx->npart>1) {
1523: PetscCall(PetscViewerASCIIPrintf(viewer," multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart));
1524: if (ctx->detect) PetscCall(PetscViewerASCIIPrintf(viewer," detecting zeros when factorizing at subinterval boundaries\n"));
1525: }
1526: /* view child KSP */
1527: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1528: PetscCall(PetscViewerASCIIPushTab(viewer));
1529: if (ctx->npart>1 && ctx->subc) {
1530: PetscCall(PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer));
1531: if (!ctx->subc->color) PetscCall(KSPView(ksp,sviewer));
1532: PetscCall(PetscViewerFlush(sviewer));
1533: PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer));
1534: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1535: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1536: } else PetscCall(KSPView(ksp,viewer));
1537: PetscCall(PetscViewerASCIIPopTab(viewer));
1538: }
1539: }
1540: }
1541: PetscFunctionReturn(PETSC_SUCCESS);
1542: }
1544: static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1545: {
1546: PetscBool isfilt;
1548: PetscFunctionBegin;
1549: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1550: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1551: PetscCall(PetscFree(eps->data));
1552: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1553: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1554: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1555: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1556: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1557: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1558: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1559: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1560: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1561: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1562: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1563: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1564: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1565: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1566: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1567: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1568: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1569: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1570: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",NULL));
1571: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",NULL));
1572: PetscFunctionReturn(PETSC_SUCCESS);
1573: }
1575: static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1576: {
1577: PetscBool isfilt;
1579: PetscFunctionBegin;
1580: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1581: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1582: PetscFunctionReturn(PETSC_SUCCESS);
1583: }
1585: static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1586: {
1587: PetscFunctionBegin;
1588: if (eps->which==EPS_ALL) {
1589: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1590: }
1591: PetscFunctionReturn(PETSC_SUCCESS);
1592: }
1594: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1595: {
1596: EPS_KRYLOVSCHUR *ctx;
1598: PetscFunctionBegin;
1599: PetscCall(PetscNew(&ctx));
1600: eps->data = (void*)ctx;
1601: ctx->lock = PETSC_TRUE;
1602: ctx->nev = 1;
1603: ctx->ncv = PETSC_DETERMINE;
1604: ctx->mpd = PETSC_DETERMINE;
1605: ctx->npart = 1;
1606: ctx->detect = PETSC_FALSE;
1607: ctx->global = PETSC_TRUE;
1609: eps->useds = PETSC_TRUE;
1611: /* solve and computevectors determined at setup */
1612: eps->ops->setup = EPSSetUp_KrylovSchur;
1613: eps->ops->setupsort = EPSSetUpSort_KrylovSchur;
1614: eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1615: eps->ops->destroy = EPSDestroy_KrylovSchur;
1616: eps->ops->reset = EPSReset_KrylovSchur;
1617: eps->ops->view = EPSView_KrylovSchur;
1618: eps->ops->backtransform = EPSBackTransform_Default;
1619: eps->ops->setdefaultst = EPSSetDefaultST_KrylovSchur;
1621: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1622: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1623: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1624: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1625: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1626: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1627: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1628: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1629: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1630: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1631: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1632: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1633: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1634: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1635: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1636: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1637: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1638: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1639: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",EPSKrylovSchurSetBSEType_KrylovSchur));
1640: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",EPSKrylovSchurGetBSEType_KrylovSchur));
1641: PetscFunctionReturn(PETSC_SUCCESS);
1642: }