Actual source code: krylovschur.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    SLEPc eigensolver: "krylovschur"

 13:    Method: Krylov-Schur

 15:    Algorithm:

 17:        Single-vector Krylov-Schur method for non-symmetric problems,
 18:        including harmonic extraction.

 20:    References:

 22:        [1] "Krylov-Schur Methods in SLEPc", SLEPc Technical Report STR-7,
 23:            available at https://slepc.upv.es.

 25:        [2] G.W. Stewart, "A Krylov-Schur Algorithm for Large Eigenproblems",
 26:            SIAM J. Matrix Anal. App. 23(3):601-614, 2001.

 28:        [3] "Practical Implementation of Harmonic Krylov-Schur", SLEPc Technical
 29:             Report STR-9, available at https://slepc.upv.es.
 30: */

 32: #include <slepc/private/epsimpl.h>
 33: #include "krylovschur.h"

 35: PetscErrorCode EPSGetArbitraryValues(EPS eps,PetscScalar *rr,PetscScalar *ri)
 36: {
 37:   PetscInt       i,newi,ld,n,l;
 38:   Vec            xr=eps->work[0],xi=eps->work[1];
 39:   PetscScalar    re,im,*Zr,*Zi,*X;

 41:   PetscFunctionBegin;
 42:   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
 43:   PetscCall(DSGetDimensions(eps->ds,&n,&l,NULL,NULL));
 44:   for (i=l;i<n;i++) {
 45:     re = eps->eigr[i];
 46:     im = eps->eigi[i];
 47:     PetscCall(STBackTransform(eps->st,1,&re,&im));
 48:     newi = i;
 49:     PetscCall(DSVectors(eps->ds,DS_MAT_X,&newi,NULL));
 50:     PetscCall(DSGetArray(eps->ds,DS_MAT_X,&X));
 51:     Zr = X+i*ld;
 52:     if (newi==i+1) Zi = X+newi*ld;
 53:     else Zi = NULL;
 54:     PetscCall(EPSComputeRitzVector(eps,Zr,Zi,eps->V,xr,xi));
 55:     PetscCall(DSRestoreArray(eps->ds,DS_MAT_X,&X));
 56:     PetscCall((*eps->arbitrary)(re,im,xr,xi,rr+i,ri+i,eps->arbitraryctx));
 57:   }
 58:   PetscFunctionReturn(PETSC_SUCCESS);
 59: }

 61: static PetscErrorCode EPSSetUp_KrylovSchur_Filter(EPS eps)
 62: {
 63:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
 64:   PetscBool       estimaterange=PETSC_TRUE;
 65:   PetscReal       rleft,rright;
 66:   Mat             A;

 68:   PetscFunctionBegin;
 69:   EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
 70:   EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
 71:   PetscCheck(eps->intb<PETSC_MAX_REAL || eps->inta>PETSC_MIN_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The defined computational interval should have at least one of their sides bounded");
 72:   EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION | EPS_FEATURE_THRESHOLD,PETSC_TRUE," with polynomial filter");
 73:   if (eps->tol==(PetscReal)PETSC_DETERMINE) eps->tol = SLEPC_DEFAULT_TOL*1e-2;  /* use tighter tolerance */
 74:   PetscCall(STFilterSetInterval(eps->st,eps->inta,eps->intb));
 75:   if (!ctx->estimatedrange) {
 76:     PetscCall(STFilterGetRange(eps->st,&rleft,&rright));
 77:     estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
 78:   }
 79:   if (estimaterange) { /* user did not set a range */
 80:     PetscCall(STGetMatrix(eps->st,0,&A));
 81:     PetscCall(MatEstimateSpectralRange_EPS(A,&rleft,&rright));
 82:     PetscCall(PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright));
 83:     PetscCall(STFilterSetRange(eps->st,rleft,rright));
 84:     ctx->estimatedrange = PETSC_TRUE;
 85:   }
 86:   if (eps->ncv==PETSC_DETERMINE && eps->nev==0) eps->nev = 40;  /* user did not provide nev estimation */
 87:   PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
 88:   PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
 89:   if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

 93: static PetscErrorCode EPSSetUp_KrylovSchur(EPS eps)
 94: {
 95:   PetscReal         eta;
 96:   PetscBool         isfilt=PETSC_FALSE;
 97:   BVOrthogType      otype;
 98:   BVOrthogBlockType obtype;
 99:   EPS_KRYLOVSCHUR   *ctx = (EPS_KRYLOVSCHUR*)eps->data;
100:   enum { EPS_KS_DEFAULT,EPS_KS_SYMM,EPS_KS_SLICE,EPS_KS_FILTER,EPS_KS_INDEF,EPS_KS_TWOSIDED } variant;

102:   PetscFunctionBegin;
103:   if (eps->which==EPS_ALL) {  /* default values in case of spectrum slicing or polynomial filter  */
104:     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
105:     if (isfilt) PetscCall(EPSSetUp_KrylovSchur_Filter(eps));
106:     else PetscCall(EPSSetUp_KrylovSchur_Slice(eps));
107:   } else if (eps->isstructured) {
108:     if (eps->problem_type==EPS_BSE) PetscCall(EPSSetUp_KrylovSchur_BSE(eps));
109:     else if (eps->problem_type==EPS_HAMILT) {
110:       PetscCheck(!PetscDefined(USE_COMPLEX),PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"The Hamiltonian Krylov-Schur eigensolver is not yet implemented for complex scalars");
111:       PetscCall(EPSSetUp_KrylovSchur_Hamilt(eps));
112:     } else SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unknown matrix structure");
113:     PetscFunctionReturn(PETSC_SUCCESS);
114:   } else {
115:     PetscCall(EPSSetDimensions_Default(eps,&eps->nev,&eps->ncv,&eps->mpd));
116:     PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
117:     if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv)*((eps->stop==EPS_STOP_THRESHOLD)?10:1);
118:     if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
119:   }
120:   PetscCheck(ctx->lock || eps->mpd>=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");

122:   EPSCheckDefiniteCondition(eps,eps->arbitrary," with arbitrary selection of eigenpairs");

124:   PetscCheck(eps->extraction==EPS_RITZ || eps->extraction==EPS_HARMONIC,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");

126:   if (!ctx->keep) ctx->keep = 0.5;

128:   PetscCall(EPSAllocateSolution(eps,1));
129:   PetscCall(EPS_SetInnerProduct(eps));
130:   if (eps->arbitrary) PetscCall(EPSSetWorkVecs(eps,2));
131:   else if (eps->ishermitian && !eps->ispositive) PetscCall(EPSSetWorkVecs(eps,1));

133:   /* dispatch solve method */
134:   if (eps->ishermitian) {
135:     if (eps->which==EPS_ALL) {
136:       EPSCheckDefiniteCondition(eps,eps->which==EPS_ALL," with spectrum slicing");
137:       variant = isfilt? EPS_KS_FILTER: EPS_KS_SLICE;
138:     } else if (eps->isgeneralized && !eps->ispositive) {
139:       variant = EPS_KS_INDEF;
140:     } else {
141:       switch (eps->extraction) {
142:         case EPS_RITZ:     variant = EPS_KS_SYMM; break;
143:         case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
144:         default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
145:       }
146:     }
147:   } else if (eps->twosided) {
148:     variant = EPS_KS_TWOSIDED;
149:   } else {
150:     switch (eps->extraction) {
151:       case EPS_RITZ:     variant = EPS_KS_DEFAULT; break;
152:       case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
153:       default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
154:     }
155:   }
156:   switch (variant) {
157:     case EPS_KS_DEFAULT:
158:       eps->ops->solve = EPSSolve_KrylovSchur_Default;
159:       eps->ops->computevectors = EPSComputeVectors_Schur;
160:       PetscCall(DSSetType(eps->ds,DSNHEP));
161:       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
162:       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
163:       break;
164:     case EPS_KS_SYMM:
165:     case EPS_KS_FILTER:
166:       eps->ops->solve = EPSSolve_KrylovSchur_Default;
167:       eps->ops->computevectors = EPSComputeVectors_Hermitian;
168:       PetscCall(DSSetType(eps->ds,DSHEP));
169:       PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
170:       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
171:       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
172:       break;
173:     case EPS_KS_SLICE:
174:       eps->ops->solve = EPSSolve_KrylovSchur_Slice;
175:       eps->ops->computevectors = EPSComputeVectors_Slice;
176:       break;
177:     case EPS_KS_INDEF:
178:       eps->ops->solve = EPSSolve_KrylovSchur_Indefinite;
179:       eps->ops->computevectors = EPSComputeVectors_Indefinite;
180:       PetscCall(DSSetType(eps->ds,DSGHIEP));
181:       PetscCall(DSSetCompact(eps->ds,PETSC_TRUE));
182:       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
183:       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
184:       /* force reorthogonalization for pseudo-Lanczos */
185:       PetscCall(BVGetOrthogonalization(eps->V,&otype,NULL,&eta,&obtype));
186:       PetscCall(BVSetOrthogonalization(eps->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype));
187:       break;
188:     case EPS_KS_TWOSIDED:
189:       eps->ops->solve = EPSSolve_KrylovSchur_TwoSided;
190:       eps->ops->computevectors = EPSComputeVectors_Schur;
191:       PetscCall(DSSetType(eps->ds,DSNHEPTS));
192:       PetscCall(DSAllocate(eps->ds,eps->ncv+1));
193:       PetscCall(DSSetExtraRow(eps->ds,PETSC_TRUE));
194:       break;
195:     default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_PLIB,"Unexpected error");
196:   }
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }

