Actual source code: krylovschur.c

slepc-3.20.2 2024-03-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_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL*1e-2;  /* use tighter tolerance */
 74:   PetscCall(STFilterSetInterval(eps->st,eps->inta,eps->intb));
 75:   if (!ctx->estimatedrange) {
 76:     PetscCall(STFilterGetRange(eps->st,&rleft,&rright));
 77:     estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
 78:   }
 79:   if (estimaterange) { /* user did not set a range */
 80:     PetscCall(STGetMatrix(eps->st,0,&A));
 81:     PetscCall(MatEstimateSpectralRange_EPS(A,&rleft,&rright));
 82:     PetscCall(PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright));
 83:     PetscCall(STFilterSetRange(eps->st,rleft,rright));
 84:     ctx->estimatedrange = PETSC_TRUE;
 85:   }
 86:   if (eps->ncv==PETSC_DEFAULT && eps->nev==1) eps->nev = 40;  /* user did not provide nev estimation */
 87:   PetscCall(EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd));
 88:   PetscCheck(eps->ncv<=eps->nev+eps->mpd,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nev+mpd");
 89:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

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

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

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

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

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

121:   PetscCall(EPSAllocateSolution(eps,1));
122:   PetscCall(EPS_SetInnerProduct(eps));
123:   if (eps->arbitrary) PetscCall(EPSSetWorkVecs(eps,2));
124:   else if (eps->ishermitian && !eps->ispositive) PetscCall(EPSSetWorkVecs(eps,1));

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

193: static PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
194: {
195:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
196:   SlepcSC         sc;
197:   PetscBool       isfilt;

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

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

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

241:   /* Get the starting Arnoldi vector */
242:   PetscCall(EPSGetStartVector(eps,0,NULL));
243:   l = 0;

245:   /* Restart loop */
246:   while (eps->reason == EPS_CONVERGED_ITERATING) {
247:     eps->its++;

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

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

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

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

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

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

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

331:   if (harmonic) PetscCall(PetscFree(g));
332:   PetscCall(DSTruncate(eps->ds,eps->nconv,PETSC_TRUE));
333:   PetscFunctionReturn(PETSC_SUCCESS);
334: }

336: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
337: {
338:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

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

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

354:    Logically Collective

356:    Input Parameters:
357: +  eps - the eigenproblem solver context
358: -  keep - the number of vectors to be kept at restart

360:    Options Database Key:
361: .  -eps_krylovschur_restart - Sets the restart parameter

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

366:    Level: advanced

368: .seealso: EPSKrylovSchurGetRestart()
369: @*/
370: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
371: {
372:   PetscFunctionBegin;
375:   PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
376:   PetscFunctionReturn(PETSC_SUCCESS);
377: }

379: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
380: {
381:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

383:   PetscFunctionBegin;
384:   *keep = ctx->keep;
385:   PetscFunctionReturn(PETSC_SUCCESS);
386: }

388: /*@
389:    EPSKrylovSchurGetRestart - Gets the restart parameter used in the
390:    Krylov-Schur method.

392:    Not Collective

394:    Input Parameter:
395: .  eps - the eigenproblem solver context

397:    Output Parameter:
398: .  keep - the restart parameter

400:    Level: advanced

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

413: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
414: {
415:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

417:   PetscFunctionBegin;
418:   ctx->lock = lock;
419:   PetscFunctionReturn(PETSC_SUCCESS);
420: }

422: /*@
423:    EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
424:    the Krylov-Schur method.

426:    Logically Collective

428:    Input Parameters:
429: +  eps  - the eigenproblem solver context
430: -  lock - true if the locking variant must be selected

432:    Options Database Key:
433: .  -eps_krylovschur_locking - Sets the locking flag

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

441:    Level: advanced

443: .seealso: EPSKrylovSchurGetLocking()
444: @*/
445: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
446: {
447:   PetscFunctionBegin;
450:   PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
451:   PetscFunctionReturn(PETSC_SUCCESS);
452: }

454: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
455: {
456:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

458:   PetscFunctionBegin;
459:   *lock = ctx->lock;
460:   PetscFunctionReturn(PETSC_SUCCESS);
461: }

463: /*@
464:    EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
465:    method.

467:    Not Collective

469:    Input Parameter:
470: .  eps - the eigenproblem solver context

472:    Output Parameter:
473: .  lock - the locking flag

475:    Level: advanced

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

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

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

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

522:    Logically Collective

524:    Input Parameters:
525: +  eps   - the eigenproblem solver context
526: -  npart - number of partitions

528:    Options Database Key:
529: .  -eps_krylovschur_partitions <npart> - Sets the number of partitions

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

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

540:    Level: advanced

542: .seealso: EPSKrylovSchurSetSubintervals(), EPSSetInterval()
543: @*/
544: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
545: {
546:   PetscFunctionBegin;
549:   PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
550:   PetscFunctionReturn(PETSC_SUCCESS);
551: }

553: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
554: {
555:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

557:   PetscFunctionBegin;
558:   *npart  = ctx->npart;
559:   PetscFunctionReturn(PETSC_SUCCESS);
560: }

562: /*@
563:    EPSKrylovSchurGetPartitions - Gets the number of partitions of the
564:    communicator in case of spectrum slicing.

566:    Not Collective

568:    Input Parameter:
569: .  eps - the eigenproblem solver context

571:    Output Parameter:
572: .  npart - number of partitions

574:    Level: advanced

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

587: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
588: {
589:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

591:   PetscFunctionBegin;
592:   ctx->detect = detect;
593:   eps->state  = EPS_STATE_INITIAL;
594:   PetscFunctionReturn(PETSC_SUCCESS);
595: }

597: /*@
598:    EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
599:    zeros during the factorizations throughout the spectrum slicing computation.

601:    Logically Collective

603:    Input Parameters:
604: +  eps    - the eigenproblem solver context
605: -  detect - check for zeros

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

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

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

619:    Level: advanced

621: .seealso: EPSKrylovSchurSetPartitions(), EPSSetInterval()
622: @*/
623: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
624: {
625:   PetscFunctionBegin;
628:   PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
629:   PetscFunctionReturn(PETSC_SUCCESS);
630: }

632: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
633: {
634:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

636:   PetscFunctionBegin;
637:   *detect = ctx->detect;
638:   PetscFunctionReturn(PETSC_SUCCESS);
639: }

641: /*@
642:    EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
643:    in spectrum slicing.

645:    Not Collective

647:    Input Parameter:
648: .  eps - the eigenproblem solver context

650:    Output Parameter:
651: .  detect - whether zeros detection is enforced during factorizations

653:    Level: advanced

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

666: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
667: {
668:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

670:   PetscFunctionBegin;
671:   PetscCheck(nev>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
672:   ctx->nev = nev;
673:   if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
674:     ctx->ncv = PETSC_DEFAULT;
675:   } else {
676:     PetscCheck(ncv>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
677:     ctx->ncv = ncv;
678:   }
679:   if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
680:     ctx->mpd = PETSC_DEFAULT;
681:   } else {
682:     PetscCheck(mpd>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
683:     ctx->mpd = mpd;
684:   }
685:   eps->state = EPS_STATE_INITIAL;
686:   PetscFunctionReturn(PETSC_SUCCESS);
687: }

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

694:    Logically Collective

696:    Input Parameters:
697: +  eps - the eigenproblem solver context
698: .  nev - number of eigenvalues to compute
699: .  ncv - the maximum dimension of the subspace to be used by the subsolve
700: -  mpd - the maximum dimension allowed for the projected problem

702:    Options Database Key:
703: +  -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
704: .  -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
705: -  -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension

707:    Level: advanced

709: .seealso: EPSKrylovSchurGetDimensions(), EPSSetDimensions(), EPSSetInterval()
710: @*/
711: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
712: {
713:   PetscFunctionBegin;
718:   PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
719:   PetscFunctionReturn(PETSC_SUCCESS);
720: }

722: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
723: {
724:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

726:   PetscFunctionBegin;
727:   if (nev) *nev = ctx->nev;
728:   if (ncv) *ncv = ctx->ncv;
729:   if (mpd) *mpd = ctx->mpd;
730:   PetscFunctionReturn(PETSC_SUCCESS);
731: }

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

737:    Not Collective

739:    Input Parameter:
740: .  eps - the eigenproblem solver context

742:    Output Parameters:
743: +  nev - number of eigenvalues to compute
744: .  ncv - the maximum dimension of the subspace to be used by the subsolve
745: -  mpd - the maximum dimension allowed for the projected problem

747:    Level: advanced

749: .seealso: EPSKrylovSchurSetDimensions()
750: @*/
751: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
752: {
753:   PetscFunctionBegin;
755:   PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
756:   PetscFunctionReturn(PETSC_SUCCESS);
757: }

759: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
760: {
761:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
762:   PetscInt        i;

764:   PetscFunctionBegin;
765:   PetscCheck(subint[0]==eps->inta && subint[ctx->npart]==eps->intb,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"First and last values must match the endpoints of EPSSetInterval()");
766:   for (i=0;i<ctx->npart;i++) PetscCheck(subint[i]<=subint[i+1],PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Array must contain values in ascending order");
767:   if (ctx->subintervals) PetscCall(PetscFree(ctx->subintervals));
768:   PetscCall(PetscMalloc1(ctx->npart+1,&ctx->subintervals));
769:   for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
770:   ctx->subintset = PETSC_TRUE;
771:   eps->state = EPS_STATE_INITIAL;
772:   PetscFunctionReturn(PETSC_SUCCESS);
773: }

775: /*@C
776:    EPSKrylovSchurSetSubintervals - Sets the points that delimit the
777:    subintervals to be used in spectrum slicing with several partitions.

779:    Logically Collective

781:    Input Parameters:
782: +  eps    - the eigenproblem solver context
783: -  subint - array of real values specifying subintervals

785:    Notes:
786:    This function must be called after EPSKrylovSchurSetPartitions(). For npart
787:    partitions, the argument subint must contain npart+1 real values sorted in
788:    ascending order, subint_0, subint_1, ..., subint_npart, where the first
789:    and last values must coincide with the interval endpoints set with
790:    EPSSetInterval().

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

795:    Level: advanced

797: .seealso: EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubintervals(), EPSSetInterval()
798: @*/
799: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal *subint)
800: {
801:   PetscFunctionBegin;
803:   PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
804:   PetscFunctionReturn(PETSC_SUCCESS);
805: }

807: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal **subint)
808: {
809:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
810:   PetscInt        i;

812:   PetscFunctionBegin;
813:   if (!ctx->subintset) {
814:     PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
815:     PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
816:   }
817:   PetscCall(PetscMalloc1(ctx->npart+1,subint));
818:   for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
819:   PetscFunctionReturn(PETSC_SUCCESS);
820: }

822: /*@C
823:    EPSKrylovSchurGetSubintervals - Returns the points that delimit the
824:    subintervals used in spectrum slicing with several partitions.

826:    Not Collective

828:    Input Parameter:
829: .  eps    - the eigenproblem solver context

831:    Output Parameter:
832: .  subint - array of real values specifying subintervals

834:    Notes:
835:    If the user passed values with EPSKrylovSchurSetSubintervals(), then the
836:    same values are returned. Otherwise, the values computed internally are
837:    obtained.

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

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

844:    Fortran Notes:
845:    The calling sequence from Fortran is
846: .vb
847:    EPSKrylovSchurGetSubintervals(eps,subint,ierr)
848:    double precision subint(npart+1) output
849: .ve

851:    Level: advanced

853: .seealso: EPSKrylovSchurSetSubintervals(), EPSKrylovSchurGetPartitions(), EPSSetInterval()
854: @*/
855: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal **subint)
856: {
857:   PetscFunctionBegin;
859:   PetscAssertPointer(subint,2);
860:   PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
861:   PetscFunctionReturn(PETSC_SUCCESS);
862: }

864: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
865: {
866:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
867:   PetscInt        i,numsh;
868:   EPS_SR          sr = ctx->sr;

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

909: /*@C
910:    EPSKrylovSchurGetInertias - Gets the values of the shifts and their
911:    corresponding inertias in case of doing spectrum slicing for a
912:    computational interval.

914:    Not Collective

916:    Input Parameter:
917: .  eps - the eigenproblem solver context

919:    Output Parameters:
920: +  n        - number of shifts, including the endpoints of the interval
921: .  shifts   - the values of the shifts used internally in the solver
922: -  inertias - the values of the inertia in each shift

924:    Notes:
925:    If called after EPSSolve(), all shifts used internally by the solver are
926:    returned (including both endpoints and any intermediate ones). If called
927:    before EPSSolve() and after EPSSetUp() then only the information of the
928:    endpoints of subintervals is available.

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

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

935:    Fortran Notes:
936:    The calling sequence from Fortran is
937: .vb
938:    EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
939:    integer n
940:    double precision shifts(*)
941:    integer inertias(*)
942: .ve
943:    The arrays should be at least of length n. The value of n can be determined
944:    by an initial call
945: .vb
946:    EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
947: .ve

949:    Level: advanced

951: .seealso: EPSSetInterval(), EPSKrylovSchurSetSubintervals()
952: @*/
953: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
954: {
955:   PetscFunctionBegin;
957:   PetscAssertPointer(n,2);
958:   PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
959:   PetscFunctionReturn(PETSC_SUCCESS);
960: }

962: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
963: {
964:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
965:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

967:   PetscFunctionBegin;
968:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
969:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
970:   if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
971:   if (n) *n = sr->numEigs;
972:   if (v) PetscCall(BVCreateVec(sr->V,v));
973:   PetscFunctionReturn(PETSC_SUCCESS);
974: }

976: /*@C
977:    EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
978:    doing spectrum slicing for a computational interval with multiple
979:    communicators.

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

983:    Input Parameter:
984: .  eps - the eigenproblem solver context

986:    Output Parameters:
987: +  k - index of the subinterval for the calling process
988: .  n - number of eigenvalues found in the k-th subinterval
989: -  v - a vector owned by processes in the subcommunicator with dimensions
990:        compatible for locally computed eigenvectors (or NULL)

992:    Notes:
993:    This function is only available for spectrum slicing runs.

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

997:    Level: advanced

999: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommPairs()
1000: @*/
1001: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1002: {
1003:   PetscFunctionBegin;
1005:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1006:   PetscFunctionReturn(PETSC_SUCCESS);
1007: }

1009: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1010: {
1011:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1012:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

1014:   PetscFunctionBegin;
1015:   EPSCheckSolved(eps,1);
1016:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1017:   PetscCheck(i>=0 && i<sr->numEigs,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1018:   if (eig) *eig = sr->eigr[sr->perm[i]];
1019:   if (v) PetscCall(BVCopyVec(sr->V,sr->perm[i],v));
1020:   PetscFunctionReturn(PETSC_SUCCESS);
1021: }

1023: /*@
1024:    EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
1025:    internally in the subcommunicator to which the calling process belongs.

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

1029:    Input Parameters:
1030: +  eps - the eigenproblem solver context
1031: -  i   - index of the solution

1033:    Output Parameters:
1034: +  eig - the eigenvalue
1035: -  v   - the eigenvector

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

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

1045:    Level: advanced

1047: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo(), EPSKrylovSchurGetSubcommMats()
1048: @*/
1049: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1050: {
1051:   PetscFunctionBegin;
1054:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1055:   PetscFunctionReturn(PETSC_SUCCESS);
1056: }

1058: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1059: {
1060:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

1062:   PetscFunctionBegin;
1063:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1064:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1065:   PetscCall(EPSGetOperators(ctx->eps,A,B));
1066:   PetscFunctionReturn(PETSC_SUCCESS);
1067: }

1069: /*@C
1070:    EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1071:    internally in the subcommunicator to which the calling process belongs.

1073:    Collective on the subcommunicator

1075:    Input Parameter:
1076: .  eps - the eigenproblem solver context

1078:    Output Parameters:
1079: +  A  - the matrix associated with the eigensystem
1080: -  B  - the second matrix in the case of generalized eigenproblems

1082:    Notes:
1083:    This is the analog of EPSGetOperators(), but returns the matrices distributed
1084:    differently (in the subcommunicator rather than in the parent communicator).

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

1088:    Level: advanced

1090: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo()
1091: @*/
1092: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1093: {
1094:   PetscFunctionBegin;
1096:   PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1097:   PetscFunctionReturn(PETSC_SUCCESS);
1098: }

1100: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1101: {
1102:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1103:   Mat             A,B=NULL,Ag,Bg=NULL;
1104:   PetscBool       reuse=PETSC_TRUE;

1106:   PetscFunctionBegin;
1107:   PetscCheck(ctx->sr,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1108:   PetscCheck(eps->state,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1109:   PetscCall(EPSGetOperators(eps,&Ag,&Bg));
1110:   PetscCall(EPSGetOperators(ctx->eps,&A,&B));

1112:   PetscCall(MatScale(A,a));
1113:   if (Au) PetscCall(MatAXPY(A,ap,Au,str));
1114:   if (B) PetscCall(MatScale(B,b));
1115:   if (Bu) PetscCall(MatAXPY(B,bp,Bu,str));
1116:   PetscCall(EPSSetOperators(ctx->eps,A,B));

1118:   /* Update stored matrix state */
1119:   subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1120:   PetscCall(PetscObjectStateGet((PetscObject)A,&subctx->Astate));
1121:   if (B) PetscCall(PetscObjectStateGet((PetscObject)B,&subctx->Bstate));

1123:   /* Update matrices in the parent communicator if requested by user */
1124:   if (globalup) {
1125:     if (ctx->npart>1) {
1126:       if (!ctx->isrow) {
1127:         PetscCall(MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol));
1128:         reuse = PETSC_FALSE;
1129:       }
1130:       if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1131:       if (ctx->submata && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submata));
1132:       PetscCall(MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata));
1133:       PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag));
1134:       if (B) {
1135:         if (ctx->submatb && !reuse) PetscCall(MatDestroyMatrices(1,&ctx->submatb));
1136:         PetscCall(MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb));
1137:         PetscCall(MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg));
1138:       }
1139:     }
1140:     PetscCall(PetscObjectStateGet((PetscObject)Ag,&ctx->Astate));
1141:     if (Bg) PetscCall(PetscObjectStateGet((PetscObject)Bg,&ctx->Bstate));
1142:   }
1143:   PetscCall(EPSSetOperators(eps,Ag,Bg));
1144:   PetscFunctionReturn(PETSC_SUCCESS);
1145: }

