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