Actual source code: krylovschur.c

slepc-main 2024-11-15
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

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

 13:    Method: Krylov-Schur

 15:    Algorithm:

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

 20:    References:

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

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

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

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

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

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

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

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

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

102:   PetscFunctionBegin;
103:   if (eps->which==EPS_ALL) {  /* default values in case of spectrum slicing or polynomial filter  */
104:     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
105:     if (isfilt) PetscCall(EPSSetUp_KrylovSchur_Filter(eps));
106:     else PetscCall(EPSSetUp_KrylovSchur_Slice(eps));
107:   } else if (eps->isstructured) {
108:     PetscCall(EPSSetUp_KrylovSchur_BSE(eps));
109:     PetscFunctionReturn(PETSC_SUCCESS);
110:   } else {
111:     PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
112:     PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
113:     if (eps->max_it==PETSC_DETERMINE) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
114:     if (!eps->which) PetscCall(EPSSetWhichEigenpairs_Default(eps));
115:   }
116:   PetscCheck(ctx->lock || eps->mpd>=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");

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

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

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

124:   PetscCall(EPSAllocateSolution(eps,1));
125:   PetscCall(EPS_SetInnerProduct(eps));
126:   if (eps->arbitrary) PetscCall(EPSSetWorkVecs(eps,2));
127:   else if (eps->ishermitian && !eps->ispositive) PetscCall(EPSSetWorkVecs(eps,1));

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

196: static PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
197: {
198:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
199:   SlepcSC         sc;
200:   PetscBool       isfilt;

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

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

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

244:   /* Get the starting Arnoldi vector */
245:   PetscCall(EPSGetStartVector(eps,0,NULL));
246:   l = 0;

248:   /* Restart loop */
249:   while (eps->reason == EPS_CONVERGED_ITERATING) {
250:     eps->its++;

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

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

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

283:     /* Check convergence */
284:     PetscCall(EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k));
285:     PetscCall((*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx));
286:     nconv = k;

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

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

329:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(eps->V,nv,k+l));
330:     eps->nconv = k;
331:     PetscCall(EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv));
332:   }

334:   if (harmonic) PetscCall(PetscFree(g));
335:   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
336:   PetscFunctionReturn(PETSC_SUCCESS);
337: }

339: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
340: {
341:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

343:   PetscFunctionBegin;
344:   if (keep==(PetscReal)PETSC_DEFAULT || keep==(PetscReal)PETSC_DECIDE) ctx->keep = 0.5;
345:   else {
346:     PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument %g must be in the range [0.1,0.9]",(double)keep);
347:     ctx->keep = keep;
348:   }
349:   PetscFunctionReturn(PETSC_SUCCESS);
350: }

352: /*@
353:    EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
354:    method, in particular the proportion of basis vectors that must be kept
355:    after restart.

357:    Logically Collective

359:    Input Parameters:
360: +  eps - the eigenproblem solver context
361: -  keep - the number of vectors to be kept at restart

363:    Options Database Key:
364: .  -eps_krylovschur_restart - Sets the restart parameter

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

369:    Level: advanced

371: .seealso: EPSKrylovSchurGetRestart()
372: @*/
373: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
374: {
375:   PetscFunctionBegin;
378:   PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
379:   PetscFunctionReturn(PETSC_SUCCESS);
380: }

382: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
383: {
384:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

386:   PetscFunctionBegin;
387:   *keep = ctx->keep;
388:   PetscFunctionReturn(PETSC_SUCCESS);
389: }