200: static PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
201: {
202:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
203:   SlepcSC         sc;
204:   PetscBool       isfilt;

206:   PetscFunctionBegin;
207:   PetscCall(EPSSetUpSort_Default(eps));
208:   if (eps->which==EPS_ALL) {
209:     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
210:     if (isfilt) {
211:       PetscCall(DSGetSlepcSC(eps->ds,&sc));
212:       sc->rg            = NULL;
213:       sc->comparison    = SlepcCompareLargestReal;
214:       sc->comparisonctx = NULL;
215:       sc->map           = NULL;
216:       sc->mapobj        = NULL;
217:     } else {
218:       if (!ctx->global && ctx->sr->numEigs>0) {
219:         PetscCall(DSGetSlepcSC(eps->ds,&sc));
220:         sc->rg            = NULL;
221:         sc->comparison    = SlepcCompareLargestMagnitude;
222:         sc->comparisonctx = NULL;
223:         sc->map           = NULL;
224:         sc->mapobj        = NULL;
225:       }
226:     }
227:   }
228:   PetscFunctionReturn(PETSC_SUCCESS);
229: }

231: PetscErrorCode EPSSolve_KrylovSchur_Default(EPS eps)
232: {
233:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
234:   PetscInt        i,j,*pj,k,l,nv,ld,nconv;
235:   Mat             U,Op,H,T;
236:   PetscScalar     *g;
237:   PetscReal       beta,gamma=1.0;
238:   PetscBool       breakdown,harmonic,hermitian;

240:   PetscFunctionBegin;
241:   PetscCall(DSGetLeadingDimension(eps->ds,&ld));
242:   harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
243:   hermitian = (eps->ishermitian && !harmonic)?PETSC_TRUE:PETSC_FALSE;
244:   if (harmonic) PetscCall(PetscMalloc1(ld,&g));
245:   if (eps->arbitrary) pj = &j;
246:   else pj = NULL;

248:   /* Get the starting Arnoldi vector */
249:   PetscCall(EPSGetStartVector(eps,0,NULL));
250:   l = 0;

252:   /* Restart loop */
253:   while (eps->reason == EPS_CONVERGED_ITERATING) {
254:     eps->its++;

256:     /* Compute an nv-step Arnoldi factorization */
257:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
258:     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
259:     PetscCall(STGetOperator(eps->st,&Op));
260:     if (hermitian) {
261:       PetscCall(DSGetMat(eps->ds,DS_MAT_T,&T));
262:       PetscCall(BVMatLanczos(eps->V,Op,T,eps->nconv+l,&nv,&beta,&breakdown));
263:       PetscCall(DSRestoreMat(eps->ds,DS_MAT_T,&T));
264:     } else {
265:       PetscCall(DSGetMat(eps->ds,DS_MAT_A,&H));
266:       PetscCall(BVMatArnoldi(eps->V,Op,H,eps->nconv+l,&nv,&beta,&breakdown));
267:       PetscCall(DSRestoreMat(eps->ds,DS_MAT_A,&H));
268:     }
269:     PetscCall(STRestoreOperator(eps->st,&Op));
270:     PetscCall(DSSetDimensions(eps->ds,nv,eps->nconv,eps->nconv+l));
271:     PetscCall(DSSetState(eps->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
272:     PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));

274:     /* Compute translation of Krylov decomposition if harmonic extraction used */
275:     if (PetscUnlikely(harmonic)) PetscCall(DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,g,&gamma));

277:     /* Solve projected problem */
278:     PetscCall(DSSolve(eps->ds,eps->eigr,eps->eigi));
279:     if (PetscUnlikely(eps->arbitrary)) {
280:       PetscCall(EPSGetArbitraryValues(eps,eps->rr,eps->ri));
281:       j=1;
282:     }
283:     PetscCall(DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj));
284:     PetscCall(DSUpdateExtraRow(eps->ds));
285:     PetscCall(DSSynchronize(eps->ds,eps->eigr,eps->eigi));

287:     /* Check convergence */
288:     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
289:     EPSSetCtxThreshold(eps,eps->eigr,eps->eigi,k);
290:     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
291:     nconv = k;

293:     /* Update l */
294:     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
295:     else {
296:       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
297:       if (!hermitian) PetscCall(DSGetTruncateSize(eps->ds,k,nv,&l));
298:     }
299:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
300:     if (l) PetscCall(PetscInfo(eps,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));

