Actual source code: krylovschur.c
slepc-main 2024-11-09
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: SLEPc eigensolver: "krylovschur"
13: Method: Krylov-Schur
15: Algorithm:
17: Single-vector Krylov-Schur method for non-symmetric problems,
18: including harmonic extraction.
20: References:
22: [1] "Krylov-Schur Methods in SLEPc", SLEPc Technical Report STR-7,
23: available at https://slepc.upv.es.
25: [2] G.W. Stewart, "A Krylov-Schur Algorithm for Large Eigenproblems",
26: SIAM J. Matrix Anal. App. 23(3):601-614, 2001.
28: [3] "Practical Implementation of Harmonic Krylov-Schur", SLEPc Technical
29: Report STR-9, available at https://slepc.upv.es.
30: */
32: #include <slepc/private/epsimpl.h>
33: #include "krylovschur.h"
35: PetscErrorCode EPSGetArbitraryValues(EPS eps,PetscScalar *rr,PetscScalar *ri)
36: {
37: PetscInt i,newi,ld,n,l;
38: Vec xr=eps->work[0],xi=eps->work[1];
39: PetscScalar re,im,*Zr,*Zi,*X;
41: PetscFunctionBegin;
42: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
43: PetscCall(DSGetDimensions(eps->ds,&n,&l,NULL,NULL));
44: for (i=l;i<n;i++) {
45: re = eps->eigr[i];
46: im = eps->eigi[i];
47: PetscCall(STBackTransform(eps->st,1,&re,&im));
48: newi = i;
49: PetscCall(DSVectors(eps->ds,DS_MAT_X,&newi,NULL));
50: PetscCall(DSGetArray(eps->ds,DS_MAT_X,&X));
51: Zr = X+i*ld;
52: if (newi==i+1) Zi = X+newi*ld;
53: else Zi = NULL;
54: PetscCall(EPSComputeRitzVector(eps,Zr,Zi,eps->V,xr,xi));
55: PetscCall(DSRestoreArray(eps->ds,DS_MAT_X,&X));
56: PetscCall((*eps->arbitrary)(re,im,xr,xi,rr+i,ri+i,eps->arbitraryctx));
57: }
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode EPSSetUp_KrylovSchur_Filter(EPS eps)
62: {
63: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
64: PetscBool estimaterange=PETSC_TRUE;
65: PetscReal rleft,rright;
66: Mat A;
68: PetscFunctionBegin;
69: EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
70: EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
71: PetscCheck(eps->intb<PETSC_MAX_REAL || eps->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The defined computational interval should have at least one of their sides bounded");
72: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
73: 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==1) 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);
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 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: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
286: nconv = k;
288: /* Update l */
289: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
290: else {
291: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
292: if (!hermitian) PetscCall(DSGetTruncateSize(eps->ds,k,nv,&l));
293: }
294: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
295: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
297: if (eps->reason == EPS_CONVERGED_ITERATING) {
298: if (PetscUnlikely(breakdown || k==nv)) {
299: /* Start a new Arnoldi factorization */
300: PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
301: if (k<eps->nev) {
302: PetscCall(EPSGetStartVector(eps,k,&breakdown));
303: if (breakdown) {
304: eps->reason = EPS_DIVERGED_BREAKDOWN;
305: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
306: }
307: }
308: } else {
309: /* Undo translation of Krylov decomposition */
310: if (PetscUnlikely(harmonic)) {
311: PetscCall(DSSetDimensions(eps->ds,nv,k,l));
312: PetscCall(DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma));
313: /* gamma u^ = u - U*g~ */
314: PetscCall(BVSetActiveColumns(eps->V,0,nv));
315: PetscCall(BVMultColumn(eps->V,-1.0,1.0,nv,g));
316: PetscCall(BVScaleColumn(eps->V,nv,1.0/gamma));
317: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
318: PetscCall(DSSetDimensions(eps->ds,nv,k,nv));
319: }
320: /* Prepare the Rayleigh quotient for restart */
321: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
322: }
323: }
324: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
325: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
326: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
327: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
329: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(eps->V,nv,k+l));
330: eps->nconv = k;
331: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
332: }
334: if (harmonic) PetscCall(PetscFree(g));
335: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
340: {
341: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
343: PetscFunctionBegin;
344: if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
345: else {
346: 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);
347: ctx->keep = keep;
348: }
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: /*@
353: EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
354: method, in particular the proportion of basis vectors that must be kept
355: after restart.
357: Logically Collective
359: Input Parameters:
360: + eps - the eigenproblem solver context
361: - keep - the number of vectors to be kept at restart
363: Options Database Key:
364: . -eps_krylovschur_restart - Sets the restart parameter
366: Notes:
367: Allowed values are in the range [0.1,0.9]. The default is 0.5.
369: Level: advanced
371: .seealso: EPSKrylovSchurGetRestart()
372: @*/
373: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
374: {
375: PetscFunctionBegin;
378: PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
379: PetscFunctionReturn(PETSC_SUCCESS);
380: }
382: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
383: {
384: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
386: PetscFunctionBegin;
387: *keep = ctx->keep;
388: PetscFunctionReturn(PETSC_SUCCESS);
389: }
391: /*@
392: EPSKrylovSchurGetRestart - Gets the restart parameter used in the
393: Krylov-Schur method.
395: Not Collective
397: Input Parameter:
398: . eps - the eigenproblem solver context
400: Output Parameter:
401: . keep - the restart parameter
403: Level: advanced
405: .seealso: EPSKrylovSchurSetRestart()
406: @*/
407: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
408: {
409: PetscFunctionBegin;
411: PetscAssertPointer(keep,2);
412: PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
417: {
418: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
420: PetscFunctionBegin;
421: ctx->lock = lock;
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*@
426: EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
427: the Krylov-Schur method.
429: Logically Collective
431: Input Parameters:
432: + eps - the eigenproblem solver context
433: - lock - true if the locking variant must be selected
435: Options Database Key:
436: . -eps_krylovschur_locking - Sets the locking flag
438: Notes:
439: The default is to lock converged eigenpairs when the method restarts.
440: This behaviour can be changed so that all directions are kept in the
441: working subspace even if already converged to working accuracy (the
442: non-locking variant).
444: Level: advanced
446: .seealso: EPSKrylovSchurGetLocking()
447: @*/
448: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
449: {
450: PetscFunctionBegin;
453: PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
458: {
459: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
461: PetscFunctionBegin;
462: *lock = ctx->lock;
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: /*@
467: EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
468: method.
470: Not Collective
472: Input Parameter:
473: . eps - the eigenproblem solver context
475: Output Parameter:
476: . lock - the locking flag
478: Level: advanced
480: .seealso: EPSKrylovSchurSetLocking()
481: @*/
482: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
483: {
484: PetscFunctionBegin;
486: PetscAssertPointer(lock,2);
487: PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }
491: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
492: {
493: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
494: PetscMPIInt size;
495: PetscInt newnpart;
497: PetscFunctionBegin;
498: if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
499: newnpart = 1;
500: } else {
501: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
502: PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
503: newnpart = npart;
504: }
505: if (ctx->npart!=newnpart) {
506: if (ctx->npart>1) {
507: PetscCall(PetscSubcommDestroy(&ctx->subc));
508: if (ctx->commset) {
509: PetscCallMPI(MPI_Comm_free(&ctx->commrank));
510: ctx->commset = PETSC_FALSE;
511: }
512: }
513: PetscCall(EPSDestroy(&ctx->eps));
514: ctx->npart = newnpart;
515: eps->state = EPS_STATE_INITIAL;
516: }
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: /*@
521: EPSKrylovSchurSetPartitions - Sets the number of partitions for the
522: case of doing spectrum slicing for a computational interval with the
523: communicator split in several sub-communicators.
525: Logically Collective
527: Input Parameters:
528: + eps - the eigenproblem solver context
529: - npart - number of partitions
531: Options Database Key:
532: . -eps_krylovschur_partitions <npart> - Sets the number of partitions
534: Notes:
535: By default, npart=1 so all processes in the communicator participate in
536: the processing of the whole interval. If npart>1 then the interval is
537: divided into npart subintervals, each of them being processed by a
538: subset of processes.
540: The interval is split proportionally unless the separation points are
541: specified with EPSKrylovSchurSetSubintervals().
543: Level: advanced
545: .seealso: EPSKrylovSchurSetSubintervals(), EPSSetInterval()
546: @*/
547: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
548: {
549: PetscFunctionBegin;
552: PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
557: {
558: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
560: PetscFunctionBegin;
561: *npart = ctx->npart;
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: /*@
566: EPSKrylovSchurGetPartitions - Gets the number of partitions of the
567: communicator in case of spectrum slicing.
569: Not Collective
571: Input Parameter:
572: . eps - the eigenproblem solver context
574: Output Parameter:
575: . npart - number of partitions
577: Level: advanced
579: .seealso: EPSKrylovSchurSetPartitions()
580: @*/
581: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
582: {
583: PetscFunctionBegin;
585: PetscAssertPointer(npart,2);
586: PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
591: {
592: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
594: PetscFunctionBegin;
595: ctx->detect = detect;
596: eps->state = EPS_STATE_INITIAL;
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: /*@
601: EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
602: zeros during the factorizations throughout the spectrum slicing computation.
604: Logically Collective
606: Input Parameters:
607: + eps - the eigenproblem solver context
608: - detect - check for zeros
610: Options Database Key:
611: . -eps_krylovschur_detect_zeros - Check for zeros; this takes an optional
612: bool value (0/1/no/yes/true/false)
614: Notes:
615: A zero in the factorization indicates that a shift coincides with an eigenvalue.
617: This flag is turned off by default, and may be necessary in some cases,
618: especially when several partitions are being used. This feature currently
619: requires an external package for factorizations with support for zero
620: detection, e.g. MUMPS.
622: Level: advanced
624: .seealso: EPSKrylovSchurSetPartitions(), EPSSetInterval()
625: @*/
626: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
627: {
628: PetscFunctionBegin;
631: PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
636: {
637: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
639: PetscFunctionBegin;
640: *detect = ctx->detect;
641: PetscFunctionReturn(PETSC_SUCCESS);
642: }
644: /*@
645: EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
646: in spectrum slicing.
648: Not Collective
650: Input Parameter:
651: . eps - the eigenproblem solver context
653: Output Parameter:
654: . detect - whether zeros detection is enforced during factorizations
656: Level: advanced
658: .seealso: EPSKrylovSchurSetDetectZeros()
659: @*/
660: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
661: {
662: PetscFunctionBegin;
664: PetscAssertPointer(detect,2);
665: PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
670: {
671: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
673: PetscFunctionBegin;
674: if (nev != PETSC_CURRENT) {
675: PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
676: ctx->nev = nev;
677: }
678: if (ncv == PETSC_DETERMINE) {
679: ctx->ncv = PETSC_DETERMINE;
680: } else if (ncv != PETSC_CURRENT) {
681: PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
682: ctx->ncv = ncv;
683: }
684: if (mpd == PETSC_DETERMINE) {
685: ctx->mpd = PETSC_DETERMINE;
686: } else if (mpd != PETSC_CURRENT) {
687: PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
688: ctx->mpd = mpd;
689: }
690: eps->state = EPS_STATE_INITIAL;
691: PetscFunctionReturn(PETSC_SUCCESS);
692: }
694: /*@
695: EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
696: step in case of doing spectrum slicing for a computational interval.
697: The meaning of the parameters is the same as in EPSSetDimensions().
699: Logically Collective
701: Input Parameters:
702: + eps - the eigenproblem solver context
703: . nev - number of eigenvalues to compute
704: . ncv - the maximum dimension of the subspace to be used by the subsolve
705: - mpd - the maximum dimension allowed for the projected problem
707: Options Database Key:
708: + -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
709: . -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
710: - -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension
712: Note:
713: Use PETSC_DETERMINE for ncv and mpd to assign a default value. For any
714: of the arguments, use PETSC_CURRENT to preserve the current value.
716: Level: advanced
718: .seealso: EPSKrylovSchurGetDimensions(), EPSSetDimensions(), EPSSetInterval()
719: @*/
720: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
721: {
722: PetscFunctionBegin;
727: PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
732: {
733: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
735: PetscFunctionBegin;
736: if (nev) *nev = ctx->nev;
737: if (ncv) *ncv = ctx->ncv;
738: if (mpd) *mpd = ctx->mpd;
739: PetscFunctionReturn(PETSC_SUCCESS);
740: }
742: /*@
743: EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
744: step in case of doing spectrum slicing for a computational interval.
746: Not Collective
748: Input Parameter:
749: . eps - the eigenproblem solver context
751: Output Parameters:
752: + nev - number of eigenvalues to compute
753: . ncv - the maximum dimension of the subspace to be used by the subsolve
754: - mpd - the maximum dimension allowed for the projected problem
756: Level: advanced
758: .seealso: EPSKrylovSchurSetDimensions()
759: @*/
760: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
761: {
762: PetscFunctionBegin;
764: PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
765: PetscFunctionReturn(PETSC_SUCCESS);
766: }
768: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
769: {
770: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
771: PetscInt i;
773: PetscFunctionBegin;
774: 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()");
775: 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");
776: if (ctx->subintervals) PetscCall(PetscFree(ctx->subintervals));
777: PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
778: for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
779: ctx->subintset = PETSC_TRUE;
780: eps->state = EPS_STATE_INITIAL;
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: /*@
785: EPSKrylovSchurSetSubintervals - Sets the points that delimit the
786: subintervals to be used in spectrum slicing with several partitions.
788: Logically Collective
790: Input Parameters:
791: + eps - the eigenproblem solver context
792: - subint - array of real values specifying subintervals
794: Notes:
795: This function must be called after EPSKrylovSchurSetPartitions(). For npart
796: partitions, the argument subint must contain npart+1 real values sorted in
797: ascending order, subint_0, subint_1, ..., subint_npart, where the first
798: and last values must coincide with the interval endpoints set with
799: EPSSetInterval().
801: The subintervals are then defined by two consecutive points [subint_0,subint_1],
802: [subint_1,subint_2], and so on.
804: Level: advanced
806: .seealso: EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubintervals(), EPSSetInterval()
807: @*/
808: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal subint[])
809: {
810: PetscFunctionBegin;
812: PetscAssertPointer(subint,2);
813: PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal **subint)
818: {
819: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
820: PetscInt i;
822: PetscFunctionBegin;
823: if (!ctx->subintset) {
824: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
825: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
826: }
827: PetscCall(PetscMalloc1(ctx->npart+1,subint));
828: for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: /*@C
833: EPSKrylovSchurGetSubintervals - Returns the points that delimit the
834: subintervals used in spectrum slicing with several partitions.
836: Not Collective
838: Input Parameter:
839: . eps - the eigenproblem solver context
841: Output Parameter:
842: . subint - array of real values specifying subintervals
844: Notes:
845: If the user passed values with EPSKrylovSchurSetSubintervals(), then the
846: same values are returned. Otherwise, the values computed internally are
847: obtained.
849: This function is only available for spectrum slicing runs.
851: The returned array has length npart+1 (see EPSKrylovSchurGetPartitions())
852: and should be freed by the user.
854: Fortran Notes:
855: The calling sequence from Fortran is
856: .vb
857: EPSKrylovSchurGetSubintervals(eps,subint,ierr)
858: double precision subint(npart+1) output
859: .ve
861: Level: advanced
863: .seealso: EPSKrylovSchurSetSubintervals(), EPSKrylovSchurGetPartitions(), EPSSetInterval()
864: @*/
865: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal **subint)
866: {
867: PetscFunctionBegin;
869: PetscAssertPointer(subint,2);
870: PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
875: {
876: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
877: PetscInt i,numsh;
878: EPS_SR sr = ctx->sr;
880: PetscFunctionBegin;
881: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
882: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
883: switch (eps->state) {
884: case EPS_STATE_INITIAL:
885: break;
886: case EPS_STATE_SETUP:
887: numsh = ctx->npart+1;
888: if (n) *n = numsh;
889: if (shifts) {
890: PetscCall(PetscMalloc1(numsh,shifts));
891: (*shifts)[0] = eps->inta;
892: if (ctx->npart==1) (*shifts)[1] = eps->intb;
893: else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
894: }
895: if (inertias) {
896: PetscCall(PetscMalloc1(numsh,inertias));
897: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
898: if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
899: else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
900: }
901: break;
902: case EPS_STATE_SOLVED:
903: case EPS_STATE_EIGENVECTORS:
904: numsh = ctx->nshifts;
905: if (n) *n = numsh;
906: if (shifts) {
907: PetscCall(PetscMalloc1(numsh,shifts));
908: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
909: }
910: if (inertias) {
911: PetscCall(PetscMalloc1(numsh,inertias));
912: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
913: }
914: break;
915: }
916: PetscFunctionReturn(PETSC_SUCCESS);
917: }
919: /*@C
920: EPSKrylovSchurGetInertias - Gets the values of the shifts and their
921: corresponding inertias in case of doing spectrum slicing for a
922: computational interval.
924: Not Collective
926: Input Parameter:
927: . eps - the eigenproblem solver context
929: Output Parameters:
930: + n - number of shifts, including the endpoints of the interval
931: . shifts - the values of the shifts used internally in the solver
932: - inertias - the values of the inertia in each shift
934: Notes:
935: If called after EPSSolve(), all shifts used internally by the solver are
936: returned (including both endpoints and any intermediate ones). If called
937: before EPSSolve() and after EPSSetUp() then only the information of the
938: endpoints of subintervals is available.
940: This function is only available for spectrum slicing runs.
942: The returned arrays should be freed by the user. Can pass NULL in any of
943: the two arrays if not required.
945: Fortran Notes:
946: The calling sequence from Fortran is
947: .vb
948: EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
949: integer n
950: double precision shifts(*)
951: integer inertias(*)
952: .ve
953: The arrays should be at least of length n. The value of n can be determined
954: by an initial call
955: .vb
956: EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr)
957: .ve
959: Level: advanced
961: .seealso: EPSSetInterval(), EPSKrylovSchurSetSubintervals()
962: @*/
963: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[])
964: {
965: PetscFunctionBegin;
967: PetscAssertPointer(n,2);
968: PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
969: PetscFunctionReturn(PETSC_SUCCESS);
970: }
972: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
973: {
974: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
975: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
977: PetscFunctionBegin;
978: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
979: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
980: if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
981: if (n) *n = sr->numEigs;
982: if (v) PetscCall(BVCreateVec(sr->V,v));
983: PetscFunctionReturn(PETSC_SUCCESS);
984: }
986: /*@
987: EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
988: doing spectrum slicing for a computational interval with multiple
989: communicators.
991: Collective on the subcommunicator (if v is given)
993: Input Parameter:
994: . eps - the eigenproblem solver context
996: Output Parameters:
997: + k - index of the subinterval for the calling process
998: . n - number of eigenvalues found in the k-th subinterval
999: - v - a vector owned by processes in the subcommunicator with dimensions
1000: compatible for locally computed eigenvectors (or NULL)
1002: Notes:
1003: This function is only available for spectrum slicing runs.
1005: The returned Vec should be destroyed by the user.
1007: Level: advanced
1009: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommPairs()
1010: @*/
1011: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1012: {
1013: PetscFunctionBegin;
1015: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1016: PetscFunctionReturn(PETSC_SUCCESS);
1017: }
1019: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1020: {
1021: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1022: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1024: PetscFunctionBegin;
1025: EPSCheckSolved(eps,1);
1026: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1027: PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1028: if (eig) *eig = sr->eigr[sr->perm[i]];
1029: if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: /*@
1034: EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
1035: internally in the subcommunicator to which the calling process belongs.
1037: Collective on the subcommunicator (if v is given)
1039: Input Parameters:
1040: + eps - the eigenproblem solver context
1041: - i - index of the solution
1043: Output Parameters:
1044: + eig - the eigenvalue
1045: - v - the eigenvector
1047: Notes:
1048: It is allowed to pass NULL for v if the eigenvector is not required.
1049: Otherwise, the caller must provide a valid Vec objects, i.e.,
1050: it must be created by the calling program with EPSKrylovSchurGetSubcommInfo().
1052: The index i should be a value between 0 and n-1, where n is the number of
1053: vectors in the local subinterval, see EPSKrylovSchurGetSubcommInfo().
1055: Level: advanced
1057: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo(), EPSKrylovSchurGetSubcommMats()
1058: @*/
1059: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1060: {
1061: PetscFunctionBegin;
1064: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1065: PetscFunctionReturn(PETSC_SUCCESS);
1066: }
1068: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1069: {
1070: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1072: PetscFunctionBegin;
1073: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1074: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1075: PetscCall(EPSGetOperators(ctx->eps,A,B));
1076: PetscFunctionReturn(PETSC_SUCCESS);
1077: }
1079: /*@
1080: EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1081: internally in the subcommunicator to which the calling process belongs.
1083: Collective on the subcommunicator
1085: Input Parameter:
1086: . eps - the eigenproblem solver context
1088: Output Parameters:
1089: + A - the matrix associated with the eigensystem
1090: - B - the second matrix in the case of generalized eigenproblems
1092: Notes:
1093: This is the analog of EPSGetOperators(), but returns the matrices distributed
1094: differently (in the subcommunicator rather than in the parent communicator).
1096: These matrices should not be modified by the user.
1098: Level: advanced
1100: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo()
1101: @*/
1102: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1103: {
1104: PetscFunctionBegin;
1106: PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1107: PetscFunctionReturn(PETSC_SUCCESS);
1108: }
1110: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1111: {
1112: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1113: Mat A,B=NULL,Ag,Bg=NULL;
1114: PetscBool reuse=PETSC_TRUE;
1116: PetscFunctionBegin;
1117: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1118: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1119: PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1120: PetscCall(EPSGetOperators(ctx->eps,&A,&B));
1122: PetscCall(MatScale(A,a));
1123: if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1124: if (B) PetscCall(MatScale(B,b));
1125: if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1126: PetscCall(EPSSetOperators(ctx->eps,A,B));
1128: /* Update stored matrix state */
1129: subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1130: PetscCall(MatGetState(A,&subctx->Astate));
1131: if (B) PetscCall(MatGetState(B,&subctx->Bstate));
1133: /* Update matrices in the parent communicator if requested by user */
1134: if (globalup) {
1135: if (ctx->npart>1) {
1136: if (!ctx->isrow) {
1137: PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1138: reuse = PETSC_FALSE;
1139: }
1140: if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1141: if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1142: PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1143: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1144: if (B) {
1145: if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1146: PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1147: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1148: }
1149: }
1150: PetscCall(MatGetState(Ag,&ctx->Astate));
1151: if (Bg) PetscCall(MatGetState(Bg,&ctx->Bstate));
1152: }
1153: PetscCall(EPSSetOperators(eps,Ag,Bg));
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: /*@
1158: EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1159: internally in the subcommunicator to which the calling process belongs.
1161: Collective
1163: Input Parameters:
1164: + eps - the eigenproblem solver context
1165: . s - scalar that multiplies the existing A matrix
1166: . a - scalar used in the axpy operation on A
1167: . Au - matrix used in the axpy operation on A
1168: . t - scalar that multiplies the existing B matrix
1169: . b - scalar used in the axpy operation on B
1170: . Bu - matrix used in the axpy operation on B
1171: . str - structure flag
1172: - globalup - flag indicating if global matrices must be updated
1174: Notes:
1175: This function modifies the eigenproblem matrices at the subcommunicator level,
1176: and optionally updates the global matrices in the parent communicator. The updates
1177: are expressed as A <-- s*A + a*Au, B <-- t*B + b*Bu.
1179: It is possible to update one of the matrices, or both.
1181: The matrices Au and Bu must be equal in all subcommunicators.
1183: The str flag is passed to the MatAXPY() operations to perform the updates.
1185: If globalup is true, communication is carried out to reconstruct the updated
1186: matrices in the parent communicator. The user must be warned that if global
1187: matrices are not in sync with subcommunicator matrices, the errors computed
1188: by EPSComputeError() will be wrong even if the computed solution is correct
1189: (the synchronization may be done only once at the end).
1191: Level: advanced
1193: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommMats()
1194: @*/
1195: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1196: {
1197: PetscFunctionBegin;
1207: PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1208: PetscFunctionReturn(PETSC_SUCCESS);
1209: }
1211: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1212: {
1213: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1214: Mat A,B=NULL,Ar=NULL,Br=NULL;
1215: PetscMPIInt rank;
1216: PetscObjectState Astate,Bstate=0;
1217: PetscObjectId Aid,Bid=0;
1218: STType sttype;
1219: PetscInt nmat;
1220: const char *prefix;
1221: MPI_Comm child;
1223: PetscFunctionBegin;
1224: PetscCall(EPSGetOperators(eps,&A,&B));
1225: if (ctx->npart==1) {
1226: if (!ctx->eps) PetscCall(EPSCreate(((PetscObject)eps)->comm,&ctx->eps));
1227: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1228: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1229: PetscCall(EPSSetOperators(ctx->eps,A,B));
1230: } else {
1231: PetscCall(MatGetState(A,&Astate));
1232: PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1233: if (B) {
1234: PetscCall(MatGetState(B,&Bstate));
1235: PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1236: }
1237: if (!ctx->subc) {
1238: /* Create context for subcommunicators */
1239: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1240: PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1241: PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1242: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1244: /* Duplicate matrices */
1245: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1246: ctx->Astate = Astate;
1247: ctx->Aid = Aid;
1248: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1249: if (B) {
1250: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1251: ctx->Bstate = Bstate;
1252: ctx->Bid = Bid;
1253: PetscCall(MatPropagateSymmetryOptions(B,Br));
1254: }
1255: } else {
1256: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1257: if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1258: PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1259: if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
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: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1271: PetscCall(MatDestroy(&Ar));
1272: PetscCall(MatDestroy(&Br));
1273: }
1274: }
1276: /* Create auxiliary EPS */
1277: if (!ctx->eps) {
1278: PetscCall(EPSCreate(child,&ctx->eps));
1279: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1280: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1281: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1282: PetscCall(MatDestroy(&Ar));
1283: PetscCall(MatDestroy(&Br));
1284: }
1285: /* Create subcommunicator grouping processes with same rank */
1286: if (!ctx->commset) {
1287: PetscCallMPI(MPI_Comm_rank(child,&rank));
1288: PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1289: ctx->commset = PETSC_TRUE;
1290: }
1291: }
1292: PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1293: PetscCall(STGetType(eps->st,&sttype));
1294: PetscCall(STSetType(ctx->eps->st,sttype));
1296: ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1297: ctx_local->npart = ctx->npart;
1298: ctx_local->global = PETSC_FALSE;
1299: ctx_local->eps = eps;
1300: ctx_local->subc = ctx->subc;
1301: ctx_local->commrank = ctx->commrank;
1302: *childeps = ctx->eps;
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1307: {
1308: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1309: ST st;
1310: PetscBool isfilt;
1312: PetscFunctionBegin;
1313: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1314: PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1315: PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1316: PetscCall(EPSGetST(ctx->eps,&st));
1317: PetscCall(STGetOperator(st,NULL));
1318: PetscCall(STGetKSP(st,ksp));
1319: PetscFunctionReturn(PETSC_SUCCESS);
1320: }
1322: /*@
1323: EPSKrylovSchurGetKSP - Retrieve the linear solver object associated with the
1324: internal EPS object in case of doing spectrum slicing for a computational interval.
1326: Collective
1328: Input Parameter:
1329: . eps - the eigenproblem solver context
1331: Output Parameter:
1332: . ksp - the internal KSP object
1334: Notes:
1335: When invoked to compute all eigenvalues in an interval with spectrum
1336: slicing, EPSKRYLOVSCHUR creates another EPS object internally that is
1337: used to compute eigenvalues by chunks near selected shifts. This function
1338: allows access to the KSP object associated to this internal EPS object.
1340: This function is only available for spectrum slicing runs. In case of
1341: having more than one partition, the returned KSP will be different
1342: in MPI processes belonging to different partitions. Hence, if required,
1343: EPSKrylovSchurSetPartitions() must be called BEFORE this function.
1345: Level: advanced
1347: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions()
1348: @*/
1349: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1350: {
1351: PetscFunctionBegin;
1353: PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: static PetscErrorCode EPSKrylovSchurSetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType bse)
1358: {
1359: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1361: PetscFunctionBegin;
1362: switch (bse) {
1363: case EPS_KRYLOVSCHUR_BSE_SHAO:
1364: case EPS_KRYLOVSCHUR_BSE_GRUNING:
1365: case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
1366: if (ctx->bse != bse) {
1367: ctx->bse = bse;
1368: eps->state = EPS_STATE_INITIAL;
1369: }
1370: break;
1371: default:
1372: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid BSE type");
1373: }
1374: PetscFunctionReturn(PETSC_SUCCESS);
1375: }
1377: /*@
1378: EPSKrylovSchurSetBSEType - Sets the method to be used for BSE structured eigenproblems
1379: in the Krylov-Schur solver.
1381: Logically Collective
1383: Input Parameters:
1384: + eps - the eigenproblem solver context
1385: - bse - the BSE method
1387: Options Database Key:
1388: . -eps_krylovschur_bse_type - Sets the BSE type (either 'shao', 'gruning', or 'projectedbse')
1390: Level: advanced
1392: .seealso: EPSKrylovSchurGetBSEType(), EPSKrylovSchurBSEType, MatCreateBSE()
1393: @*/
1394: PetscErrorCode EPSKrylovSchurSetBSEType(EPS eps,EPSKrylovSchurBSEType bse)
1395: {
1396: PetscFunctionBegin;
1399: PetscTryMethod(eps,"EPSKrylovSchurSetBSEType_C",(EPS,EPSKrylovSchurBSEType),(eps,bse));
1400: PetscFunctionReturn(PETSC_SUCCESS);
1401: }
1403: static PetscErrorCode EPSKrylovSchurGetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType *bse)
1404: {
1405: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1407: PetscFunctionBegin;
1408: *bse = ctx->bse;
1409: PetscFunctionReturn(PETSC_SUCCESS);
1410: }
1412: /*@
1413: EPSKrylovSchurGetBSEType - Gets the method used for BSE structured eigenproblems
1414: in the Krylov-Schur solver.
1416: Not Collective
1418: Input Parameter:
1419: . eps - the eigenproblem solver context
1421: Output Parameter:
1422: . bse - the BSE method
1424: Level: advanced
1426: .seealso: EPSKrylovSchurSetBSEType(), EPSKrylovSchurBSEType, MatCreateBSE()
1427: @*/
1428: PetscErrorCode EPSKrylovSchurGetBSEType(EPS eps,EPSKrylovSchurBSEType *bse)
1429: {
1430: PetscFunctionBegin;
1432: PetscAssertPointer(bse,2);
1433: PetscUseMethod(eps,"EPSKrylovSchurGetBSEType_C",(EPS,EPSKrylovSchurBSEType*),(eps,bse));
1434: PetscFunctionReturn(PETSC_SUCCESS);
1435: }
1437: static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems *PetscOptionsObject)
1438: {
1439: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1440: PetscBool flg,lock,b,f1,f2,f3,isfilt;
1441: PetscReal keep;
1442: PetscInt i,j,k;
1443: KSP ksp;
1444: EPSKrylovSchurBSEType bse;
1446: PetscFunctionBegin;
1447: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");
1449: PetscCall(PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg));
1450: if (flg) PetscCall(EPSKrylovSchurSetRestart(eps,keep));
1452: PetscCall(PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg));
1453: if (flg) PetscCall(EPSKrylovSchurSetLocking(eps,lock));
1455: i = ctx->npart;
1456: PetscCall(PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg));
1457: if (flg) PetscCall(EPSKrylovSchurSetPartitions(eps,i));
1459: b = ctx->detect;
1460: PetscCall(PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg));
1461: if (flg) PetscCall(EPSKrylovSchurSetDetectZeros(eps,b));
1463: i = 1;
1464: j = k = PETSC_DECIDE;
1465: PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1466: PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1467: PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1468: if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));
1470: PetscCall(PetscOptionsEnum("-eps_krylovschur_bse_type","Method for BSE structured eigenproblems","EPSKrylovSchurSetBSEType",EPSKrylovSchurBSETypes,(PetscEnum)ctx->bse,(PetscEnum*)&bse,&flg));
1471: if (flg) PetscCall(EPSKrylovSchurSetBSEType(eps,bse));
1473: PetscOptionsHeadEnd();
1475: /* set options of child KSP in spectrum slicing */
1476: if (eps->which==EPS_ALL) {
1477: if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1478: PetscCall(EPSSetDefaultST(eps));
1479: PetscCall(STSetFromOptions(eps->st)); /* need to advance this to check ST type */
1480: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1481: if (!isfilt) {
1482: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1483: PetscCall(KSPSetFromOptions(ksp));
1484: }
1485: }
1486: PetscFunctionReturn(PETSC_SUCCESS);
1487: }
1489: static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1490: {
1491: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1492: PetscBool isascii,isfilt;
1493: KSP ksp;
1494: PetscViewer sviewer;
1496: PetscFunctionBegin;
1497: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1498: if (isascii) {
1499: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1500: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1501: if (eps->problem_type==EPS_BSE) PetscCall(PetscViewerASCIIPrintf(viewer," BSE method: %s\n",EPSKrylovSchurBSETypes[ctx->bse]));
1502: if (eps->which==EPS_ALL) {
1503: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1504: if (isfilt) PetscCall(PetscViewerASCIIPrintf(viewer," using filtering to extract all eigenvalues in an interval\n"));
1505: else {
1506: PetscCall(PetscViewerASCIIPrintf(viewer," doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd));
1507: if (ctx->npart>1) {
1508: PetscCall(PetscViewerASCIIPrintf(viewer," multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart));
1509: if (ctx->detect) PetscCall(PetscViewerASCIIPrintf(viewer," detecting zeros when factorizing at subinterval boundaries\n"));
1510: }
1511: /* view child KSP */
1512: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1513: PetscCall(PetscViewerASCIIPushTab(viewer));
1514: if (ctx->npart>1 && ctx->subc) {
1515: PetscCall(PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer));
1516: if (!ctx->subc->color) PetscCall(KSPView(ksp,sviewer));
1517: PetscCall(PetscViewerFlush(sviewer));
1518: PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer));
1519: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1520: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1521: } else PetscCall(KSPView(ksp,viewer));
1522: PetscCall(PetscViewerASCIIPopTab(viewer));
1523: }
1524: }
1525: }
1526: PetscFunctionReturn(PETSC_SUCCESS);
1527: }
1529: static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1530: {
1531: PetscBool isfilt;
1533: PetscFunctionBegin;
1534: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1535: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1536: PetscCall(PetscFree(eps->data));
1537: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1538: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1539: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1540: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1541: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1542: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1543: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1544: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1545: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1546: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1547: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1548: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1549: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1550: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1551: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1552: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1553: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1554: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1555: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",NULL));
1556: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",NULL));
1557: PetscFunctionReturn(PETSC_SUCCESS);
1558: }
1560: static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1561: {
1562: PetscBool isfilt;
1564: PetscFunctionBegin;
1565: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1566: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1567: PetscFunctionReturn(PETSC_SUCCESS);
1568: }
1570: static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1571: {
1572: PetscFunctionBegin;
1573: if (eps->which==EPS_ALL) {
1574: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1575: }
1576: PetscFunctionReturn(PETSC_SUCCESS);
1577: }
1579: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1580: {
1581: EPS_KRYLOVSCHUR *ctx;
1583: PetscFunctionBegin;
1584: PetscCall(PetscNew(&ctx));
1585: eps->data = (void*)ctx;
1586: ctx->lock = PETSC_TRUE;
1587: ctx->nev = 1;
1588: ctx->ncv = PETSC_DETERMINE;
1589: ctx->mpd = PETSC_DETERMINE;
1590: ctx->npart = 1;
1591: ctx->detect = PETSC_FALSE;
1592: ctx->global = PETSC_TRUE;
1594: eps->useds = PETSC_TRUE;
1596: /* solve and computevectors determined at setup */
1597: eps->ops->setup = EPSSetUp_KrylovSchur;
1598: eps->ops->setupsort = EPSSetUpSort_KrylovSchur;
1599: eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1600: eps->ops->destroy = EPSDestroy_KrylovSchur;
1601: eps->ops->reset = EPSReset_KrylovSchur;
1602: eps->ops->view = EPSView_KrylovSchur;
1603: eps->ops->backtransform = EPSBackTransform_Default;
1604: eps->ops->setdefaultst = EPSSetDefaultST_KrylovSchur;
1606: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1607: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1608: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1609: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1610: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1611: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1612: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1613: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1614: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1615: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1616: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1617: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1618: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1619: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1620: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1621: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1622: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1623: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1624: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",EPSKrylovSchurSetBSEType_KrylovSchur));
1625: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",EPSKrylovSchurGetBSEType_KrylovSchur));
1626: PetscFunctionReturn(PETSC_SUCCESS);
1627: }