391: /*@
392:    EPSKrylovSchurGetRestart - Gets the restart parameter used in the
393:    Krylov-Schur method.

395:    Not Collective

397:    Input Parameter:
398: .  eps - the eigenproblem solver context

400:    Output Parameter:
401: .  keep - the restart parameter

403:    Level: advanced

405: .seealso: EPSKrylovSchurSetRestart()
406: @*/
407: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
408: {
409:   PetscFunctionBegin;
411:   PetscAssertPointer(keep,2);
412:   PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

416: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
417: {
418:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

420:   PetscFunctionBegin;
421:   ctx->lock = lock;
422:   PetscFunctionReturn(PETSC_SUCCESS);
423: }

425: /*@
426:    EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
427:    the Krylov-Schur method.

429:    Logically Collective

431:    Input Parameters:
432: +  eps  - the eigenproblem solver context
433: -  lock - true if the locking variant must be selected

435:    Options Database Key:
436: .  -eps_krylovschur_locking - Sets the locking flag

438:    Notes:
439:    The default is to lock converged eigenpairs when the method restarts.
440:    This behaviour can be changed so that all directions are kept in the
441:    working subspace even if already converged to working accuracy (the
442:    non-locking variant).

444:    Level: advanced

446: .seealso: EPSKrylovSchurGetLocking()
447: @*/
448: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
449: {
450:   PetscFunctionBegin;
453:   PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
454:   PetscFunctionReturn(PETSC_SUCCESS);
455: }

457: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
458: {
459:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

461:   PetscFunctionBegin;
462:   *lock = ctx->lock;
463:   PetscFunctionReturn(PETSC_SUCCESS);
464: }

466: /*@
467:    EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
468:    method.

470:    Not Collective

472:    Input Parameter:
473: .  eps - the eigenproblem solver context

475:    Output Parameter:
476: .  lock - the locking flag

478:    Level: advanced

480: .seealso: EPSKrylovSchurSetLocking()
481: @*/
482: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
483: {
484:   PetscFunctionBegin;
486:   PetscAssertPointer(lock,2);
487:   PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
488:   PetscFunctionReturn(PETSC_SUCCESS);
489: }

491: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
492: {
493:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
494:   PetscMPIInt     size;
495:   PetscInt        newnpart;

497:   PetscFunctionBegin;
498:   if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
499:     newnpart = 1;
500:   } else {
501:     PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size));
502:     PetscCheck(npart>0 && npart<=size,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
503:     newnpart = npart;
504:   }
505:   if (ctx->npart!=newnpart) {
506:     if (ctx->npart>1) {
507:       PetscCall(PetscSubcommDestroy(&ctx->subc));
508:       if (ctx->commset) {
509:         PetscCallMPI(MPI_Comm_free(&ctx->commrank));
510:         ctx->commset = PETSC_FALSE;
511:       }
512:     }
513:     PetscCall(EPSDestroy(&ctx->eps));
514:     ctx->npart = newnpart;
515:     eps->state = EPS_STATE_INITIAL;
516:   }
517:   PetscFunctionReturn(PETSC_SUCCESS);
518: }

520: /*@
521:    EPSKrylovSchurSetPartitions - Sets the number of partitions for the
522:    case of doing spectrum slicing for a computational interval with the
523:    communicator split in several sub-communicators.

525:    Logically Collective

527:    Input Parameters:
528: +  eps   - the eigenproblem solver context
529: -  npart - number of partitions

531:    Options Database Key:
532: .  -eps_krylovschur_partitions <npart> - Sets the number of partitions

534:    Notes:
535:    By default, npart=1 so all processes in the communicator participate in
536:    the processing of the whole interval. If npart>1 then the interval is
537:    divided into npart subintervals, each of them being processed by a
538:    subset of processes.

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

543:    Level: advanced

545: .seealso: EPSKrylovSchurSetSubintervals(), EPSSetInterval()
546: @*/
547: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
548: {
549:   PetscFunctionBegin;
552:   PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
553:   PetscFunctionReturn(PETSC_SUCCESS);
554: }

556: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
557: {
558:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

560:   PetscFunctionBegin;
561:   *npart  = ctx->npart;
562:   PetscFunctionReturn(PETSC_SUCCESS);
563: }

565: /*@
566:    EPSKrylovSchurGetPartitions - Gets the number of partitions of the
567:    communicator in case of spectrum slicing.

569:    Not Collective

571:    Input Parameter:
572: .  eps - the eigenproblem solver context

574:    Output Parameter:
575: .  npart - number of partitions

577:    Level: advanced

579: .seealso: EPSKrylovSchurSetPartitions()
580: @*/
581: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
582: {
583:   PetscFunctionBegin;
585:   PetscAssertPointer(npart,2);
586:   PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
587:   PetscFunctionReturn(PETSC_SUCCESS);
588: }

590: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
591: {
592:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

594:   PetscFunctionBegin;
595:   ctx->detect = detect;
596:   eps->state  = EPS_STATE_INITIAL;
597:   PetscFunctionReturn(PETSC_SUCCESS);
598: }

600: /*@
601:    EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
602:    zeros during the factorizations throughout the spectrum slicing computation.

604:    Logically Collective

606:    Input Parameters:
607: +  eps    - the eigenproblem solver context
608: -  detect - check for zeros

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

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

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

622:    Level: advanced

624: .seealso: EPSKrylovSchurSetPartitions(), EPSSetInterval()
625: @*/
626: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
627: {
628:   PetscFunctionBegin;
631:   PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
632:   PetscFunctionReturn(PETSC_SUCCESS);
633: }