302:     if (eps->reason == EPS_CONVERGED_ITERATING) {
303:       if (PetscUnlikely(breakdown || k==nv)) {
304:         /* Start a new Arnoldi factorization */
305:         PetscCall(PetscInfo(eps,"Breakdown in Krylov-Schur method (it=%" PetscInt_FMT " norm=%g)\n",eps->its,(double)beta));
306:         if (k<eps->nev) {
307:           PetscCall(EPSGetStartVector(eps,k,&breakdown));
308:           if (breakdown) {
309:             eps->reason = EPS_DIVERGED_BREAKDOWN;
310:             PetscCall(PetscInfo(eps,"Unable to generate more start vectors\n"));
311:           }
312:         }
313:       } else {
314:         /* Undo translation of Krylov decomposition */
315:         if (PetscUnlikely(harmonic)) {
316:           PetscCall(DSSetDimensions(eps->ds,nv,k,l));
317:           PetscCall(DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma));
318:           /* gamma u^ = u - U*g~ */
319:           PetscCall(BVSetActiveColumns(eps->V,0,nv));
320:           PetscCall(BVMultColumn(eps->V,-1.0,1.0,nv,g));
321:           PetscCall(BVScaleColumn(eps->V,nv,1.0/gamma));
322:           PetscCall(BVSetActiveColumns(eps->V,eps->nconv,nv));
323:           PetscCall(DSSetDimensions(eps->ds,nv,k,nv));
324:         }
325:         /* Prepare the Rayleigh quotient for restart */
326:         PetscCall(DSTruncate(eps->ds,k+l,PETSC_FALSE));
327:       }
328:     }
329:     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
330:     PetscCall(DSGetMat(eps->ds,DS_MAT_Q,&U));
331:     PetscCall(BVMultInPlace(eps->V,U,eps->nconv,k+l));
332:     PetscCall(DSRestoreMat(eps->ds,DS_MAT_Q,&U));

334:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
335:       PetscCall(BVCopyColumn(eps->V,nv,k+l));  /* copy restart vector from the last column */
336:       if (eps->stop==EPS_STOP_THRESHOLD && nv-k<5) {  /* reallocate */
337:         eps->ncv = eps->mpd+k;
338:         PetscCall(EPSReallocateSolution(eps,eps->ncv+1));
339:         for (i=nv;i<eps->ncv;i++) eps->perm[i] = i;
340:         PetscCall(DSReallocate(eps->ds,eps->ncv+1));
341:         PetscCall(DSGetLeadingDimension(eps->ds,&ld));
342:       }
343:     }

345:     eps->nconv = k;
346:     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
347:   }

349:   if (harmonic) PetscCall(PetscFree(g));
350:   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
351:   PetscFunctionReturn(PETSC_SUCCESS);
352: }

354: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
355: {
356:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

358:   PetscFunctionBegin;
359:   if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
360:   else {
361:     PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument %g must be in the range [0.1,0.9]",(double)keep);
362:     ctx->keep = keep;
363:   }
364:   PetscFunctionReturn(PETSC_SUCCESS);
365: }

367: /*@
368:    EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
369:    method, in particular the proportion of basis vectors that must be kept
370:    after restart.

372:    Logically Collective

374:    Input Parameters:
375: +  eps - the eigenproblem solver context
376: -  keep - the number of vectors to be kept at restart

378:    Options Database Key:
379: .  -eps_krylovschur_restart - Sets the restart parameter

381:    Notes:
382:    Allowed values are in the range [0.1,0.9]. The default is 0.5.

384:    Level: advanced

386: .seealso: `EPSKrylovSchurGetRestart()`
387: @*/
388: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
389: {
390:   PetscFunctionBegin;
393:   PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
394:   PetscFunctionReturn(PETSC_SUCCESS);
395: }

397: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
398: {
399:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

401:   PetscFunctionBegin;
402:   *keep = ctx->keep;
403:   PetscFunctionReturn(PETSC_SUCCESS);
404: }

406: /*@
407:    EPSKrylovSchurGetRestart - Gets the restart parameter used in the
408:    Krylov-Schur method.

410:    Not Collective

412:    Input Parameter:
413: .  eps - the eigenproblem solver context

415:    Output Parameter:
416: .  keep - the restart parameter

418:    Level: advanced

420: .seealso: `EPSKrylovSchurSetRestart()`
421: @*/
422: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
423: {
424:   PetscFunctionBegin;
426:   PetscAssertPointer(keep,2);
427:   PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
428:   PetscFunctionReturn(PETSC_SUCCESS);
429: }

431: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
432: {
433:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

435:   PetscFunctionBegin;
436:   ctx->lock = lock;
437:   PetscFunctionReturn(PETSC_SUCCESS);
438: }

440: /*@
441:    EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
442:    the Krylov-Schur method.

444:    Logically Collective

446:    Input Parameters:
447: +  eps  - the eigenproblem solver context
448: -  lock - true if the locking variant must be selected

450:    Options Database Key:
451: .  -eps_krylovschur_locking - Sets the locking flag

453:    Notes:
454:    The default is to lock converged eigenpairs when the method restarts.
455:    This behaviour can be changed so that all directions are kept in the
456:    working subspace even if already converged to working accuracy (the
457:    non-locking variant).

459:    Level: advanced

461: .seealso: `EPSKrylovSchurGetLocking()`
462: @*/
463: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
464: {
465:   PetscFunctionBegin;
468:   PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
469:   PetscFunctionReturn(PETSC_SUCCESS);
470: }

472: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
473: {
474:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

476:   PetscFunctionBegin;
477:   *lock = ctx->lock;
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: /*@
482:    EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
483:    method.

485:    Not Collective

487:    Input Parameter:
488: .  eps - the eigenproblem solver context

490:    Output Parameter:
491: .  lock - the locking flag

493:    Level: advanced

495: .seealso: `EPSKrylovSchurSetLocking()`
496: @*/
497: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
498: {
499:   PetscFunctionBegin;
501:   PetscAssertPointer(lock,2);
502:   PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
503:   PetscFunctionReturn(PETSC_SUCCESS);
504: }

506: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
507: {
508:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
509:   PetscMPIInt     size;
510:   PetscInt        newnpart;

512:   PetscFunctionBegin;
513:   if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
514:     newnpart = 1;
515:   } else {
516:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
517:     PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
518:     newnpart = npart;
519:   }
520:   if (ctx->npart!=newnpart) {
521:     if (ctx->npart>1) {
522:       PetscCall(PetscSubcommDestroy(&ctx->subc));
523:       if (ctx->commset) {
524:         PetscCallMPI(MPI_Comm_free(&ctx->commrank));
525:         ctx->commset = PETSC_FALSE;
526:       }
527:     }
528:     PetscCall(EPSDestroy(&ctx->eps));
529:     ctx->npart = newnpart;
530:     eps->state = EPS_STATE_INITIAL;
531:   }
532:   PetscFunctionReturn(PETSC_SUCCESS);
533: }