1147: /*@
1148:    EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1149:    internally in the subcommunicator to which the calling process belongs.

1151:    Collective

1153:    Input Parameters:
1154: +  eps - the eigenproblem solver context
1155: .  s   - scalar that multiplies the existing A matrix
1156: .  a   - scalar used in the axpy operation on A
1157: .  Au  - matrix used in the axpy operation on A
1158: .  t   - scalar that multiplies the existing B matrix
1159: .  b   - scalar used in the axpy operation on B
1160: .  Bu  - matrix used in the axpy operation on B
1161: .  str - structure flag
1162: -  globalup - flag indicating if global matrices must be updated

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

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

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

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

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

1181:    Level: advanced

1183: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommMats()
1184: @*/
1185: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1186: {
1187:   PetscFunctionBegin;
1197:   PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1198:   PetscFunctionReturn(PETSC_SUCCESS);
1199: }

1201: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1202: {
1203:   EPS_KRYLOVSCHUR  *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1204:   Mat              A,B=NULL,Ar=NULL,Br=NULL;
1205:   PetscMPIInt      rank;
1206:   PetscObjectState Astate,Bstate=0;
1207:   PetscObjectId    Aid,Bid=0;
1208:   STType           sttype;
1209:   PetscInt         nmat;
1210:   const char       *prefix;
1211:   MPI_Comm         child;

1213:   PetscFunctionBegin;
1214:   PetscCall(EPSGetOperators(eps,&A,&B));
1215:   if (ctx->npart==1) {
1216:     if (!ctx->eps) PetscCall(EPSCreate(((PetscObject)eps)->comm,&ctx->eps));
1217:     PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1218:     PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1219:     PetscCall(EPSSetOperators(ctx->eps,A,B));
1220:   } else {
1221:     PetscCall(PetscObjectStateGet((PetscObject)A,&Astate));
1222:     PetscCall(PetscObjectGetId((PetscObject)A,&Aid));
1223:     if (B) {
1224:       PetscCall(PetscObjectStateGet((PetscObject)B,&Bstate));
1225:       PetscCall(PetscObjectGetId((PetscObject)B,&Bid));
1226:     }
1227:     if (!ctx->subc) {
1228:       /* Create context for subcommunicators */
1229:       PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc));
1230:       PetscCall(PetscSubcommSetNumber(ctx->subc,ctx->npart));
1231:       PetscCall(PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS));
1232:       PetscCall(PetscSubcommGetChild(ctx->subc,&child));