635: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
636: {
637:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

639:   PetscFunctionBegin;
640:   *detect = ctx->detect;
641:   PetscFunctionReturn(PETSC_SUCCESS);
642: }

644: /*@
645:    EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
646:    in spectrum slicing.

648:    Not Collective

650:    Input Parameter:
651: .  eps - the eigenproblem solver context

653:    Output Parameter:
654: .  detect - whether zeros detection is enforced during factorizations

656:    Level: advanced

658: .seealso: EPSKrylovSchurSetDetectZeros()
659: @*/
660: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
661: {
662:   PetscFunctionBegin;
664:   PetscAssertPointer(detect,2);
665:   PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
666:   PetscFunctionReturn(PETSC_SUCCESS);
667: }

669: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
670: {
671:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

673:   PetscFunctionBegin;
674:   if (nev != PETSC_CURRENT) {
675:     PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
676:     ctx->nev = nev;
677:   }
678:   if (ncv == PETSC_DETERMINE) {
679:     ctx->ncv = PETSC_DETERMINE;
680:   } else if (ncv != PETSC_CURRENT) {
681:     PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
682:     ctx->ncv = ncv;
683:   }
684:   if (mpd == PETSC_DETERMINE) {
685:     ctx->mpd = PETSC_DETERMINE;
686:   } else if (mpd != PETSC_CURRENT) {
687:     PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
688:     ctx->mpd = mpd;
689:   }
690:   eps->state = EPS_STATE_INITIAL;
691:   PetscFunctionReturn(PETSC_SUCCESS);
692: }

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

699:    Logically Collective

701:    Input Parameters:
702: +  eps - the eigenproblem solver context
703: .  nev - number of eigenvalues to compute
704: .  ncv - the maximum dimension of the subspace to be used by the subsolve
705: -  mpd - the maximum dimension allowed for the projected problem

707:    Options Database Key:
708: +  -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
709: .  -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
710: -  -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension

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

716:    Level: advanced

718: .seealso: EPSKrylovSchurGetDimensions(), EPSSetDimensions(), EPSSetInterval()
719: @*/
720: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
721: {
722:   PetscFunctionBegin;
727:   PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
728:   PetscFunctionReturn(PETSC_SUCCESS);
729: }

731: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
732: {
733:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

735:   PetscFunctionBegin;
736:   if (nev) *nev = ctx->nev;
737:   if (ncv) *ncv = ctx->ncv;
738:   if (mpd) *mpd = ctx->mpd;
739:   PetscFunctionReturn(PETSC_SUCCESS);
740: }

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

746:    Not Collective

748:    Input Parameter:
749: .  eps - the eigenproblem solver context

751:    Output Parameters:
752: +  nev - number of eigenvalues to compute
753: .  ncv - the maximum dimension of the subspace to be used by the subsolve
754: -  mpd - the maximum dimension allowed for the projected problem

756:    Level: advanced

758: .seealso: EPSKrylovSchurSetDimensions()
759: @*/
760: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
761: {
762:   PetscFunctionBegin;
764:   PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
765:   PetscFunctionReturn(PETSC_SUCCESS);
766: }

768: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
769: {
770:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
771:   PetscInt        i;

773:   PetscFunctionBegin;
774:   PetscCheck(subint[0]==eps->inta && subint[ctx->npart]==eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"First and last values must match the endpoints of EPSSetInterval()");
775:   for (i=0;i<ctx->npart;i++) PetscCheck(subint[i]<=subint[i+1],PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Array must contain values in ascending order");
776:   if (ctx->subintervals) PetscCall(PetscFree(ctx->subintervals));
777:   PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
778:   for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
779:   ctx->subintset = PETSC_TRUE;
780:   eps->state = EPS_STATE_INITIAL;
781:   PetscFunctionReturn(PETSC_SUCCESS);
782: }

784: /*@
785:    EPSKrylovSchurSetSubintervals - Sets the points that delimit the
786:    subintervals to be used in spectrum slicing with several partitions.

788:    Logically Collective

790:    Input Parameters:
791: +  eps    - the eigenproblem solver context
792: -  subint - array of real values specifying subintervals

794:    Notes:
795:    This function must be called after EPSKrylovSchurSetPartitions(). For npart
796:    partitions, the argument subint must contain npart+1 real values sorted in
797:    ascending order, subint_0, subint_1, ..., subint_npart, where the first
798:    and last values must coincide with the interval endpoints set with
799:    EPSSetInterval().

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

804:    Level: advanced

806: .seealso: EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubintervals(), EPSSetInterval()
807: @*/
808: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal subint[])
809: {
810:   PetscFunctionBegin;
812:   PetscAssertPointer(subint,2);
813:   PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
814:   PetscFunctionReturn(PETSC_SUCCESS);
815: }