535: /*@
536:    EPSKrylovSchurSetPartitions - Sets the number of partitions for the
537:    case of doing spectrum slicing for a computational interval with the
538:    communicator split in several sub-communicators.

540:    Logically Collective

542:    Input Parameters:
543: +  eps   - the eigenproblem solver context
544: -  npart - number of partitions

546:    Options Database Key:
547: .  -eps_krylovschur_partitions <npart> - Sets the number of partitions

549:    Notes:
550:    By default, npart=1 so all processes in the communicator participate in
551:    the processing of the whole interval. If npart>1 then the interval is
552:    divided into npart subintervals, each of them being processed by a
553:    subset of processes.

555:    The interval is split proportionally unless the separation points are
556:    specified with EPSKrylovSchurSetSubintervals().

558:    Level: advanced

560: .seealso: `EPSKrylovSchurSetSubintervals()`, `EPSSetInterval()`
561: @*/
562: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
563: {
564:   PetscFunctionBegin;
567:   PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
568:   PetscFunctionReturn(PETSC_SUCCESS);
569: }

571: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
572: {
573:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

575:   PetscFunctionBegin;
576:   *npart  = ctx->npart;
577:   PetscFunctionReturn(PETSC_SUCCESS);
578: }

580: /*@
581:    EPSKrylovSchurGetPartitions - Gets the number of partitions of the
582:    communicator in case of spectrum slicing.

584:    Not Collective

586:    Input Parameter:
587: .  eps - the eigenproblem solver context

589:    Output Parameter:
590: .  npart - number of partitions

592:    Level: advanced

594: .seealso: `EPSKrylovSchurSetPartitions()`
595: @*/
596: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
597: {
598:   PetscFunctionBegin;
600:   PetscAssertPointer(npart,2);
601:   PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
602:   PetscFunctionReturn(PETSC_SUCCESS);
603: }

605: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
606: {
607:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

609:   PetscFunctionBegin;
610:   ctx->detect = detect;
611:   eps->state  = EPS_STATE_INITIAL;
612:   PetscFunctionReturn(PETSC_SUCCESS);
613: }

615: /*@
616:    EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
617:    zeros during the factorizations throughout the spectrum slicing computation.

619:    Logically Collective

621:    Input Parameters:
622: +  eps    - the eigenproblem solver context
623: -  detect - check for zeros

625:    Options Database Key:
626: .  -eps_krylovschur_detect_zeros - Check for zeros; this takes an optional
627:    bool value (0/1/no/yes/true/false)

629:    Notes:
630:    A zero in the factorization indicates that a shift coincides with an eigenvalue.

632:    This flag is turned off by default, and may be necessary in some cases,
633:    especially when several partitions are being used. This feature currently
634:    requires an external package for factorizations with support for zero
635:    detection, e.g. MUMPS.

637:    Level: advanced

639: .seealso: `EPSKrylovSchurSetPartitions()`, `EPSSetInterval()`
640: @*/
641: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
642: {
643:   PetscFunctionBegin;
646:   PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
647:   PetscFunctionReturn(PETSC_SUCCESS);
648: }

650: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
651: {
652:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

654:   PetscFunctionBegin;
655:   *detect = ctx->detect;
656:   PetscFunctionReturn(PETSC_SUCCESS);
657: }

659: /*@
660:    EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
661:    in spectrum slicing.

663:    Not Collective

665:    Input Parameter:
666: .  eps - the eigenproblem solver context

668:    Output Parameter:
669: .  detect - whether zeros detection is enforced during factorizations

671:    Level: advanced

673: .seealso: `EPSKrylovSchurSetDetectZeros()`
674: @*/
675: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
676: {
677:   PetscFunctionBegin;
679:   PetscAssertPointer(detect,2);
680:   PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
681:   PetscFunctionReturn(PETSC_SUCCESS);
682: }

684: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
685: {
686:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

688:   PetscFunctionBegin;
689:   if (nev != PETSC_CURRENT) {
690:     PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
691:     ctx->nev = nev;
692:   }
693:   if (ncv == PETSC_DETERMINE) {
694:     ctx->ncv = PETSC_DETERMINE;
695:   } else if (ncv != PETSC_CURRENT) {
696:     PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
697:     ctx->ncv = ncv;
698:   }
699:   if (mpd == PETSC_DETERMINE) {
700:     ctx->mpd = PETSC_DETERMINE;
701:   } else if (mpd != PETSC_CURRENT) {
702:     PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
703:     ctx->mpd = mpd;
704:   }
705:   eps->state = EPS_STATE_INITIAL;
706:   PetscFunctionReturn(PETSC_SUCCESS);
707: }

709: /*@
710:    EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
711:    step in case of doing spectrum slicing for a computational interval.
712:    The meaning of the parameters is the same as in EPSSetDimensions().

714:    Logically Collective

716:    Input Parameters:
717: +  eps - the eigenproblem solver context
718: .  nev - number of eigenvalues to compute
719: .  ncv - the maximum dimension of the subspace to be used by the subsolve
720: -  mpd - the maximum dimension allowed for the projected problem

722:    Options Database Key:
723: +  -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
724: .  -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
725: -  -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension

727:    Note:
728:    Use PETSC_DETERMINE for ncv and mpd to assign a default value. For any
729:    of the arguments, use PETSC_CURRENT to preserve the current value.

731:    Level: advanced

733: .seealso: `EPSKrylovSchurGetDimensions()`, `EPSSetDimensions()`, `EPSSetInterval()`
734: @*/
735: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
736: {
737:   PetscFunctionBegin;
742:   PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
743:   PetscFunctionReturn(PETSC_SUCCESS);
744: }

746: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
747: {
748:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

750:   PetscFunctionBegin;
751:   if (nev) *nev = ctx->nev;
752:   if (ncv) *ncv = ctx->ncv;
753:   if (mpd) *mpd = ctx->mpd;
754:   PetscFunctionReturn(PETSC_SUCCESS);
755: }

757: /*@
758:    EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
759:    step in case of doing spectrum slicing for a computational interval.

761:    Not Collective

763:    Input Parameter:
764: .  eps - the eigenproblem solver context

766:    Output Parameters:
767: +  nev - number of eigenvalues to compute
768: .  ncv - the maximum dimension of the subspace to be used by the subsolve
769: -  mpd - the maximum dimension allowed for the projected problem

771:    Level: advanced

773: .seealso: `EPSKrylovSchurSetDimensions()`
774: @*/
775: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
776: {
777:   PetscFunctionBegin;
779:   PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
780:   PetscFunctionReturn(PETSC_SUCCESS);
781: }

