Actual source code: krylovschur.c
slepc-3.20.2 2024-03-15
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_DEFAULT) 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_DEFAULT && 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_DEFAULT) 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 {
108: PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
109: 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");
110: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
111: if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
112: }
113: PetscCheck(ctx->lock || eps->mpd>=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
115: EPSCheckDefiniteCondition(eps,eps->arbitrary," with arbitrary selection of eigenpairs");
117: PetscCheck(eps->extraction==EPS_RITZ || eps->extraction==EPS_HARMONIC,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
119: if (!ctx->keep) ctx->keep = 0.5;
121: PetscCall(EPSAllocateSolution(eps,1));
122: PetscCall(EPS_SetInnerProduct(eps));
123: if (eps->arbitrary) PetscCall(EPSSetWorkVecs(eps,2));
124: else if (eps->ishermitian && !eps->ispositive) PetscCall(EPSSetWorkVecs(eps,1));
126: /* dispatch solve method */
127: if (eps->ishermitian) {
128: if (eps->which==EPS_ALL) {
129: EPSCheckDefiniteCondition(eps,eps->which==EPS_ALL," with spectrum slicing");
130: variant = isfilt? EPS_KS_FILTER: EPS_KS_SLICE;
131: } else if (eps->isgeneralized && !eps->ispositive) {
132: variant = EPS_KS_INDEF;
133: } else {
134: switch (eps->extraction) {
135: case EPS_RITZ: variant = EPS_KS_SYMM; break;
136: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
137: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
138: }
139: }
140: } else if (eps->twosided) {
141: variant = EPS_KS_TWOSIDED;
142: } else {
143: switch (eps->extraction) {
144: case EPS_RITZ: variant = EPS_KS_DEFAULT; break;
145: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
146: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
147: }
148: }
149: switch (variant) {
150: case EPS_KS_DEFAULT:
151: eps->ops->solve = EPSSolve_KrylovSchur_Default;
152: eps->ops->computevectors = EPSComputeVectors_Schur;
153: PetscCall(DSSetType(eps->ds,DSNHEP));
154: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
155: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
156: break;
157: case EPS_KS_SYMM:
158: case EPS_KS_FILTER:
159: eps->ops->solve = EPSSolve_KrylovSchur_Default;
160: eps->ops->computevectors = EPSComputeVectors_Hermitian;
161: PetscCall(DSSetType(eps->ds,DSHEP));
162: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
163: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
164: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
165: break;
166: case EPS_KS_SLICE:
167: eps->ops->solve = EPSSolve_KrylovSchur_Slice;
168: eps->ops->computevectors = EPSComputeVectors_Slice;
169: break;
170: case EPS_KS_INDEF:
171: eps->ops->solve = EPSSolve_KrylovSchur_Indefinite;
172: eps->ops->computevectors = EPSComputeVectors_Indefinite;
173: PetscCall(DSSetType(eps->ds,DSGHIEP));
174: PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
175: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
176: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
177: /* force reorthogonalization for pseudo-Lanczos */
178: PetscCall(BVGetOrthogonalization(eps->V,&otype,NULL,&eta,&obtype));
179: PetscCall(BVSetOrthogonalization(eps->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
180: break;
181: case EPS_KS_TWOSIDED:
182: eps->ops->solve = EPSSolve_KrylovSchur_TwoSided;
183: eps->ops->computevectors = EPSComputeVectors_Schur;
184: PetscCall(DSSetType(eps->ds,DSNHEPTS));
185: PetscCall(DSAllocate(eps->ds,eps->ncv+1));
186: PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
187: break;
188: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
189: }
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: static PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
194: {
195: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
196: SlepcSC sc;
197: PetscBool isfilt;
199: PetscFunctionBegin;
200: PetscCall(EPSSetUpSort_Default(eps));
201: if (eps->which==EPS_ALL) {
202: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
203: if (isfilt) {
204: PetscCall(DSGetSlepcSC(eps->ds,&sc));
205: sc->rg = NULL;
206: sc->comparison = SlepcCompareLargestReal;
207: sc->comparisonctx = NULL;
208: sc->map = NULL;
209: sc->mapobj = NULL;
210: } else {
211: if (!ctx->global && ctx->sr->numEigs>0) {
212: PetscCall(DSGetSlepcSC(eps->ds,&sc));
213: sc->rg = NULL;
214: sc->comparison = SlepcCompareLargestMagnitude;
215: sc->comparisonctx = NULL;
216: sc->map = NULL;
217: sc->mapobj = NULL;
218: }
219: }
220: }
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: PetscErrorCode EPSSolve_KrylovSchur_Default(EPS eps)
225: {
226: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
227: PetscInt j,*pj,k,l,nv,ld,nconv;
228: Mat U,Op,H,T;
229: PetscScalar *g;
230: PetscReal beta,gamma=1.0;
231: PetscBool breakdown,harmonic,hermitian;
233: PetscFunctionBegin;
234: PetscCall(DSGetLeadingDimension(eps->ds,&ld));
235: harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
236: hermitian = (eps->ishermitian && !harmonic)?PETSC_TRUE:PETSC_FALSE;
237: if (harmonic) PetscCall(PetscMalloc1(ld,&g));
238: if (eps->arbitrary) pj = &j;
239: else pj = NULL;
241: /* Get the starting Arnoldi vector */
242: PetscCall(EPSGetStartVector(eps,0,NULL));
243: l = 0;
245: /* Restart loop */
246: while (eps->reason == EPS_CONVERGED_ITERATING) {
247: eps->its++;
249: /* Compute an nv-step Arnoldi factorization */
250: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
251: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
252: PetscCall(STGetOperator(eps->st,&Op));
253: if (hermitian) {
254: PetscCall(DSGetMat(eps->ds,DS_MAT_T,&T));
255: PetscCall(BVMatLanczos(eps->V,Op,T,eps->nconv+l,&nv,&beta,&breakdown));
256: PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&T));
257: } else {
258: PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
259: PetscCall(BVMatArnoldi(eps->V,Op,H,eps->nconv+l,&nv,&beta,&breakdown));
260: PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
261: }
262: PetscCall(STRestoreOperator(eps->st,&Op));
263: PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
264: PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
265: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
267: /* Compute translation of Krylov decomposition if harmonic extraction used */
268: if (PetscUnlikely(harmonic)) PetscCall(DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,g,&gamma));
270: /* Solve projected problem */
271: PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
272: if (PetscUnlikely(eps->arbitrary)) {
273: PetscCall(EPSGetArbitraryValues(eps,eps->rr,eps->ri));
274: j=1;
275: }
276: PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj));
277: PetscCall(DSUpdateExtraRow(eps->ds));
278: PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));
280: /* Check convergence */
281: PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
282: PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
283: nconv = k;
285: /* Update l */
286: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
287: else {
288: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
289: if (!hermitian) PetscCall(DSGetTruncateSize(eps->ds,k,nv,&l));
290: }
291: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
292: if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
294: if (eps->reason == EPS_CONVERGED_ITERATING) {
295: if (PetscUnlikely(breakdown || k==nv)) {
296: /* Start a new Arnoldi factorization */
297: PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
298: if (k<eps->nev) {
299: PetscCall(EPSGetStartVector(eps,k,&breakdown));
300: if (breakdown) {
301: eps->reason = EPS_DIVERGED_BREAKDOWN;
302: PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
303: }
304: }
305: } else {
306: /* Undo translation of Krylov decomposition */
307: if (PetscUnlikely(harmonic)) {
308: PetscCall(DSSetDimensions(eps->ds,nv,k,l));
309: PetscCall(DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma));
310: /* gamma u^ = u - U*g~ */
311: PetscCall(BVSetActiveColumns(eps->V,0,nv));
312: PetscCall(BVMultColumn(eps->V,-1.0,1.0,nv,g));
313: PetscCall(BVScaleColumn(eps->V,nv,1.0/gamma));
314: PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
315: PetscCall(DSSetDimensions(eps->ds,nv,k,nv));
316: }
317: /* Prepare the Rayleigh quotient for restart */
318: PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
319: }
320: }
321: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
322: PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
323: PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
324: PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));
326: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(eps->V,nv,k+l));
327: eps->nconv = k;
328: PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
329: }
331: if (harmonic) PetscCall(PetscFree(g));
332: PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
337: {
338: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
340: PetscFunctionBegin;
341: if (keep==(PetscReal)PETSC_DEFAULT) ctx->keep = 0.5;
342: else {
343: 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);
344: ctx->keep = keep;
345: }
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*@
350: EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
351: method, in particular the proportion of basis vectors that must be kept
352: after restart.
354: Logically Collective
356: Input Parameters:
357: + eps - the eigenproblem solver context
358: - keep - the number of vectors to be kept at restart
360: Options Database Key:
361: . -eps_krylovschur_restart - Sets the restart parameter
363: Notes:
364: Allowed values are in the range [0.1,0.9]. The default is 0.5.
366: Level: advanced
368: .seealso: EPSKrylovSchurGetRestart()
369: @*/
370: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
371: {
372: PetscFunctionBegin;
375: PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
376: PetscFunctionReturn(PETSC_SUCCESS);
377: }
379: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
380: {
381: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
383: PetscFunctionBegin;
384: *keep = ctx->keep;
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /*@
389: EPSKrylovSchurGetRestart - Gets the restart parameter used in the
390: Krylov-Schur method.
392: Not Collective
394: Input Parameter:
395: . eps - the eigenproblem solver context
397: Output Parameter:
398: . keep - the restart parameter
400: Level: advanced
402: .seealso: EPSKrylovSchurSetRestart()
403: @*/
404: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
405: {
406: PetscFunctionBegin;
408: PetscAssertPointer(keep,2);
409: PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
414: {
415: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
417: PetscFunctionBegin;
418: ctx->lock = lock;
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: /*@
423: EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
424: the Krylov-Schur method.
426: Logically Collective
428: Input Parameters:
429: + eps - the eigenproblem solver context
430: - lock - true if the locking variant must be selected
432: Options Database Key:
433: . -eps_krylovschur_locking - Sets the locking flag
435: Notes:
436: The default is to lock converged eigenpairs when the method restarts.
437: This behaviour can be changed so that all directions are kept in the
438: working subspace even if already converged to working accuracy (the
439: non-locking variant).
441: Level: advanced
443: .seealso: EPSKrylovSchurGetLocking()
444: @*/
445: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
446: {
447: PetscFunctionBegin;
450: PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
454: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
455: {
456: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
458: PetscFunctionBegin;
459: *lock = ctx->lock;
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: /*@
464: EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
465: method.
467: Not Collective
469: Input Parameter:
470: . eps - the eigenproblem solver context
472: Output Parameter:
473: . lock - the locking flag
475: Level: advanced
477: .seealso: EPSKrylovSchurSetLocking()
478: @*/
479: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
480: {
481: PetscFunctionBegin;
483: PetscAssertPointer(lock,2);
484: PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
489: {
490: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
491: PetscMPIInt size;
492: PetscInt newnpart;
494: PetscFunctionBegin;
495: if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
496: newnpart = 1;
497: } else {
498: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
499: PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
500: newnpart = npart;
501: }
502: if (ctx->npart!=newnpart) {
503: if (ctx->npart>1) {
504: PetscCall(PetscSubcommDestroy(&ctx->subc));
505: if (ctx->commset) {
506: PetscCallMPI(MPI_Comm_free(&ctx->commrank));
507: ctx->commset = PETSC_FALSE;
508: }
509: }
510: PetscCall(EPSDestroy(&ctx->eps));
511: ctx->npart = newnpart;
512: eps->state = EPS_STATE_INITIAL;
513: }
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: /*@
518: EPSKrylovSchurSetPartitions - Sets the number of partitions for the
519: case of doing spectrum slicing for a computational interval with the
520: communicator split in several sub-communicators.
522: Logically Collective
524: Input Parameters:
525: + eps - the eigenproblem solver context
526: - npart - number of partitions
528: Options Database Key:
529: . -eps_krylovschur_partitions <npart> - Sets the number of partitions
531: Notes:
532: By default, npart=1 so all processes in the communicator participate in
533: the processing of the whole interval. If npart>1 then the interval is
534: divided into npart subintervals, each of them being processed by a
535: subset of processes.
537: The interval is split proportionally unless the separation points are
538: specified with EPSKrylovSchurSetSubintervals().
540: Level: advanced
542: .seealso: EPSKrylovSchurSetSubintervals(), EPSSetInterval()
543: @*/
544: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
545: {
546: PetscFunctionBegin;
549: PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
554: {
555: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
557: PetscFunctionBegin;
558: *npart = ctx->npart;
559: PetscFunctionReturn(PETSC_SUCCESS);
560: }
562: /*@
563: EPSKrylovSchurGetPartitions - Gets the number of partitions of the
564: communicator in case of spectrum slicing.
566: Not Collective
568: Input Parameter:
569: . eps - the eigenproblem solver context
571: Output Parameter:
572: . npart - number of partitions
574: Level: advanced
576: .seealso: EPSKrylovSchurSetPartitions()
577: @*/
578: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
579: {
580: PetscFunctionBegin;
582: PetscAssertPointer(npart,2);
583: PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
587: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
588: {
589: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
591: PetscFunctionBegin;
592: ctx->detect = detect;
593: eps->state = EPS_STATE_INITIAL;
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: /*@
598: EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
599: zeros during the factorizations throughout the spectrum slicing computation.
601: Logically Collective
603: Input Parameters:
604: + eps - the eigenproblem solver context
605: - detect - check for zeros
607: Options Database Key:
608: . -eps_krylovschur_detect_zeros - Check for zeros; this takes an optional
609: bool value (0/1/no/yes/true/false)
611: Notes:
612: A zero in the factorization indicates that a shift coincides with an eigenvalue.
614: This flag is turned off by default, and may be necessary in some cases,
615: especially when several partitions are being used. This feature currently
616: requires an external package for factorizations with support for zero
617: detection, e.g. MUMPS.
619: Level: advanced
621: .seealso: EPSKrylovSchurSetPartitions(), EPSSetInterval()
622: @*/
623: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
624: {
625: PetscFunctionBegin;
628: PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
633: {
634: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
636: PetscFunctionBegin;
637: *detect = ctx->detect;
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: /*@
642: EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
643: in spectrum slicing.
645: Not Collective
647: Input Parameter:
648: . eps - the eigenproblem solver context
650: Output Parameter:
651: . detect - whether zeros detection is enforced during factorizations
653: Level: advanced
655: .seealso: EPSKrylovSchurSetDetectZeros()
656: @*/
657: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
658: {
659: PetscFunctionBegin;
661: PetscAssertPointer(detect,2);
662: PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
667: {
668: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
670: PetscFunctionBegin;
671: PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
672: ctx->nev = nev;
673: if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
674: ctx->ncv = PETSC_DEFAULT;
675: } else {
676: PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
677: ctx->ncv = ncv;
678: }
679: if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
680: ctx->mpd = PETSC_DEFAULT;
681: } else {
682: PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
683: ctx->mpd = mpd;
684: }
685: eps->state = EPS_STATE_INITIAL;
686: PetscFunctionReturn(PETSC_SUCCESS);
687: }
689: /*@
690: EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
691: step in case of doing spectrum slicing for a computational interval.
692: The meaning of the parameters is the same as in EPSSetDimensions().
694: Logically Collective
696: Input Parameters:
697: + eps - the eigenproblem solver context
698: . nev - number of eigenvalues to compute
699: . ncv - the maximum dimension of the subspace to be used by the subsolve
700: - mpd - the maximum dimension allowed for the projected problem
702: Options Database Key:
703: + -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
704: . -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
705: - -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension
707: Level: advanced
709: .seealso: EPSKrylovSchurGetDimensions(), EPSSetDimensions(), EPSSetInterval()
710: @*/
711: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
712: {
713: PetscFunctionBegin;
718: PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
719: PetscFunctionReturn(PETSC_SUCCESS);
720: }
722: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
723: {
724: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
726: PetscFunctionBegin;
727: if (nev) *nev = ctx->nev;
728: if (ncv) *ncv = ctx->ncv;
729: if (mpd) *mpd = ctx->mpd;
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: /*@
734: EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
735: step in case of doing spectrum slicing for a computational interval.
737: Not Collective
739: Input Parameter:
740: . eps - the eigenproblem solver context
742: Output Parameters:
743: + nev - number of eigenvalues to compute
744: . ncv - the maximum dimension of the subspace to be used by the subsolve
745: - mpd - the maximum dimension allowed for the projected problem
747: Level: advanced
749: .seealso: EPSKrylovSchurSetDimensions()
750: @*/
751: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
752: {
753: PetscFunctionBegin;
755: PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
756: PetscFunctionReturn(PETSC_SUCCESS);
757: }
759: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
760: {
761: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
762: PetscInt i;
764: PetscFunctionBegin;
765: 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()");
766: 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");
767: if (ctx->subintervals) PetscCall(PetscFree(ctx->subintervals));
768: PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
769: for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
770: ctx->subintset = PETSC_TRUE;
771: eps->state = EPS_STATE_INITIAL;
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: /*@C
776: EPSKrylovSchurSetSubintervals - Sets the points that delimit the
777: subintervals to be used in spectrum slicing with several partitions.
779: Logically Collective
781: Input Parameters:
782: + eps - the eigenproblem solver context
783: - subint - array of real values specifying subintervals
785: Notes:
786: This function must be called after EPSKrylovSchurSetPartitions(). For npart
787: partitions, the argument subint must contain npart+1 real values sorted in
788: ascending order, subint_0, subint_1, ..., subint_npart, where the first
789: and last values must coincide with the interval endpoints set with
790: EPSSetInterval().
792: The subintervals are then defined by two consecutive points [subint_0,subint_1],
793: [subint_1,subint_2], and so on.
795: Level: advanced
797: .seealso: EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubintervals(), EPSSetInterval()
798: @*/
799: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal *subint)
800: {
801: PetscFunctionBegin;
803: PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
804: PetscFunctionReturn(PETSC_SUCCESS);
805: }
807: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal **subint)
808: {
809: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
810: PetscInt i;
812: PetscFunctionBegin;
813: if (!ctx->subintset) {
814: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
815: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
816: }
817: PetscCall(PetscMalloc1(ctx->npart+1,subint));
818: for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: /*@C
823: EPSKrylovSchurGetSubintervals - Returns the points that delimit the
824: subintervals used in spectrum slicing with several partitions.
826: Not Collective
828: Input Parameter:
829: . eps - the eigenproblem solver context
831: Output Parameter:
832: . subint - array of real values specifying subintervals
834: Notes:
835: If the user passed values with EPSKrylovSchurSetSubintervals(), then the
836: same values are returned. Otherwise, the values computed internally are
837: obtained.
839: This function is only available for spectrum slicing runs.
841: The returned array has length npart+1 (see EPSKrylovSchurGetPartitions())
842: and should be freed by the user.
844: Fortran Notes:
845: The calling sequence from Fortran is
846: .vb
847: EPSKrylovSchurGetSubintervals(eps,subint,ierr)
848: double precision subint(npart+1) output
849: .ve
851: Level: advanced
853: .seealso: EPSKrylovSchurSetSubintervals(), EPSKrylovSchurGetPartitions(), EPSSetInterval()
854: @*/
855: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal **subint)
856: {
857: PetscFunctionBegin;
859: PetscAssertPointer(subint,2);
860: PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
865: {
866: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
867: PetscInt i,numsh;
868: EPS_SR sr = ctx->sr;
870: PetscFunctionBegin;
871: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
872: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
873: switch (eps->state) {
874: case EPS_STATE_INITIAL:
875: break;
876: case EPS_STATE_SETUP:
877: numsh = ctx->npart+1;
878: if (n) *n = numsh;
879: if (shifts) {
880: PetscCall(PetscMalloc1(numsh,shifts));
881: (*shifts)[0] = eps->inta;
882: if (ctx->npart==1) (*shifts)[1] = eps->intb;
883: else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
884: }
885: if (inertias) {
886: PetscCall(PetscMalloc1(numsh,inertias));
887: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
888: if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
889: else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
890: }
891: break;
892: case EPS_STATE_SOLVED:
893: case EPS_STATE_EIGENVECTORS:
894: numsh = ctx->nshifts;
895: if (n) *n = numsh;
896: if (shifts) {
897: PetscCall(PetscMalloc1(numsh,shifts));
898: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
899: }
900: if (inertias) {
901: PetscCall(PetscMalloc1(numsh,inertias));
902: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
903: }
904: break;
905: }
906: PetscFunctionReturn(PETSC_SUCCESS);
907: }
909: /*@C
910: EPSKrylovSchurGetInertias - Gets the values of the shifts and their
911: corresponding inertias in case of doing spectrum slicing for a
912: computational interval.
914: Not Collective
916: Input Parameter:
917: . eps - the eigenproblem solver context
919: Output Parameters:
920: + n - number of shifts, including the endpoints of the interval
921: . shifts - the values of the shifts used internally in the solver
922: - inertias - the values of the inertia in each shift
924: Notes:
925: If called after EPSSolve(), all shifts used internally by the solver are
926: returned (including both endpoints and any intermediate ones). If called
927: before EPSSolve() and after EPSSetUp() then only the information of the
928: endpoints of subintervals is available.
930: This function is only available for spectrum slicing runs.
932: The returned arrays should be freed by the user. Can pass NULL in any of
933: the two arrays if not required.
935: Fortran Notes:
936: The calling sequence from Fortran is
937: .vb
938: EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
939: integer n
940: double precision shifts(*)
941: integer inertias(*)
942: .ve
943: The arrays should be at least of length n. The value of n can be determined
944: by an initial call
945: .vb
946: EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
947: .ve
949: Level: advanced
951: .seealso: EPSSetInterval(), EPSKrylovSchurSetSubintervals()
952: @*/
953: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
954: {
955: PetscFunctionBegin;
957: PetscAssertPointer(n,2);
958: PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }
962: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
963: {
964: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
965: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
967: PetscFunctionBegin;
968: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
969: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
970: if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
971: if (n) *n = sr->numEigs;
972: if (v) PetscCall(BVCreateVec(sr->V,v));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: /*@C
977: EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
978: doing spectrum slicing for a computational interval with multiple
979: communicators.
981: Collective on the subcommunicator (if v is given)
983: Input Parameter:
984: . eps - the eigenproblem solver context
986: Output Parameters:
987: + k - index of the subinterval for the calling process
988: . n - number of eigenvalues found in the k-th subinterval
989: - v - a vector owned by processes in the subcommunicator with dimensions
990: compatible for locally computed eigenvectors (or NULL)
992: Notes:
993: This function is only available for spectrum slicing runs.
995: The returned Vec should be destroyed by the user.
997: Level: advanced
999: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommPairs()
1000: @*/
1001: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1002: {
1003: PetscFunctionBegin;
1005: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1006: PetscFunctionReturn(PETSC_SUCCESS);
1007: }
1009: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1010: {
1011: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1012: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1014: PetscFunctionBegin;
1015: EPSCheckSolved(eps,1);
1016: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1017: PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1018: if (eig) *eig = sr->eigr[sr->perm[i]];
1019: if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1020: PetscFunctionReturn(PETSC_SUCCESS);
1021: }
1023: /*@
1024: EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
1025: internally in the subcommunicator to which the calling process belongs.
1027: Collective on the subcommunicator (if v is given)
1029: Input Parameters:
1030: + eps - the eigenproblem solver context
1031: - i - index of the solution
1033: Output Parameters:
1034: + eig - the eigenvalue
1035: - v - the eigenvector
1037: Notes:
1038: It is allowed to pass NULL for v if the eigenvector is not required.
1039: Otherwise, the caller must provide a valid Vec objects, i.e.,
1040: it must be created by the calling program with EPSKrylovSchurGetSubcommInfo().
1042: The index i should be a value between 0 and n-1, where n is the number of
1043: vectors in the local subinterval, see EPSKrylovSchurGetSubcommInfo().
1045: Level: advanced
1047: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo(), EPSKrylovSchurGetSubcommMats()
1048: @*/
1049: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1050: {
1051: PetscFunctionBegin;
1054: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1059: {
1060: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1062: PetscFunctionBegin;
1063: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1064: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1065: PetscCall(EPSGetOperators(ctx->eps,A,B));
1066: PetscFunctionReturn(PETSC_SUCCESS);
1067: }
1069: /*@C
1070: EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1071: internally in the subcommunicator to which the calling process belongs.
1073: Collective on the subcommunicator
1075: Input Parameter:
1076: . eps - the eigenproblem solver context
1078: Output Parameters:
1079: + A - the matrix associated with the eigensystem
1080: - B - the second matrix in the case of generalized eigenproblems
1082: Notes:
1083: This is the analog of EPSGetOperators(), but returns the matrices distributed
1084: differently (in the subcommunicator rather than in the parent communicator).
1086: These matrices should not be modified by the user.
1088: Level: advanced
1090: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo()
1091: @*/
1092: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1093: {
1094: PetscFunctionBegin;
1096: PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1097: PetscFunctionReturn(PETSC_SUCCESS);
1098: }
1100: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1101: {
1102: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1103: Mat A,B=NULL,Ag,Bg=NULL;
1104: PetscBool reuse=PETSC_TRUE;
1106: PetscFunctionBegin;
1107: PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1108: PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1109: PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1110: PetscCall(EPSGetOperators(ctx->eps,&A,&B));
1112: PetscCall(MatScale(A,a));
1113: if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1114: if (B) PetscCall(MatScale(B,b));
1115: if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1116: PetscCall(EPSSetOperators(ctx->eps,A,B));
1118: /* Update stored matrix state */
1119: subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1120: PetscCall(PetscObjectStateGet((PetscObject)A,&subctx->Astate));
1121: if (B) PetscCall(PetscObjectStateGet((PetscObject)B,&subctx->Bstate));
1123: /* Update matrices in the parent communicator if requested by user */
1124: if (globalup) {
1125: if (ctx->npart>1) {
1126: if (!ctx->isrow) {
1127: PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1128: reuse = PETSC_FALSE;
1129: }
1130: if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1131: if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1132: PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1133: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1134: if (B) {
1135: if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1136: PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1137: PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1138: }
1139: }
1140: PetscCall(PetscObjectStateGet((PetscObject)Ag,&ctx->Astate));
1141: if (Bg) PetscCall(PetscObjectStateGet((PetscObject)Bg,&ctx->Bstate));
1142: }
1143: PetscCall(EPSSetOperators(eps,Ag,Bg));
1144: PetscFunctionReturn(PETSC_SUCCESS);
1145: }
1147: /*@
1148: EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1149: internally in the subcommunicator to which the calling process belongs.
1151: Collective
1153: Input Parameters:
1154: + eps - the eigenproblem solver context
1155: . s - scalar that multiplies the existing A matrix
1156: . a - scalar used in the axpy operation on A
1157: . Au - matrix used in the axpy operation on A
1158: . t - scalar that multiplies the existing B matrix
1159: . b - scalar used in the axpy operation on B
1160: . Bu - matrix used in the axpy operation on B
1161: . str - structure flag
1162: - globalup - flag indicating if global matrices must be updated
1164: Notes:
1165: This function modifies the eigenproblem matrices at the subcommunicator level,
1166: and optionally updates the global matrices in the parent communicator. The updates
1167: are expressed as A <-- s*A + a*Au, B <-- t*B + b*Bu.
1169: It is possible to update one of the matrices, or both.
1171: The matrices Au and Bu must be equal in all subcommunicators.
1173: The str flag is passed to the MatAXPY() operations to perform the updates.
1175: If globalup is true, communication is carried out to reconstruct the updated
1176: matrices in the parent communicator. The user must be warned that if global
1177: matrices are not in sync with subcommunicator matrices, the errors computed
1178: by EPSComputeError() will be wrong even if the computed solution is correct
1179: (the synchronization may be done only once at the end).
1181: Level: advanced
1183: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommMats()
1184: @*/
1185: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1186: {
1187: PetscFunctionBegin;
1197: PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1198: PetscFunctionReturn(PETSC_SUCCESS);
1199: }
1201: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1202: {
1203: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1204: Mat A,B=NULL,Ar=NULL,Br=NULL;
1205: PetscMPIInt rank;
1206: PetscObjectState Astate,Bstate=0;
1207: PetscObjectId Aid,Bid=0;
1208: STType sttype;
1209: PetscInt nmat;
1210: const char *prefix;
1211: MPI_Comm child;
1213: PetscFunctionBegin;
1214: PetscCall(EPSGetOperators(eps,&A,&B));
1215: if (ctx->npart==1) {
1216: if (!ctx->eps) PetscCall(EPSCreate(((PetscObject)eps)->comm,&ctx->eps));
1217: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1218: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1219: PetscCall(EPSSetOperators(ctx->eps,A,B));
1220: } else {
1221: PetscCall(PetscObjectStateGet((PetscObject)A,&Astate));
1222: PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1223: if (B) {
1224: PetscCall(PetscObjectStateGet((PetscObject)B,&Bstate));
1225: PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1226: }
1227: if (!ctx->subc) {
1228: /* Create context for subcommunicators */
1229: PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1230: PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1231: PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1232: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1234: /* Duplicate matrices */
1235: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1236: ctx->Astate = Astate;
1237: ctx->Aid = Aid;
1238: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1239: if (B) {
1240: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1241: ctx->Bstate = Bstate;
1242: ctx->Bid = Bid;
1243: PetscCall(MatPropagateSymmetryOptions(B,Br));
1244: }
1245: } else {
1246: PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1247: if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1248: PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1249: if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
1250: PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1251: ctx->Astate = Astate;
1252: ctx->Aid = Aid;
1253: PetscCall(MatPropagateSymmetryOptions(A,Ar));
1254: if (B) {
1255: PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1256: ctx->Bstate = Bstate;
1257: ctx->Bid = Bid;
1258: PetscCall(MatPropagateSymmetryOptions(B,Br));
1259: }
1260: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1261: PetscCall(MatDestroy(&Ar));
1262: PetscCall(MatDestroy(&Br));
1263: }
1264: }
1266: /* Create auxiliary EPS */
1267: if (!ctx->eps) {
1268: PetscCall(EPSCreate(child,&ctx->eps));
1269: PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1270: PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1271: PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1272: PetscCall(MatDestroy(&Ar));
1273: PetscCall(MatDestroy(&Br));
1274: }
1275: /* Create subcommunicator grouping processes with same rank */
1276: if (!ctx->commset) {
1277: PetscCallMPI(MPI_Comm_rank(child,&rank));
1278: PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1279: ctx->commset = PETSC_TRUE;
1280: }
1281: }
1282: PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1283: PetscCall(STGetType(eps->st,&sttype));
1284: PetscCall(STSetType(ctx->eps->st,sttype));
1286: ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1287: ctx_local->npart = ctx->npart;
1288: ctx_local->global = PETSC_FALSE;
1289: ctx_local->eps = eps;
1290: ctx_local->subc = ctx->subc;
1291: ctx_local->commrank = ctx->commrank;
1292: *childeps = ctx->eps;
1293: PetscFunctionReturn(PETSC_SUCCESS);
1294: }
1296: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1297: {
1298: EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1299: ST st;
1300: PetscBool isfilt;
1302: PetscFunctionBegin;
1303: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1304: PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1305: PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1306: PetscCall(EPSGetST(ctx->eps,&st));
1307: PetscCall(STGetOperator(st,NULL));
1308: PetscCall(STGetKSP(st,ksp));
1309: PetscFunctionReturn(PETSC_SUCCESS);
1310: }
1312: /*@
1313: EPSKrylovSchurGetKSP - Retrieve the linear solver object associated with the
1314: internal EPS object in case of doing spectrum slicing for a computational interval.
1316: Collective
1318: Input Parameter:
1319: . eps - the eigenproblem solver context
1321: Output Parameter:
1322: . ksp - the internal KSP object
1324: Notes:
1325: When invoked to compute all eigenvalues in an interval with spectrum
1326: slicing, EPSKRYLOVSCHUR creates another EPS object internally that is
1327: used to compute eigenvalues by chunks near selected shifts. This function
1328: allows access to the KSP object associated to this internal EPS object.
1330: This function is only available for spectrum slicing runs. In case of
1331: having more than one partition, the returned KSP will be different
1332: in MPI processes belonging to different partitions. Hence, if required,
1333: EPSKrylovSchurSetPartitions() must be called BEFORE this function.
1335: Level: advanced
1337: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions()
1338: @*/
1339: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1340: {
1341: PetscFunctionBegin;
1343: PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1344: PetscFunctionReturn(PETSC_SUCCESS);
1345: }
1347: static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems *PetscOptionsObject)
1348: {
1349: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1350: PetscBool flg,lock,b,f1,f2,f3,isfilt;
1351: PetscReal keep;
1352: PetscInt i,j,k;
1353: KSP ksp;
1355: PetscFunctionBegin;
1356: PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");
1358: PetscCall(PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg));
1359: if (flg) PetscCall(EPSKrylovSchurSetRestart(eps,keep));
1361: PetscCall(PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg));
1362: if (flg) PetscCall(EPSKrylovSchurSetLocking(eps,lock));
1364: i = ctx->npart;
1365: PetscCall(PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg));
1366: if (flg) PetscCall(EPSKrylovSchurSetPartitions(eps,i));
1368: b = ctx->detect;
1369: PetscCall(PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg));
1370: if (flg) PetscCall(EPSKrylovSchurSetDetectZeros(eps,b));
1372: i = 1;
1373: j = k = PETSC_DECIDE;
1374: PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1375: PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1376: PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1377: if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));
1379: PetscOptionsHeadEnd();
1381: /* set options of child KSP in spectrum slicing */
1382: if (eps->which==EPS_ALL) {
1383: if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1384: PetscCall(EPSSetDefaultST(eps));
1385: PetscCall(STSetFromOptions(eps->st)); /* need to advance this to check ST type */
1386: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1387: if (!isfilt) {
1388: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1389: PetscCall(KSPSetFromOptions(ksp));
1390: }
1391: }
1392: PetscFunctionReturn(PETSC_SUCCESS);
1393: }
1395: static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1396: {
1397: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1398: PetscBool isascii,isfilt;
1399: KSP ksp;
1400: PetscViewer sviewer;
1402: PetscFunctionBegin;
1403: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1404: if (isascii) {
1405: PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1406: PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-"));
1407: if (eps->which==EPS_ALL) {
1408: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1409: if (isfilt) PetscCall(PetscViewerASCIIPrintf(viewer," using filtering to extract all eigenvalues in an interval\n"));
1410: else {
1411: PetscCall(PetscViewerASCIIPrintf(viewer," doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd));
1412: if (ctx->npart>1) {
1413: PetscCall(PetscViewerASCIIPrintf(viewer," multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart));
1414: if (ctx->detect) PetscCall(PetscViewerASCIIPrintf(viewer," detecting zeros when factorizing at subinterval boundaries\n"));
1415: }
1416: /* view child KSP */
1417: PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1418: PetscCall(PetscViewerASCIIPushTab(viewer));
1419: if (ctx->npart>1 && ctx->subc) {
1420: PetscCall(PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer));
1421: if (!ctx->subc->color) PetscCall(KSPView(ksp,sviewer));
1422: PetscCall(PetscViewerFlush(sviewer));
1423: PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer));
1424: PetscCall(PetscViewerFlush(viewer));
1425: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1426: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1427: } else PetscCall(KSPView(ksp,viewer));
1428: PetscCall(PetscViewerASCIIPopTab(viewer));
1429: }
1430: }
1431: }
1432: PetscFunctionReturn(PETSC_SUCCESS);
1433: }
1435: static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1436: {
1437: PetscBool isfilt;
1439: PetscFunctionBegin;
1440: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1441: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1442: PetscCall(PetscFree(eps->data));
1443: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1444: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1445: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1446: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1447: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1448: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1449: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1450: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1451: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1452: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1453: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1454: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1455: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1456: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1457: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1458: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1459: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1460: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1461: PetscFunctionReturn(PETSC_SUCCESS);
1462: }
1464: static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1465: {
1466: PetscBool isfilt;
1468: PetscFunctionBegin;
1469: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1470: if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1471: PetscFunctionReturn(PETSC_SUCCESS);
1472: }
1474: static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1475: {
1476: PetscFunctionBegin;
1477: if (eps->which==EPS_ALL) {
1478: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1479: }
1480: PetscFunctionReturn(PETSC_SUCCESS);
1481: }
1483: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1484: {
1485: EPS_KRYLOVSCHUR *ctx;
1487: PetscFunctionBegin;
1488: PetscCall(PetscNew(&ctx));
1489: eps->data = (void*)ctx;
1490: ctx->lock = PETSC_TRUE;
1491: ctx->nev = 1;
1492: ctx->ncv = PETSC_DEFAULT;
1493: ctx->mpd = PETSC_DEFAULT;
1494: ctx->npart = 1;
1495: ctx->detect = PETSC_FALSE;
1496: ctx->global = PETSC_TRUE;
1498: eps->useds = PETSC_TRUE;
1500: /* solve and computevectors determined at setup */
1501: eps->ops->setup = EPSSetUp_KrylovSchur;
1502: eps->ops->setupsort = EPSSetUpSort_KrylovSchur;
1503: eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1504: eps->ops->destroy = EPSDestroy_KrylovSchur;
1505: eps->ops->reset = EPSReset_KrylovSchur;
1506: eps->ops->view = EPSView_KrylovSchur;
1507: eps->ops->backtransform = EPSBackTransform_Default;
1508: eps->ops->setdefaultst = EPSSetDefaultST_KrylovSchur;
1510: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1511: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1512: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1513: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1514: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1515: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1516: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1517: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1518: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1519: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1520: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1521: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1522: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1523: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1524: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1525: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1526: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1527: PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1528: PetscFunctionReturn(PETSC_SUCCESS);
1529: }