1234:       /* Duplicate matrices */
1235:       PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1236:       ctx->Astate = Astate;
1237:       ctx->Aid = Aid;
1238:       PetscCall(MatPropagateSymmetryOptions(A,Ar));
1239:       if (B) {
1240:         PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1241:         ctx->Bstate = Bstate;
1242:         ctx->Bid = Bid;
1243:         PetscCall(MatPropagateSymmetryOptions(B,Br));
1244:       }
1245:     } else {
1246:       PetscCall(PetscSubcommGetChild(ctx->subc,&child));
1247:       if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1248:         PetscCall(STGetNumMatrices(ctx->eps->st,&nmat));
1249:         if (nmat) PetscCall(EPSGetOperators(ctx->eps,&Ar,&Br));
1250:         PetscCall(MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar));
1251:         ctx->Astate = Astate;
1252:         ctx->Aid = Aid;
1253:         PetscCall(MatPropagateSymmetryOptions(A,Ar));
1254:         if (B) {
1255:           PetscCall(MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br));
1256:           ctx->Bstate = Bstate;
1257:           ctx->Bid = Bid;
1258:           PetscCall(MatPropagateSymmetryOptions(B,Br));
1259:         }
1260:         PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1261:         PetscCall(MatDestroy(&Ar));
1262:         PetscCall(MatDestroy(&Br));
1263:       }
1264:     }