817: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal **subint)
818: {
819:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
820:   PetscInt        i;

822:   PetscFunctionBegin;
823:   if (!ctx->subintset) {
824:     PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
825:     PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
826:   }
827:   PetscCall(PetscMalloc1(ctx->npart+1,subint));
828:   for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
829:   PetscFunctionReturn(PETSC_SUCCESS);
830: }

832: /*@C
833:    EPSKrylovSchurGetSubintervals - Returns the points that delimit the
834:    subintervals used in spectrum slicing with several partitions.

836:    Not Collective

838:    Input Parameter:
839: .  eps    - the eigenproblem solver context

841:    Output Parameter:
842: .  subint - array of real values specifying subintervals

844:    Notes:
845:    If the user passed values with EPSKrylovSchurSetSubintervals(), then the
846:    same values are returned. Otherwise, the values computed internally are
847:    obtained.

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

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

854:    Fortran Notes:
855:    The calling sequence from Fortran is
856: .vb
857:    EPSKrylovSchurGetSubintervals(eps,subint,ierr)
858:    double precision subint(npart+1) output
859: .ve

861:    Level: advanced

863: .seealso: EPSKrylovSchurSetSubintervals(), EPSKrylovSchurGetPartitions(), EPSSetInterval()
864: @*/
865: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal **subint)
866: {
867:   PetscFunctionBegin;
869:   PetscAssertPointer(subint,2);
870:   PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
871:   PetscFunctionReturn(PETSC_SUCCESS);
872: }

874: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
875: {
876:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
877:   PetscInt        i,numsh;
878:   EPS_SR          sr = ctx->sr;

880:   PetscFunctionBegin;
881:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
882:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
883:   switch (eps->state) {
884:   case EPS_STATE_INITIAL:
885:     break;
886:   case EPS_STATE_SETUP:
887:     numsh = ctx->npart+1;
888:     if (n) *n = numsh;
889:     if (shifts) {
890:       PetscCall(PetscMalloc1(numsh,shifts));
891:       (*shifts)[0] = eps->inta;
892:       if (ctx->npart==1) (*shifts)[1] = eps->intb;
893:       else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
894:     }
895:     if (inertias) {
896:       PetscCall(PetscMalloc1(numsh,inertias));
897:       (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
898:       if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
899:       else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
900:     }
901:     break;
902:   case EPS_STATE_SOLVED:
903:   case EPS_STATE_EIGENVECTORS:
904:     numsh = ctx->nshifts;
905:     if (n) *n = numsh;
906:     if (shifts) {
907:       PetscCall(PetscMalloc1(numsh,shifts));
908:       for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
909:     }
910:     if (inertias) {
911:       PetscCall(PetscMalloc1(numsh,inertias));
912:       for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
913:     }
914:     break;
915:   }
916:   PetscFunctionReturn(PETSC_SUCCESS);
917: }

919: /*@C
920:    EPSKrylovSchurGetInertias - Gets the values of the shifts and their
921:    corresponding inertias in case of doing spectrum slicing for a
922:    computational interval.

924:    Not Collective

926:    Input Parameter:
927: .  eps - the eigenproblem solver context

929:    Output Parameters:
930: +  n        - number of shifts, including the endpoints of the interval
931: .  shifts   - the values of the shifts used internally in the solver
932: -  inertias - the values of the inertia in each shift

934:    Notes:
935:    If called after EPSSolve(), all shifts used internally by the solver are
936:    returned (including both endpoints and any intermediate ones). If called
937:    before EPSSolve() and after EPSSetUp() then only the information of the
938:    endpoints of subintervals is available.

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

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

945:    Fortran Notes:
946:    The calling sequence from Fortran is
947: .vb
948:    EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
949:    integer n
950:    double precision shifts(*)
951:    integer inertias(*)
952: .ve
953:    The arrays should be at least of length n. The value of n can be determined
954:    by an initial call
955: .vb
956:    EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL_ARRAY,PETSC_NULL_INTEGER_ARRAY,ierr)
957: .ve

959:    Level: advanced

961: .seealso: EPSSetInterval(), EPSKrylovSchurSetSubintervals()
962: @*/
963: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal *shifts[],PetscInt *inertias[])
964: {
965:   PetscFunctionBegin;
967:   PetscAssertPointer(n,2);
968:   PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
969:   PetscFunctionReturn(PETSC_SUCCESS);
970: }