783: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
784: {
785:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
786:   PetscInt        i;

788:   PetscFunctionBegin;
789:   PetscCheck(subint[0]==eps->inta && subint[ctx->npart]==eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"First and last values must match the endpoints of EPSSetInterval()");
790:   for (i=0;i<ctx->npart;i++) PetscCheck(subint[i]<=subint[i+1],PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Array must contain values in ascending order");
791:   if (ctx->subintervals) PetscCall(PetscFree(ctx->subintervals));
792:   PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
793:   for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
794:   ctx->subintset = PETSC_TRUE;
795:   eps->state = EPS_STATE_INITIAL;
796:   PetscFunctionReturn(PETSC_SUCCESS);
797: }

799: /*@
800:    EPSKrylovSchurSetSubintervals - Sets the points that delimit the
801:    subintervals to be used in spectrum slicing with several partitions.

803:    Logically Collective

805:    Input Parameters:
806: +  eps    - the eigenproblem solver context
807: -  subint - array of real values specifying subintervals

809:    Notes:
810:    This function must be called after EPSKrylovSchurSetPartitions(). For npart
811:    partitions, the argument subint must contain npart+1 real values sorted in
812:    ascending order, subint_0, subint_1, ..., subint_npart, where the first
813:    and last values must coincide with the interval endpoints set with
814:    EPSSetInterval().

816:    The subintervals are then defined by two consecutive points [subint_0,subint_1],
817:    [subint_1,subint_2], and so on.

819:    Level: advanced

821: .seealso: `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubintervals()`, `EPSSetInterval()`
822: @*/
823: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal subint[])
824: {
825:   PetscFunctionBegin;
827:   PetscAssertPointer(subint,2);
828:   PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
829:   PetscFunctionReturn(PETSC_SUCCESS);
830: }

832: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal *subint[])
833: {
834:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
835:   PetscInt        i;

837:   PetscFunctionBegin;
838:   if (!ctx->subintset) {
839:     PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
840:     PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
841:   }
842:   PetscCall(PetscMalloc1(ctx->npart+1,subint));
843:   for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
844:   PetscFunctionReturn(PETSC_SUCCESS);
845: }

847: /*@C
848:    EPSKrylovSchurGetSubintervals - Returns the points that delimit the
849:    subintervals used in spectrum slicing with several partitions.

851:    Not Collective

853:    Input Parameter:
854: .  eps    - the eigenproblem solver context

856:    Output Parameter:
857: .  subint - array of real values specifying subintervals

859:    Notes:
860:    If the user passed values with EPSKrylovSchurSetSubintervals(), then the
861:    same values are returned. Otherwise, the values computed internally are
862:    obtained.

864:    This function is only available for spectrum slicing runs.

866:    The returned array has length npart+1 (see EPSKrylovSchurGetPartitions())
867:    and should be freed by the user.

869:    Fortran Notes:
870:    The calling sequence from Fortran is
871: .vb
872:    EPSKrylovSchurGetSubintervals(eps,subint,ierr)
873:    double precision subint(npart+1) output
874: .ve

876:    Level: advanced

878: .seealso: `EPSKrylovSchurSetSubintervals()`, `EPSKrylovSchurGetPartitions()`, `EPSSetInterval()`
879: @*/
880: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal *subint[]) PeNS
881: {
882:   PetscFunctionBegin;
884:   PetscAssertPointer(subint,2);
885:   PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
886:   PetscFunctionReturn(PETSC_SUCCESS);
887: }

889: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[])
890: {
891:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
892:   PetscInt        i,numsh;
893:   EPS_SR          sr = ctx->sr;

895:   PetscFunctionBegin;
896:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
897:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
898:   switch (eps->state) {
899:   case EPS_STATE_INITIAL:
900:     break;
901:   case EPS_STATE_SETUP:
902:     numsh = ctx->npart+1;
903:     if (n) *n = numsh;
904:     if (shifts) {
905:       PetscCall(PetscMalloc1(numsh,shifts));
906:       (*shifts)[0] = eps->inta;
907:       if (ctx->npart==1) (*shifts)[1] = eps->intb;
908:       else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
909:     }
910:     if (inertias) {
911:       PetscCall(PetscMalloc1(numsh,inertias));
912:       (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
913:       if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
914:       else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
915:     }
916:     break;
917:   case EPS_STATE_SOLVED:
918:   case EPS_STATE_EIGENVECTORS:
919:     numsh = ctx->nshifts;
920:     if (n) *n = numsh;
921:     if (shifts) {
922:       PetscCall(PetscMalloc1(numsh,shifts));
923:       for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
924:     }
925:     if (inertias) {
926:       PetscCall(PetscMalloc1(numsh,inertias));
927:       for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
928:     }
929:     break;
930:   }
931:   PetscFunctionReturn(PETSC_SUCCESS);
932: }

934: /*@C
935:    EPSKrylovSchurGetInertias - Gets the values of the shifts and their
936:    corresponding inertias in case of doing spectrum slicing for a
937:    computational interval.

939:    Not Collective

941:    Input Parameter:
942: .  eps - the eigenproblem solver context

944:    Output Parameters:
945: +  n        - number of shifts, including the endpoints of the interval
946: .  shifts   - the values of the shifts used internally in the solver
947: -  inertias - the values of the inertia in each shift

949:    Notes:
950:    If called after EPSSolve(), all shifts used internally by the solver are
951:    returned (including both endpoints and any intermediate ones). If called
952:    before EPSSolve() and after EPSSetUp() then only the information of the
953:    endpoints of subintervals is available.

955:    This function is only available for spectrum slicing runs.

957:    The returned arrays should be freed by the user. Can pass NULL in any of
958:    the two arrays if not required.

960:    Fortran Notes:
961:    The calling sequence from Fortran is
962: .vb
963:    EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
964:    integer n
965:    double precision shifts(*)
966:    integer inertias(*)
967: .ve
968:    The arrays should be at least of length n. The value of n can be determined
969:    by an initial call
970: .vb
971:    EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr)
972: .ve

974:    Level: advanced

976: .seealso: `EPSSetInterval()`, `EPSKrylovSchurSetSubintervals()`
977: @*/
978: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[]) PeNS
979: {
980:   PetscFunctionBegin;
982:   PetscAssertPointer(n,2);
983:   PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
984:   PetscFunctionReturn(PETSC_SUCCESS);
985: }

987: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
988: {
989:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
990:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

992:   PetscFunctionBegin;
993:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
994:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
995:   if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
996:   if (n) *n = sr->numEigs;
997:   if (v) PetscCall(BVCreateVec(sr->V,v));
998:   PetscFunctionReturn(PETSC_SUCCESS);
999: }

1001: /*@
1002:    EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
1003:    doing spectrum slicing for a computational interval with multiple
1004:    communicators.

1006:    Collective on the subcommunicator (if v is given)

1008:    Input Parameter:
1009: .  eps - the eigenproblem solver context

1011:    Output Parameters:
1012: +  k - index of the subinterval for the calling process
1013: .  n - number of eigenvalues found in the k-th subinterval
1014: -  v - a vector owned by processes in the subcommunicator with dimensions
1015:        compatible for locally computed eigenvectors (or NULL)

1017:    Notes:
1018:    This function is only available for spectrum slicing runs.

1020:    The returned Vec should be destroyed by the user.

1022:    Level: advanced

1024: .seealso: `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommPairs()`
1025: @*/
1026: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1027: {
1028:   PetscFunctionBegin;
1030:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1031:   PetscFunctionReturn(PETSC_SUCCESS);
1032: }