1266:     /* Create auxiliary EPS */
1267:     if (!ctx->eps) {
1268:       PetscCall(EPSCreate(child,&ctx->eps));
1269:       PetscCall(EPSGetOptionsPrefix(eps,&prefix));
1270:       PetscCall(EPSSetOptionsPrefix(ctx->eps,prefix));
1271:       PetscCall(EPSSetOperators(ctx->eps,Ar,Br));
1272:       PetscCall(MatDestroy(&Ar));
1273:       PetscCall(MatDestroy(&Br));
1274:     }
1275:     /* Create subcommunicator grouping processes with same rank */
1276:     if (!ctx->commset) {
1277:       PetscCallMPI(MPI_Comm_rank(child,&rank));
1278:       PetscCallMPI(MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank));
1279:       ctx->commset = PETSC_TRUE;
1280:     }
1281:   }
1282:   PetscCall(EPSSetType(ctx->eps,((PetscObject)eps)->type_name));
1283:   PetscCall(STGetType(eps->st,&sttype));
1284:   PetscCall(STSetType(ctx->eps->st,sttype));

1286:   ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1287:   ctx_local->npart = ctx->npart;
1288:   ctx_local->global = PETSC_FALSE;
1289:   ctx_local->eps = eps;
1290:   ctx_local->subc = ctx->subc;
1291:   ctx_local->commrank = ctx->commrank;
1292:   *childeps = ctx->eps;
1293:   PetscFunctionReturn(PETSC_SUCCESS);
1294: }