972: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
973: {
974:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
975:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

977:   PetscFunctionBegin;
978:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
979:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
980:   if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
981:   if (n) *n = sr->numEigs;
982:   if (v) PetscCall(BVCreateVec(sr->V,v));
983:   PetscFunctionReturn(PETSC_SUCCESS);
984: }

986: /*@
987:    EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
988:    doing spectrum slicing for a computational interval with multiple
989:    communicators.

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

993:    Input Parameter:
994: .  eps - the eigenproblem solver context

996:    Output Parameters:
997: +  k - index of the subinterval for the calling process
998: .  n - number of eigenvalues found in the k-th subinterval
999: -  v - a vector owned by processes in the subcommunicator with dimensions
1000:        compatible for locally computed eigenvectors (or NULL)

1002:    Notes:
1003:    This function is only available for spectrum slicing runs.

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

1007:    Level: advanced

1009: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommPairs()
1010: @*/
1011: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1012: {
1013:   PetscFunctionBegin;
1015:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1016:   PetscFunctionReturn(PETSC_SUCCESS);
1017: }

1019: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1020: {
1021:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1022:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

1024:   PetscFunctionBegin;
1025:   EPSCheckSolved(eps,1);
1026:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1027:   PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1028:   if (eig) *eig = sr->eigr[sr->perm[i]];
1029:   if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1030:   PetscFunctionReturn(PETSC_SUCCESS);
1031: }

1033: /*@
1034:    EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
1035:    internally in the subcommunicator to which the calling process belongs.

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

1039:    Input Parameters:
1040: +  eps - the eigenproblem solver context
1041: -  i   - index of the solution

1043:    Output Parameters:
1044: +  eig - the eigenvalue
1045: -  v   - the eigenvector

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

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

1055:    Level: advanced

1057: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo(), EPSKrylovSchurGetSubcommMats()
1058: @*/
1059: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1060: {
1061:   PetscFunctionBegin;
1064:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1065:   PetscFunctionReturn(PETSC_SUCCESS);
1066: }

1068: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1069: {
1070:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

1072:   PetscFunctionBegin;
1073:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1074:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1075:   PetscCall(EPSGetOperators(ctx->eps,A,B));
1076:   PetscFunctionReturn(PETSC_SUCCESS);
1077: }

1079: /*@
1080:    EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1081:    internally in the subcommunicator to which the calling process belongs.

1083:    Collective on the subcommunicator

1085:    Input Parameter:
1086: .  eps - the eigenproblem solver context

1088:    Output Parameters:
1089: +  A  - the matrix associated with the eigensystem
1090: -  B  - the second matrix in the case of generalized eigenproblems

1092:    Notes:
1093:    This is the analog of EPSGetOperators(), but returns the matrices distributed
1094:    differently (in the subcommunicator rather than in the parent communicator).

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

1098:    Level: advanced

1100: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo()
1101: @*/
1102: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1103: {
1104:   PetscFunctionBegin;
1106:   PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1107:   PetscFunctionReturn(PETSC_SUCCESS);
1108: }

1110: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1111: {
1112:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1113:   Mat             A,B=NULL,Ag,Bg=NULL;
1114:   PetscBool       reuse=PETSC_TRUE;

1116:   PetscFunctionBegin;
1117:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1118:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1119:   PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1120:   PetscCall(EPSGetOperators(ctx->eps,&A,&B));

1122:   PetscCall(MatScale(A,a));
1123:   if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1124:   if (B) PetscCall(MatScale(B,b));
1125:   if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1126:   PetscCall(EPSSetOperators(ctx->eps,A,B));

1128:   /* Update stored matrix state */
1129:   subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1130:   PetscCall(MatGetState(A,&subctx->Astate));
1131:   if (B) PetscCall(MatGetState(B,&subctx->Bstate));

1133:   /* Update matrices in the parent communicator if requested by user */
1134:   if (globalup) {
1135:     if (ctx->npart>1) {
1136:       if (!ctx->isrow) {
1137:         PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1138:         reuse = PETSC_FALSE;
1139:       }
1140:       if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1141:       if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1142:       PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1143:       PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1144:       if (B) {
1145:         if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1146:         PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1147:         PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1148:       }
1149:     }
1150:     PetscCall(MatGetState(Ag,&ctx->Astate));
1151:     if (Bg) PetscCall(MatGetState(Bg,&ctx->Bstate));
1152:   }
1153:   PetscCall(EPSSetOperators(eps,Ag,Bg));
1154:   PetscFunctionReturn(PETSC_SUCCESS);
1155: }