1034: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1035: {
1036:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1037:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

1039:   PetscFunctionBegin;
1040:   EPSCheckSolved(eps,1);
1041:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1042:   PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1043:   if (eig) *eig = sr->eigr[sr->perm[i]];
1044:   if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1045:   PetscFunctionReturn(PETSC_SUCCESS);
1046: }

1048: /*@
1049:    EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
1050:    internally in the subcommunicator to which the calling process belongs.

1052:    Collective on the subcommunicator (if v is given)

1054:    Input Parameters:
1055: +  eps - the eigenproblem solver context
1056: -  i   - index of the solution

1058:    Output Parameters:
1059: +  eig - the eigenvalue
1060: -  v   - the eigenvector

1062:    Notes:
1063:    It is allowed to pass NULL for v if the eigenvector is not required.
1064:    Otherwise, the caller must provide a valid Vec objects, i.e.,
1065:    it must be created by the calling program with EPSKrylovSchurGetSubcommInfo().

1067:    The index i should be a value between 0 and n-1, where n is the number of
1068:    vectors in the local subinterval, see EPSKrylovSchurGetSubcommInfo().

1070:    Level: advanced

1072: .seealso: `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`, `EPSKrylovSchurGetSubcommMats()`
1073: @*/
1074: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1075: {
1076:   PetscFunctionBegin;
1079:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1080:   PetscFunctionReturn(PETSC_SUCCESS);
1081: }

1083: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1084: {
1085:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

1087:   PetscFunctionBegin;
1088:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1089:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1090:   PetscCall(EPSGetOperators(ctx->eps,A,B));
1091:   PetscFunctionReturn(PETSC_SUCCESS);
1092: }

1094: /*@
1095:    EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1096:    internally in the subcommunicator to which the calling process belongs.

1098:    Collective on the subcommunicator

1100:    Input Parameter:
1101: .  eps - the eigenproblem solver context

1103:    Output Parameters:
1104: +  A  - the matrix associated with the eigensystem
1105: -  B  - the second matrix in the case of generalized eigenproblems

1107:    Notes:
1108:    This is the analog of EPSGetOperators(), but returns the matrices distributed
1109:    differently (in the subcommunicator rather than in the parent communicator).

1111:    These matrices should not be modified by the user.

1113:    Level: advanced

1115: .seealso: `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommInfo()`
1116: @*/
1117: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1118: {
1119:   PetscFunctionBegin;
1121:   PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1122:   PetscFunctionReturn(PETSC_SUCCESS);
1123: }

1125: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1126: {
1127:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1128:   Mat             A,B=NULL,Ag,Bg=NULL;
1129:   PetscBool       reuse=PETSC_TRUE;

1131:   PetscFunctionBegin;
1132:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1133:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1134:   PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1135:   PetscCall(EPSGetOperators(ctx->eps,&A,&B));

1137:   PetscCall(MatScale(A,a));
1138:   if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1139:   if (B) PetscCall(MatScale(B,b));
1140:   if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1141:   PetscCall(EPSSetOperators(ctx->eps,A,B));

1143:   /* Update stored matrix state */
1144:   subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1145:   PetscCall(MatGetState(A,&subctx->Astate));
1146:   if (B) PetscCall(MatGetState(B,&subctx->Bstate));

1148:   /* Update matrices in the parent communicator if requested by user */
1149:   if (globalup) {
1150:     if (ctx->npart>1) {
1151:       if (!ctx->isrow) {
1152:         PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1153:         reuse = PETSC_FALSE;
1154:       }
1155:       if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1156:       if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1157:       PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1158:       PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1159:       if (B) {
1160:         if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1161:         PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1162:         PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1163:       }
1164:     }
1165:     PetscCall(MatGetState(Ag,&ctx->Astate));
1166:     if (Bg) PetscCall(MatGetState(Bg,&ctx->Bstate));
1167:   }
1168:   PetscCall(EPSSetOperators(eps,Ag,Bg));
1169:   PetscFunctionReturn(PETSC_SUCCESS);
1170: }

1172: /*@
1173:    EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1174:    internally in the subcommunicator to which the calling process belongs.

1176:    Collective

1178:    Input Parameters:
1179: +  eps - the eigenproblem solver context
1180: .  s   - scalar that multiplies the existing A matrix
1181: .  a   - scalar used in the axpy operation on A
1182: .  Au  - matrix used in the axpy operation on A
1183: .  t   - scalar that multiplies the existing B matrix
1184: .  b   - scalar used in the axpy operation on B
1185: .  Bu  - matrix used in the axpy operation on B
1186: .  str - structure flag
1187: -  globalup - flag indicating if global matrices must be updated

1189:    Notes:
1190:    This function modifies the eigenproblem matrices at the subcommunicator level,
1191:    and optionally updates the global matrices in the parent communicator. The updates
1192:    are expressed as A <-- s*A + a*Au,  B <-- t*B + b*Bu.

1194:    It is possible to update one of the matrices, or both.

1196:    The matrices Au and Bu must be equal in all subcommunicators.

1198:    The str flag is passed to the MatAXPY() operations to perform the updates.

1200:    If globalup is true, communication is carried out to reconstruct the updated
1201:    matrices in the parent communicator. The user must be warned that if global
1202:    matrices are not in sync with subcommunicator matrices, the errors computed
1203:    by EPSComputeError() will be wrong even if the computed solution is correct
1204:    (the synchronization may be done only once at the end).

1206:    Level: advanced

1208: .seealso: `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`, `EPSKrylovSchurGetSubcommMats()`
1209: @*/
1210: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1211: {
1212:   PetscFunctionBegin;
1222:   PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1223:   PetscFunctionReturn(PETSC_SUCCESS);
1224: }