1296: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1297: {
1298:   EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1299:   ST              st;
1300:   PetscBool       isfilt;

1302:   PetscFunctionBegin;
1303:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1304:   PetscCheck(eps->which==EPS_ALL && !isfilt,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations with spectrum slicing");
1305:   PetscCall(EPSKrylovSchurGetChildEPS(eps,&ctx->eps));
1306:   PetscCall(EPSGetST(ctx->eps,&st));
1307:   PetscCall(STGetOperator(st,NULL));
1308:   PetscCall(STGetKSP(st,ksp));
1309:   PetscFunctionReturn(PETSC_SUCCESS);
1310: }

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

1316:    Collective

1318:    Input Parameter:
1319: .  eps - the eigenproblem solver context

1321:    Output Parameter:
1322: .  ksp - the internal KSP object

1324:    Notes:
1325:    When invoked to compute all eigenvalues in an interval with spectrum
1326:    slicing, EPSKRYLOVSCHUR creates another EPS object internally that is
1327:    used to compute eigenvalues by chunks near selected shifts. This function
1328:    allows access to the KSP object associated to this internal EPS object.

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

1335:    Level: advanced

1337: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions()
1338: @*/
1339: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1340: {
1341:   PetscFunctionBegin;
1343:   PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1344:   PetscFunctionReturn(PETSC_SUCCESS);
1345: }

1347: static PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems *PetscOptionsObject)
1348: {
1349:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1350:   PetscBool       flg,lock,b,f1,f2,f3,isfilt;
1351:   PetscReal       keep;
1352:   PetscInt        i,j,k;
1353:   KSP             ksp;

1355:   PetscFunctionBegin;
1356:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");

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

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

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

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

1372:     i = 1;
1373:     j = k = PETSC_DECIDE;
1374:     PetscCall(PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1));
1375:     PetscCall(PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2));
1376:     PetscCall(PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3));
1377:     if (f1 || f2 || f3) PetscCall(EPSKrylovSchurSetDimensions(eps,i,j,k));