1157: /*@
1158:    EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1159:    internally in the subcommunicator to which the calling process belongs.

1161:    Collective

1163:    Input Parameters:
1164: +  eps - the eigenproblem solver context
1165: .  s   - scalar that multiplies the existing A matrix
1166: .  a   - scalar used in the axpy operation on A
1167: .  Au  - matrix used in the axpy operation on A
1168: .  t   - scalar that multiplies the existing B matrix
1169: .  b   - scalar used in the axpy operation on B
1170: .  Bu  - matrix used in the axpy operation on B
1171: .  str - structure flag
1172: -  globalup - flag indicating if global matrices must be updated

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

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

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

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

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

1191:    Level: advanced

1193: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommMats()
1194: @*/
1195: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1196: {
1197:   PetscFunctionBegin;
1207:   PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1208:   PetscFunctionReturn(PETSC_SUCCESS);
1209: }

1211: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1212: {
1213:   EPS_KRYLOVSCHUR  *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1214:   Mat              A,B=NULL,Ar=NULL,Br=NULL;
1215:   PetscMPIInt      rank;
1216:   PetscObjectState Astate,Bstate=0;
1217:   PetscObjectId    Aid,Bid=0;
1218:   STType           sttype;
1219:   PetscInt         nmat;
1220:   const char       *prefix;
1221:   MPI_Comm         child;

1223:   PetscFunctionBegin;
1224:   PetscCall(EPSGetOperators(eps,&A,&B));
1225:   if (ctx->npart==1) {
1226:     if (!ctx->eps) PetscCall(EPSCreate(((PetscObject)eps)->comm,&ctx->eps));
1227:     PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1228:     PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1229:     PetscCall(EPSSetOperators(ctx->eps,A,B));
1230:   } else {
1231:     PetscCall(MatGetState(A,&Astate));
1232:     PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1233:     if (B) {
1234:       PetscCall(MatGetState(B,&Bstate));
1235:       PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1236:     }
1237:     if (!ctx->subc) {
1238:       /* Create context for subcommunicators */
1239:       PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1240:       PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1241:       PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1242:       PetscCall(PetscSubcommGetChild(ctx->subc,&child));

1244:       /* Duplicate matrices */
1245:       PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1246:       ctx->Astate = Astate;
1247:       ctx->Aid = Aid;
1248:       PetscCall(MatPropagateSymmetryOptions(A,Ar));
1249:       if (B) {
1250:         PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1251:         ctx->Bstate = Bstate;
1252:         ctx->Bid = Bid;
1253:         PetscCall(MatPropagateSymmetryOptions(B,Br));
1254:       }
1255:     } else {
1256:       PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1257:       if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1258:         PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1259:         if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
1260:         PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1261:         ctx->Astate = Astate;
1262:         ctx->Aid = Aid;
1263:         PetscCall(MatPropagateSymmetryOptions(A,Ar));
1264:         if (B) {
1265:           PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1266:           ctx->Bstate = Bstate;
1267:           ctx->Bid = Bid;
1268:           PetscCall(MatPropagateSymmetryOptions(B,Br));
1269:         }
1270:         PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1271:         PetscCall(MatDestroy(&Ar));
1272:         PetscCall(MatDestroy(&Br));
1273:       }
1274:     }

1276:     /* Create auxiliary EPS */
1277:     if (!ctx->eps) {
1278:       PetscCall(EPSCreate(child,&ctx->eps));
1279:       PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1280:       PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1281:       PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1282:       PetscCall(MatDestroy(&Ar));
1283:       PetscCall(MatDestroy(&Br));
1284:     }
1285:     /* Create subcommunicator grouping processes with same rank */
1286:     if (!ctx->commset) {
1287:       PetscCallMPI(MPI_Comm_rank(child,&rank));
1288:       PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1289:       ctx->commset = PETSC_TRUE;
1290:     }
1291:   }
1292:   PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1293:   PetscCall(STGetType(eps->st,&sttype));
1294:   PetscCall(STSetType(ctx->eps->st,sttype));