1226: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1227: {
1228:   EPS_KRYLOVSCHUR  *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1229:   Mat              A,B=NULL,Ar=NULL,Br=NULL;
1230:   PetscMPIInt      rank;
1231:   PetscObjectState Astate,Bstate=0;
1232:   PetscObjectId    Aid,Bid=0;
1233:   STType           sttype;
1234:   PetscInt         nmat;
1235:   const char       *prefix;
1236:   MPI_Comm         child;

1238:   PetscFunctionBegin;
1239:   PetscCall(EPSGetOperators(eps,&A,&B));
1240:   if (ctx->npart==1) {
1241:     if (!ctx->eps) PetscCall(EPSCreate(((PetscObject)eps)->comm,&ctx->eps));
1242:     PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1243:     PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1244:     PetscCall(EPSSetOperators(ctx->eps,A,B));
1245:   } else {
1246:     PetscCall(MatGetState(A,&Astate));
1247:     PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1248:     if (B) {
1249:       PetscCall(MatGetState(B,&Bstate));
1250:       PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1251:     }
1252:     if (!ctx->subc) {
1253:       /* Create context for subcommunicators */
1254:       PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1255:       PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1256:       PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1257:       PetscCall(PetscSubcommGetChild(ctx->subc,&child));

1259:       /* Duplicate matrices */
1260:       PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1261:       ctx->Astate = Astate;
1262:       ctx->Aid = Aid;
1263:       PetscCall(MatPropagateSymmetryOptions(A,Ar));
1264:       if (B) {
1265:         PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1266:         ctx->Bstate = Bstate;
1267:         ctx->Bid = Bid;
1268:         PetscCall(MatPropagateSymmetryOptions(B,Br));
1269:       }
1270:     } else {
1271:       PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1272:       if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1273:         PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1274:         if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
1275:         PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1276:         ctx->Astate = Astate;
1277:         ctx->Aid = Aid;
1278:         PetscCall(MatPropagateSymmetryOptions(A,Ar));
1279:         if (B) {
1280:           PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1281:           ctx->Bstate = Bstate;
1282:           ctx->Bid = Bid;
1283:           PetscCall(MatPropagateSymmetryOptions(B,Br));
1284:         }
1285:         PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1286:         PetscCall(MatDestroy(&Ar));
1287:         PetscCall(MatDestroy(&Br));
1288:       }
1289:     }

1291:     /* Create auxiliary EPS */
1292:     if (!ctx->eps) {
1293:       PetscCall(EPSCreate(child,&ctx->eps));
1294:       PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1295:       PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1296:       PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1297:       PetscCall(MatDestroy(&Ar));
1298:       PetscCall(MatDestroy(&Br));
1299:     }
1300:     /* Create subcommunicator grouping processes with same rank */
1301:     if (!ctx->commset) {
1302:       PetscCallMPI(MPI_Comm_rank(child,&rank));
1303:       PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1304:       ctx->commset = PETSC_TRUE;
1305:     }
1306:   }
1307:   PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1308:   PetscCall(STGetType(eps->st,&sttype));
1309:   PetscCall(STSetType(ctx->eps->st,sttype));

1311:   ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1312:   ctx_local->npart = ctx->npart;
1313:   ctx_local->global = PETSC_FALSE;
1314:   ctx_local->eps = eps;
1315:   ctx_local->subc = ctx->subc;
1316:   ctx_local->commrank = ctx->commrank;
1317:   *childeps = ctx->eps;
1318:   PetscFunctionReturn(PETSC_SUCCESS);
1319: }

1321: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1322: {
1323:   EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1324:   ST              st;
1325:   PetscBool       isfilt;

1327:   PetscFunctionBegin;
1328:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1329:   PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1330:   PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1331:   PetscCall(EPSGetST(ctx->eps,&st));
1332:   PetscCall(STGetOperator(st,NULL));
1333:   PetscCall(STGetKSP(st,ksp));
1334:   PetscFunctionReturn(PETSC_SUCCESS);
1335: }

1337: /*@
1338:    EPSKrylovSchurGetKSP - Retrieve the linear solver object associated with the
1339:    internal EPS object in case of doing spectrum slicing for a computational interval.

1341:    Collective

1343:    Input Parameter:
1344: .  eps - the eigenproblem solver context

1346:    Output Parameter:
1347: .  ksp - the internal KSP object

1349:    Notes:
1350:    When invoked to compute all eigenvalues in an interval with spectrum
1351:    slicing, EPSKRYLOVSCHUR creates another EPS object internally that is
1352:    used to compute eigenvalues by chunks near selected shifts. This function
1353:    allows access to the KSP object associated to this internal EPS object.

1355:    This function is only available for spectrum slicing runs. In case of
1356:    having more than one partition, the returned KSP will be different
1357:    in MPI processes belonging to different partitions. Hence, if required,
1358:    EPSKrylovSchurSetPartitions() must be called BEFORE this function.

1360:    Level: advanced

1362: .seealso: `EPSSetInterval()`, `EPSKrylovSchurSetPartitions()`
1363: @*/
1364: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1365: {
1366:   PetscFunctionBegin;
1368:   PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1369:   PetscFunctionReturn(PETSC_SUCCESS);
1370: }

1372: static PetscErrorCode EPSKrylovSchurSetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType bse)
1373: {
1374:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

1376:   PetscFunctionBegin;
1377:   switch (bse) {
1378:     case EPS_KRYLOVSCHUR_BSE_SHAO:
1379:     case EPS_KRYLOVSCHUR_BSE_GRUNING:
1380:     case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
1381:       if (ctx->bse != bse) {
1382:         ctx->bse = bse;
1383:         eps->state = EPS_STATE_INITIAL;
1384:       }
1385:       break;
1386:     default:
1387:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid BSE type");
1388:   }
1389:   PetscFunctionReturn(PETSC_SUCCESS);
1390: }

1392: /*@
1393:    EPSKrylovSchurSetBSEType - Sets the method to be used for BSE structured eigenproblems
1394:    in the Krylov-Schur solver.

1396:    Logically Collective

1398:    Input Parameters:
1399: +  eps - the eigenproblem solver context
1400: -  bse - the BSE method

1402:    Options Database Key:
1403: .  -eps_krylovschur_bse_type - Sets the BSE type (either 'shao', 'gruning', or 'projectedbse')

1405:    Level: advanced

1407: .seealso: `EPSKrylovSchurGetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1408: @*/
1409: PetscErrorCode EPSKrylovSchurSetBSEType(EPS eps,EPSKrylovSchurBSEType bse)
1410: {
1411:   PetscFunctionBegin;
1414:   PetscTryMethod(eps,"EPSKrylovSchurSetBSEType_C",(EPS,EPSKrylovSchurBSEType),(eps,bse));
1415:   PetscFunctionReturn(PETSC_SUCCESS);
1416: }

1418: static PetscErrorCode EPSKrylovSchurGetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType *bse)
1419: {
1420:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

1422:   PetscFunctionBegin;
1423:   *bse = ctx->bse;
1424:   PetscFunctionReturn(PETSC_SUCCESS);
1425: }

1427: /*@
1428:    EPSKrylovSchurGetBSEType - Gets the method used for BSE structured eigenproblems
1429:    in the Krylov-Schur solver.

1431:    Not Collective

1433:    Input Parameter:
1434: .  eps - the eigenproblem solver context

1436:    Output Parameter:
1437: .  bse - the BSE method

1439:    Level: advanced

1441: .seealso: `EPSKrylovSchurSetBSEType()`, `EPSKrylovSchurBSEType`, `MatCreateBSE()`
1442: @*/
1443: PetscErrorCode EPSKrylovSchurGetBSEType(EPS eps,EPSKrylovSchurBSEType *bse)
1444: {
1445:   PetscFunctionBegin;
1447:   PetscAssertPointer(bse,2);
1448:   PetscUseMethod(eps,"EPSKrylovSchurGetBSEType_C",(EPS,EPSKrylovSchurBSEType*),(eps,bse));
1449:   PetscFunctionReturn(PETSC_SUCCESS);
1450: }