1379:   PetscOptionsHeadEnd();

1381:   /* set options of child KSP in spectrum slicing */
1382:   if (eps->which==EPS_ALL) {
1383:     if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
1384:     PetscCall(EPSSetDefaultST(eps));
1385:     PetscCall(STSetFromOptions(eps->st));  /* need to advance this to check ST type */
1386:     PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1387:     if (!isfilt) {
1388:       PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1389:       PetscCall(KSPSetFromOptions(ksp));
1390:     }
1391:   }
1392:   PetscFunctionReturn(PETSC_SUCCESS);
1393: }

1395: static PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1396: {
1397:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1398:   PetscBool       isascii,isfilt;
1399:   KSP             ksp;
1400:   PetscViewer     sviewer;

1402:   PetscFunctionBegin;
1403:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
1404:   if (isascii) {
1405:     PetscCall(PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep)));
1406:     PetscCall(PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-"));
1407:     if (eps->which==EPS_ALL) {
1408:       PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1409:       if (isfilt) PetscCall(PetscViewerASCIIPrintf(viewer,"  using filtering to extract all eigenvalues in an interval\n"));
1410:       else {
1411:         PetscCall(PetscViewerASCIIPrintf(viewer,"  doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd));
1412:         if (ctx->npart>1) {
1413:           PetscCall(PetscViewerASCIIPrintf(viewer,"  multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart));
1414:           if (ctx->detect) PetscCall(PetscViewerASCIIPrintf(viewer,"  detecting zeros when factorizing at subinterval boundaries\n"));
1415:         }
1416:         /* view child KSP */
1417:         PetscCall(EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp));
1418:         PetscCall(PetscViewerASCIIPushTab(viewer));
1419:         if (ctx->npart>1 && ctx->subc) {
1420:           PetscCall(PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer));
1421:           if (!ctx->subc->color) PetscCall(KSPView(ksp,sviewer));
1422:           PetscCall(PetscViewerFlush(sviewer));
1423:           PetscCall(PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer));
1424:           PetscCall(PetscViewerFlush(viewer));
1425:           /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1426:           PetscCall(PetscViewerASCIIPopSynchronized(viewer));
1427:         } else PetscCall(KSPView(ksp,viewer));
1428:         PetscCall(PetscViewerASCIIPopTab(viewer));
1429:       }
1430:     }
1431:   }
1432:   PetscFunctionReturn(PETSC_SUCCESS);
1433: }