1296:   ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1297:   ctx_local->npart = ctx->npart;
1298:   ctx_local->global = PETSC_FALSE;
1299:   ctx_local->eps = eps;
1300:   ctx_local->subc = ctx->subc;
1301:   ctx_local->commrank = ctx->commrank;
1302:   *childeps = ctx->eps;
1303:   PetscFunctionReturn(PETSC_SUCCESS);
1304: }

1306: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1307: {
1308:   EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1309:   ST              st;
1310:   PetscBool       isfilt;

1312:   PetscFunctionBegin;
1313:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1314:   PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1315:   PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1316:   PetscCall(EPSGetST(ctx->eps,&st));
1317:   PetscCall(STGetOperator(st,NULL));
1318:   PetscCall(STGetKSP(st,ksp));
1319:   PetscFunctionReturn(PETSC_SUCCESS);
1320: }

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

1326:    Collective

1328:    Input Parameter:
1329: .  eps - the eigenproblem solver context

1331:    Output Parameter:
1332: .  ksp - the internal KSP object

1334:    Notes:
1335:    When invoked to compute all eigenvalues in an interval with spectrum
1336:    slicing, EPSKRYLOVSCHUR creates another EPS object internally that is
1337:    used to compute eigenvalues by chunks near selected shifts. This function
1338:    allows access to the KSP object associated to this internal EPS object.

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

1345:    Level: advanced

1347: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions()
1348: @*/
1349: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1350: {
1351:   PetscFunctionBegin;
1353:   PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1354:   PetscFunctionReturn(PETSC_SUCCESS);
1355: }

1357: static PetscErrorCode EPSKrylovSchurSetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType bse)
1358: {
1359:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

1361:   PetscFunctionBegin;
1362:   switch (bse) {
1363:     case EPS_KRYLOVSCHUR_BSE_SHAO:
1364:     case EPS_KRYLOVSCHUR_BSE_GRUNING:
1365:     case EPS_KRYLOVSCHUR_BSE_PROJECTEDBSE:
1366:       if (ctx->bse != bse) {
1367:         ctx->bse = bse;
1368:         eps->state = EPS_STATE_INITIAL;
1369:       }
1370:       break;
1371:     default:
1372:       SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Invalid BSE type");
1373:   }
1374:   PetscFunctionReturn(PETSC_SUCCESS);
1375: }

1377: /*@
1378:    EPSKrylovSchurSetBSEType - Sets the method to be used for BSE structured eigenproblems
1379:    in the Krylov-Schur solver.

1381:    Logically Collective

1383:    Input Parameters:
1384: +  eps - the eigenproblem solver context
1385: -  bse - the BSE method

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

1390:    Level: advanced

1392: .seealso: EPSKrylovSchurGetBSEType(), EPSKrylovSchurBSEType, MatCreateBSE()
1393: @*/
1394: PetscErrorCode EPSKrylovSchurSetBSEType(EPS eps,EPSKrylovSchurBSEType bse)
1395: {
1396:   PetscFunctionBegin;
1399:   PetscTryMethod(eps,"EPSKrylovSchurSetBSEType_C",(EPS,EPSKrylovSchurBSEType),(eps,bse));
1400:   PetscFunctionReturn(PETSC_SUCCESS);
1401: }

1403: static PetscErrorCode EPSKrylovSchurGetBSEType_KrylovSchur(EPS eps,EPSKrylovSchurBSEType *bse)
1404: {
1405:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

1407:   PetscFunctionBegin;
1408:   *bse = ctx->bse;
1409:   PetscFunctionReturn(PETSC_SUCCESS);
1410: }

1412: /*@
1413:    EPSKrylovSchurGetBSEType - Gets the method used for BSE structured eigenproblems
1414:    in the Krylov-Schur solver.

1416:    Not Collective

1418:    Input Parameter:
1419: .  eps - the eigenproblem solver context

1421:    Output Parameter:
1422: .  bse - the BSE method

1424:    Level: advanced

1426: .seealso: EPSKrylovSchurSetBSEType(), EPSKrylovSchurBSEType, MatCreateBSE()
1427: @*/
1428: PetscErrorCode EPSKrylovSchurGetBSEType(EPS eps,EPSKrylovSchurBSEType *bse)
1429: {
1430:   PetscFunctionBegin;
1432:   PetscAssertPointer(bse,2);
1433:   PetscUseMethod(eps,"EPSKrylovSchurGetBSEType_C",(EPS,EPSKrylovSchurBSEType*),(eps,bse));
1434:   PetscFunctionReturn(PETSC_SUCCESS);
1435: }