1452: static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems PetscOptionsObject)
1453: {
1454:   EPS_KRYLOVSCHUR       *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1455:   PetscBool             flg,lock,b,f1,f2,f3,isfilt;
1456:   PetscReal             keep;
1457:   PetscInt              i,j,k;
1458:   KSP                   ksp;
1459:   EPSKrylovSchurBSEType bse;

1461:   PetscFunctionBegin;
1462:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");

1464:     PetscCall(PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg));
1465:     if (flg) PetscCall(EPSKrylovSchurSetRestart(eps,keep));

1467:     PetscCall(PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg));
1468:     if (flg) PetscCall(EPSKrylovSchurSetLocking(eps,lock));

1470:     i = ctx->npart;
1471:     PetscCall(PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg));
1472:     if (flg) PetscCall(EPSKrylovSchurSetPartitions(eps,i));

1474:     b = ctx->detect;
1475:     PetscCall(PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg));
1476:     if (flg) PetscCall(EPSKrylovSchurSetDetectZeros(eps,b));

1478:     i = 1;
1479:     j = k = PETSC_DECIDE;
1480:     PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1481:     PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1482:     PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1483:     if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));

1485:     PetscCall(PetscOptionsEnum("-eps_krylovschur_bse_type","Method for BSE structured eigenproblems","EPSKrylovSchurSetBSEType",EPSKrylovSchurBSETypes,(PetscEnum)ctx->bse,(PetscEnum*)&bse,&flg));
1486:     if (flg) PetscCall(EPSKrylovSchurSetBSEType(eps,bse));

1488:   PetscOptionsHeadEnd();

1490:   /* set options of child KSP in spectrum slicing */
1491:   if (eps->which==EPS_ALL) {
1492:     if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1493:     PetscCall(EPSSetDefaultST(eps));
1494:     PetscCall(STSetFromOptions(eps->st));  /* need to advance this to check ST type */
1495:     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1496:     if (!isfilt) {
1497:       PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1498:       PetscCall(KSPSetFromOptions(ksp));
1499:     }
1500:   }
1501:   PetscFunctionReturn(PETSC_SUCCESS);
1502: }

1504: static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1505: {
1506:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1507:   PetscBool       isascii,isfilt;
1508:   KSP             ksp;
1509:   PetscViewer     sviewer;

1511:   PetscFunctionBegin;
1512:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1513:   if (isascii) {
1514:     PetscCall(PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1515:     PetscCall(PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-"));
1516:     if (eps->problem_type==EPS_BSE) PetscCall(PetscViewerASCIIPrintf(viewer,"  BSE method: %s\n",EPSKrylovSchurBSETypes[ctx->bse]));
1517:     if (eps->which==EPS_ALL) {
1518:       PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1519:       if (isfilt) PetscCall(PetscViewerASCIIPrintf(viewer,"  using filtering to extract all eigenvalues in an interval\n"));
1520:       else {
1521:         PetscCall(PetscViewerASCIIPrintf(viewer,"  doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd));
1522:         if (ctx->npart>1) {
1523:           PetscCall(PetscViewerASCIIPrintf(viewer,"  multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart));
1524:           if (ctx->detect) PetscCall(PetscViewerASCIIPrintf(viewer,"  detecting zeros when factorizing at subinterval boundaries\n"));
1525:         }
1526:         /* view child KSP */
1527:         PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1528:         PetscCall(PetscViewerASCIIPushTab(viewer));
1529:         if (ctx->npart>1 && ctx->subc) {
1530:           PetscCall(PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer));
1531:           if (!ctx->subc->color) PetscCall(KSPView(ksp,sviewer));
1532:           PetscCall(PetscViewerFlush(sviewer));
1533:           PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer));
1534:           /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1535:           PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1536:         } else PetscCall(KSPView(ksp,viewer));
1537:         PetscCall(PetscViewerASCIIPopTab(viewer));
1538:       }
1539:     }
1540:   }
1541:   PetscFunctionReturn(PETSC_SUCCESS);
1542: }

1544: static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1545: {
1546:   PetscBool      isfilt;

1548:   PetscFunctionBegin;
1549:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1550:   if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1551:   PetscCall(PetscFree(eps->data));
1552:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1553:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1554:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1555:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1556:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1557:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1558:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1559:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1560:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1561:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1562:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1563:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1564:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1565:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1566:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1567:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1568:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1569:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1570:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",NULL));
1571:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",NULL));
1572:   PetscFunctionReturn(PETSC_SUCCESS);
1573: }

1575: static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1576: {
1577:   PetscBool      isfilt;

1579:   PetscFunctionBegin;
1580:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1581:   if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1582:   PetscFunctionReturn(PETSC_SUCCESS);
1583: }

1585: static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1586: {
1587:   PetscFunctionBegin;
1588:   if (eps->which==EPS_ALL) {
1589:     if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1590:   }
1591:   PetscFunctionReturn(PETSC_SUCCESS);
1592: }

1594: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1595: {
1596:   EPS_KRYLOVSCHUR *ctx;

1598:   PetscFunctionBegin;
1599:   PetscCall(PetscNew(&ctx));
1600:   eps->data   = (void*)ctx;
1601:   ctx->lock   = PETSC_TRUE;
1602:   ctx->nev    = 1;
1603:   ctx->ncv    = PETSC_DETERMINE;
1604:   ctx->mpd    = PETSC_DETERMINE;
1605:   ctx->npart  = 1;
1606:   ctx->detect = PETSC_FALSE;
1607:   ctx->global = PETSC_TRUE;

1609:   eps->useds = PETSC_TRUE;

1611:   /* solve and computevectors determined at setup */
1612:   eps->ops->setup          = EPSSetUp_KrylovSchur;
1613:   eps->ops->setupsort      = EPSSetUpSort_KrylovSchur;
1614:   eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1615:   eps->ops->destroy        = EPSDestroy_KrylovSchur;
1616:   eps->ops->reset          = EPSReset_KrylovSchur;
1617:   eps->ops->view           = EPSView_KrylovSchur;
1618:   eps->ops->backtransform  = EPSBackTransform_Default;
1619:   eps->ops->setdefaultst   = EPSSetDefaultST_KrylovSchur;

1621:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1622:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1623:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1624:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1625:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1626:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1627:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1628:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1629:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1630:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1631:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1632:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1633:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1634:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1635:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1636:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1637:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1638:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1639:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",EPSKrylovSchurSetBSEType_KrylovSchur));
1640:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",EPSKrylovSchurGetBSEType_KrylovSchur));
1641:   PetscFunctionReturn(PETSC_SUCCESS);
1642: }