1435: static PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1436: {
1437:   PetscBool      isfilt;

1439:   PetscFunctionBegin;
1440:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1441:   if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSDestroy_KrylovSchur_Slice(eps));
1442:   PetscCall(PetscFree(eps->data));
1443:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL));
1444:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL));
1445:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL));
1446:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL));
1447:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL));
1448:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL));
1449:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL));
1450:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL));
1451:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL));
1452:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL));
1453:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL));
1454:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL));
1455:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL));
1456:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL));
1457:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL));
1458:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL));
1459:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL));
1460:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL));
1461:   PetscFunctionReturn(PETSC_SUCCESS);
1462: }

1464: static PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1465: {
1466:   PetscBool      isfilt;

1468:   PetscFunctionBegin;
1469:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt));
1470:   if (eps->which==EPS_ALL && !isfilt) PetscCall(EPSReset_KrylovSchur_Slice(eps));
1471:   PetscFunctionReturn(PETSC_SUCCESS);
1472: }

1474: static PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1475: {
1476:   PetscFunctionBegin;
1477:   if (eps->which==EPS_ALL) {
1478:     if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSINVERT));
1479:   }
1480:   PetscFunctionReturn(PETSC_SUCCESS);
1481: }

1483: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1484: {
1485:   EPS_KRYLOVSCHUR *ctx;

1487:   PetscFunctionBegin;
1488:   PetscCall(PetscNew(&ctx));
1489:   eps->data   = (void*)ctx;
1490:   ctx->lock   = PETSC_TRUE;
1491:   ctx->nev    = 1;
1492:   ctx->ncv    = PETSC_DEFAULT;
1493:   ctx->mpd    = PETSC_DEFAULT;
1494:   ctx->npart  = 1;
1495:   ctx->detect = PETSC_FALSE;
1496:   ctx->global = PETSC_TRUE;

1498:   eps->useds = PETSC_TRUE;

1500:   /* solve and computevectors determined at setup */
1501:   eps->ops->setup          = EPSSetUp_KrylovSchur;
1502:   eps->ops->setupsort      = EPSSetUpSort_KrylovSchur;
1503:   eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1504:   eps->ops->destroy        = EPSDestroy_KrylovSchur;
1505:   eps->ops->reset          = EPSReset_KrylovSchur;
1506:   eps->ops->view           = EPSView_KrylovSchur;
1507:   eps->ops->backtransform  = EPSBackTransform_Default;
1508:   eps->ops->setdefaultst   = EPSSetDefaultST_KrylovSchur;

1510:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur));
1511:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur));
1512:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur));
1513:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur));
1514:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur));
1515:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur));
1516:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur));
1517:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur));
1518:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur));
1519:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur));
1520:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur));
1521:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur));
1522:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur));
1523:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur));
1524:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur));
1525:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur));
1526:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur));
1527:   PetscCall(PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur));
1528:   PetscFunctionReturn(PETSC_SUCCESS);
1529: }