1437: static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems *PetscOptionsObject)
1438: {
1439:   EPS_KRYLOVSCHUR       *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1440:   PetscBool             flg,lock,b,f1,f2,f3,isfilt;
1441:   PetscReal             keep;
1442:   PetscInt              i,j,k;
1443:   KSP                   ksp;
1444:   EPSKrylovSchurBSEType bse;

1446:   PetscFunctionBegin;
1447:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");

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

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

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

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

1463:     i = 1;
1464:     j = k = PETSC_DECIDE;
1465:     PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1466:     PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1467:     PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1468:     if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));

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

1473:   PetscOptionsHeadEnd();

1475:   /* set options of child KSP in spectrum slicing */
1476:   if (eps->which==EPS_ALL) {
1477:     if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1478:     PetscCall(EPSSetDefaultST(eps));
1479:     PetscCall(STSetFromOptions(eps->st));  /* need to advance this to check ST type */
1480:     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1481:     if (!isfilt) {
1482:       PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1483:       PetscCall(KSPSetFromOptions(ksp));
1484:     }
1485:   }
1486:   PetscFunctionReturn(PETSC_SUCCESS);
1487: }

1489: static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1490: {
1491:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1492:   PetscBool       isascii,isfilt;
1493:   KSP             ksp;
1494:   PetscViewer     sviewer;

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

1529: static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1530: {
1531:   PetscBool      isfilt;

1533:   PetscFunctionBegin;
1534:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1535:   if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1536:   PetscCall(PetscFree(eps->data));
1537:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1538:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1539:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1540:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1541:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1542:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1543:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1544:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1545:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1546:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1547:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1548:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1549:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1550:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1551:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1552:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1553:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1554:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1555:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",NULL));
1556:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",NULL));
1557:   PetscFunctionReturn(PETSC_SUCCESS);
1558: }

1560: static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1561: {
1562:   PetscBool      isfilt;

1564:   PetscFunctionBegin;
1565:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1566:   if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1567:   PetscFunctionReturn(PETSC_SUCCESS);
1568: }

1570: static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1571: {
1572:   PetscFunctionBegin;
1573:   if (eps->which==EPS_ALL) {
1574:     if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1575:   }
1576:   PetscFunctionReturn(PETSC_SUCCESS);
1577: }

1579: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1580: {
1581:   EPS_KRYLOVSCHUR *ctx;

1583:   PetscFunctionBegin;
1584:   PetscCall(PetscNew(&ctx));
1585:   eps->data   = (void*)ctx;
1586:   ctx->lock   = PETSC_TRUE;
1587:   ctx->nev    = 1;
1588:   ctx->ncv    = PETSC_DETERMINE;
1589:   ctx->mpd    = PETSC_DETERMINE;
1590:   ctx->npart  = 1;
1591:   ctx->detect = PETSC_FALSE;
1592:   ctx->global = PETSC_TRUE;

1594:   eps->useds = PETSC_TRUE;

1596:   /* solve and computevectors determined at setup */
1597:   eps->ops->setup          = EPSSetUp_KrylovSchur;
1598:   eps->ops->setupsort      = EPSSetUpSort_KrylovSchur;
1599:   eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1600:   eps->ops->destroy        = EPSDestroy_KrylovSchur;
1601:   eps->ops->reset          = EPSReset_KrylovSchur;
1602:   eps->ops->view           = EPSView_KrylovSchur;
1603:   eps->ops->backtransform  = EPSBackTransform_Default;
1604:   eps->ops->setdefaultst   = EPSSetDefaultST_KrylovSchur;

1606:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1607:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1608:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1609:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1610:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1611:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1612:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1613:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1614:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1615:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1616:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1617:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1618:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1619:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1620:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1621:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1622:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1623:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1624:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetBSEType_C",EPSKrylovSchurSetBSEType_KrylovSchur));
1625:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetBSEType_C",EPSKrylovSchurGetBSEType_KrylovSchur));
1626:   PetscFunctionReturn(PETSC_SUCCESS);
1627: }