Actual source code: krylovschur.c

slepc-3.18.0 2022-10-01
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:   DSGetLeadingDimension(eps->ds,&ld);
 42:   DSGetDimensions(eps->ds,&n,&l,NULL,NULL);
 43:   for (i=l;i<n;i++) {
 44:     re = eps->eigr[i];
 45:     im = eps->eigi[i];
 46:     STBackTransform(eps->st,1,&re,&im);
 47:     newi = i;
 48:     DSVectors(eps->ds,DS_MAT_X,&newi,NULL);
 49:     DSGetArray(eps->ds,DS_MAT_X,&X);
 50:     Zr = X+i*ld;
 51:     if (newi==i+1) Zi = X+newi*ld;
 52:     else Zi = NULL;
 53:     EPSComputeRitzVector(eps,Zr,Zi,eps->V,xr,xi);
 54:     DSRestoreArray(eps->ds,DS_MAT_X,&X);
 55:     (*eps->arbitrary)(re,im,xr,xi,rr+i,ri+i,eps->arbitraryctx);
 56:   }
 57:   return 0;
 58: }

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

 67:   EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
 68:   EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
 70:   EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
 71:   if (eps->tol==PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL*1e-2;  /* use tighter tolerance */
 72:   STFilterSetInterval(eps->st,eps->inta,eps->intb);
 73:   if (!ctx->estimatedrange) {
 74:     STFilterGetRange(eps->st,&rleft,&rright);
 75:     estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
 76:   }
 77:   if (estimaterange) { /* user did not set a range */
 78:     STGetMatrix(eps->st,0,&A);
 79:     MatEstimateSpectralRange_EPS(A,&rleft,&rright);
 80:     PetscInfo(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright);
 81:     STFilterSetRange(eps->st,rleft,rright);
 82:     ctx->estimatedrange = PETSC_TRUE;
 83:   }
 84:   if (eps->ncv==PETSC_DEFAULT && eps->nev==1) eps->nev = 40;  /* user did not provide nev estimation */
 85:   EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
 87:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
 88:   return 0;
 89: }

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

100:   if (eps->which==EPS_ALL) {  /* default values in case of spectrum slicing or polynomial filter  */
101:     PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
102:     if (isfilt) EPSSetUp_KrylovSchur_Filter(eps);
103:     else EPSSetUp_KrylovSchur_Slice(eps);
104:   } else {
105:     EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
107:     if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
108:     if (!eps->which) EPSSetWhichEigenpairs_Default(eps);
109:   }

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


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

118:   EPSAllocateSolution(eps,1);
119:   EPS_SetInnerProduct(eps);
120:   if (eps->arbitrary) EPSSetWorkVecs(eps,2);
121:   else if (eps->ishermitian && !eps->ispositive) EPSSetWorkVecs(eps,1);

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

190: PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
191: {
192:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
193:   SlepcSC         sc;
194:   PetscBool       isfilt;

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

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

229:   DSGetLeadingDimension(eps->ds,&ld);
230:   harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
231:   hermitian = (eps->ishermitian && !harmonic)?PETSC_TRUE:PETSC_FALSE;
232:   if (harmonic) PetscMalloc1(ld,&g);
233:   if (eps->arbitrary) pj = &j;
234:   else pj = NULL;

236:   /* Get the starting Arnoldi vector */
237:   EPSGetStartVector(eps,0,NULL);
238:   l = 0;

240:   /* Restart loop */
241:   while (eps->reason == EPS_CONVERGED_ITERATING) {
242:     eps->its++;

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

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

265:     /* Solve projected problem */
266:     DSSolve(eps->ds,eps->eigr,eps->eigi);
267:     if (PetscUnlikely(eps->arbitrary)) {
268:       EPSGetArbitraryValues(eps,eps->rr,eps->ri);
269:       j=1;
270:     }
271:     DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj);
272:     DSUpdateExtraRow(eps->ds);
273:     DSSynchronize(eps->ds,eps->eigr,eps->eigi);

275:     /* Check convergence */
276:     EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k);
277:     (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
278:     nconv = k;

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

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

321:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) BVCopyColumn(eps->V,nv,k+l);
322:     eps->nconv = k;
323:     EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv);
324:   }

326:   if (harmonic) PetscFree(g);
327:   DSTruncate(eps->ds,eps->nconv,PETSC_TRUE);
328:   return 0;
329: }

331: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
332: {
333:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

335:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
336:   else {
338:     ctx->keep = keep;
339:   }
340:   return 0;
341: }

343: /*@
344:    EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
345:    method, in particular the proportion of basis vectors that must be kept
346:    after restart.

348:    Logically Collective on eps

350:    Input Parameters:
351: +  eps - the eigenproblem solver context
352: -  keep - the number of vectors to be kept at restart

354:    Options Database Key:
355: .  -eps_krylovschur_restart - Sets the restart parameter

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

360:    Level: advanced

362: .seealso: EPSKrylovSchurGetRestart()
363: @*/
364: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
365: {
368:   PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
369:   return 0;
370: }

372: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
373: {
374:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

376:   *keep = ctx->keep;
377:   return 0;
378: }

380: /*@
381:    EPSKrylovSchurGetRestart - Gets the restart parameter used in the
382:    Krylov-Schur method.

384:    Not Collective

386:    Input Parameter:
387: .  eps - the eigenproblem solver context

389:    Output Parameter:
390: .  keep - the restart parameter

392:    Level: advanced

394: .seealso: EPSKrylovSchurSetRestart()
395: @*/
396: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
397: {
400:   PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
401:   return 0;
402: }

404: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
405: {
406:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

408:   ctx->lock = lock;
409:   return 0;
410: }

412: /*@
413:    EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
414:    the Krylov-Schur method.

416:    Logically Collective on eps

418:    Input Parameters:
419: +  eps  - the eigenproblem solver context
420: -  lock - true if the locking variant must be selected

422:    Options Database Key:
423: .  -eps_krylovschur_locking - Sets the locking flag

425:    Notes:
426:    The default is to lock converged eigenpairs when the method restarts.
427:    This behaviour can be changed so that all directions are kept in the
428:    working subspace even if already converged to working accuracy (the
429:    non-locking variant).

431:    Level: advanced

433: .seealso: EPSKrylovSchurGetLocking()
434: @*/
435: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
436: {
439:   PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
440:   return 0;
441: }

443: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
444: {
445:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

447:   *lock = ctx->lock;
448:   return 0;
449: }

451: /*@
452:    EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
453:    method.

455:    Not Collective

457:    Input Parameter:
458: .  eps - the eigenproblem solver context

460:    Output Parameter:
461: .  lock - the locking flag

463:    Level: advanced

465: .seealso: EPSKrylovSchurSetLocking()
466: @*/
467: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
468: {
471:   PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
472:   return 0;
473: }

475: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
476: {
477:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
478:   PetscMPIInt     size;
479:   PetscInt        newnpart;

481:   if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
482:     newnpart = 1;
483:   } else {
484:     MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
486:     newnpart = npart;
487:   }
488:   if (ctx->npart!=newnpart) {
489:     if (ctx->npart>1) {
490:       PetscSubcommDestroy(&ctx->subc);
491:       if (ctx->commset) {
492:         MPI_Comm_free(&ctx->commrank);
493:         ctx->commset = PETSC_FALSE;
494:       }
495:     }
496:     EPSDestroy(&ctx->eps);
497:     ctx->npart = newnpart;
498:     eps->state = EPS_STATE_INITIAL;
499:   }
500:   return 0;
501: }

503: /*@
504:    EPSKrylovSchurSetPartitions - Sets the number of partitions for the
505:    case of doing spectrum slicing for a computational interval with the
506:    communicator split in several sub-communicators.

508:    Logically Collective on eps

510:    Input Parameters:
511: +  eps   - the eigenproblem solver context
512: -  npart - number of partitions

514:    Options Database Key:
515: .  -eps_krylovschur_partitions <npart> - Sets the number of partitions

517:    Notes:
518:    By default, npart=1 so all processes in the communicator participate in
519:    the processing of the whole interval. If npart>1 then the interval is
520:    divided into npart subintervals, each of them being processed by a
521:    subset of processes.

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

526:    Level: advanced

528: .seealso: EPSKrylovSchurSetSubintervals(), EPSSetInterval()
529: @*/
530: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
531: {
534:   PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
535:   return 0;
536: }

538: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
539: {
540:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

542:   *npart  = ctx->npart;
543:   return 0;
544: }

546: /*@
547:    EPSKrylovSchurGetPartitions - Gets the number of partitions of the
548:    communicator in case of spectrum slicing.

550:    Not Collective

552:    Input Parameter:
553: .  eps - the eigenproblem solver context

555:    Output Parameter:
556: .  npart - number of partitions

558:    Level: advanced

560: .seealso: EPSKrylovSchurSetPartitions()
561: @*/
562: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
563: {
566:   PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
567:   return 0;
568: }

570: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
571: {
572:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

574:   ctx->detect = detect;
575:   eps->state  = EPS_STATE_INITIAL;
576:   return 0;
577: }

579: /*@
580:    EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
581:    zeros during the factorizations throughout the spectrum slicing computation.

583:    Logically Collective on eps

585:    Input Parameters:
586: +  eps    - the eigenproblem solver context
587: -  detect - check for zeros

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

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

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

601:    Level: advanced

603: .seealso: EPSKrylovSchurSetPartitions(), EPSSetInterval()
604: @*/
605: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
606: {
609:   PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
610:   return 0;
611: }

613: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
614: {
615:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

617:   *detect = ctx->detect;
618:   return 0;
619: }

621: /*@
622:    EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
623:    in spectrum slicing.

625:    Not Collective

627:    Input Parameter:
628: .  eps - the eigenproblem solver context

630:    Output Parameter:
631: .  detect - whether zeros detection is enforced during factorizations

633:    Level: advanced

635: .seealso: EPSKrylovSchurSetDetectZeros()
636: @*/
637: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
638: {
641:   PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
642:   return 0;
643: }

645: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
646: {
647:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

650:   ctx->nev = nev;
651:   if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
652:     ctx->ncv = PETSC_DEFAULT;
653:   } else {
655:     ctx->ncv = ncv;
656:   }
657:   if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
658:     ctx->mpd = PETSC_DEFAULT;
659:   } else {
661:     ctx->mpd = mpd;
662:   }
663:   eps->state = EPS_STATE_INITIAL;
664:   return 0;
665: }

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

672:    Logically Collective on eps

674:    Input Parameters:
675: +  eps - the eigenproblem solver context
676: .  nev - number of eigenvalues to compute
677: .  ncv - the maximum dimension of the subspace to be used by the subsolve
678: -  mpd - the maximum dimension allowed for the projected problem

680:    Options Database Key:
681: +  -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
682: .  -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
683: -  -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension

685:    Level: advanced

687: .seealso: EPSKrylovSchurGetDimensions(), EPSSetDimensions(), EPSSetInterval()
688: @*/
689: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
690: {
695:   PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
696:   return 0;
697: }

699: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
700: {
701:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

703:   if (nev) *nev = ctx->nev;
704:   if (ncv) *ncv = ctx->ncv;
705:   if (mpd) *mpd = ctx->mpd;
706:   return 0;
707: }

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

713:    Not Collective

715:    Input Parameter:
716: .  eps - the eigenproblem solver context

718:    Output Parameters:
719: +  nev - number of eigenvalues to compute
720: .  ncv - the maximum dimension of the subspace to be used by the subsolve
721: -  mpd - the maximum dimension allowed for the projected problem

723:    Level: advanced

725: .seealso: EPSKrylovSchurSetDimensions()
726: @*/
727: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
728: {
730:   PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
731:   return 0;
732: }

734: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
735: {
736:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
737:   PetscInt        i;

741:   if (ctx->subintervals) PetscFree(ctx->subintervals);
742:   PetscMalloc1(ctx->npart+1,&ctx->subintervals);
743:   for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
744:   ctx->subintset = PETSC_TRUE;
745:   eps->state = EPS_STATE_INITIAL;
746:   return 0;
747: }

749: /*@C
750:    EPSKrylovSchurSetSubintervals - Sets the points that delimit the
751:    subintervals to be used in spectrum slicing with several partitions.

753:    Logically Collective on eps

755:    Input Parameters:
756: +  eps    - the eigenproblem solver context
757: -  subint - array of real values specifying subintervals

759:    Notes:
760:    This function must be called after EPSKrylovSchurSetPartitions(). For npart
761:    partitions, the argument subint must contain npart+1 real values sorted in
762:    ascending order, subint_0, subint_1, ..., subint_npart, where the first
763:    and last values must coincide with the interval endpoints set with
764:    EPSSetInterval().

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

769:    Level: advanced

771: .seealso: EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubintervals(), EPSSetInterval()
772: @*/
773: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal *subint)
774: {
776:   PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
777:   return 0;
778: }

780: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal **subint)
781: {
782:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
783:   PetscInt        i;

785:   if (!ctx->subintset) {
788:   }
789:   PetscMalloc1(ctx->npart+1,subint);
790:   for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
791:   return 0;
792: }

794: /*@C
795:    EPSKrylovSchurGetSubintervals - Returns the points that delimit the
796:    subintervals used in spectrum slicing with several partitions.

798:    Logically Collective on eps

800:    Input Parameter:
801: .  eps    - the eigenproblem solver context

803:    Output Parameter:
804: .  subint - array of real values specifying subintervals

806:    Notes:
807:    If the user passed values with EPSKrylovSchurSetSubintervals(), then the
808:    same values are returned. Otherwise, the values computed internally are
809:    obtained.

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

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

816:    Fortran Notes:
817:    The calling sequence from Fortran is
818: .vb
819:    EPSKrylovSchurGetSubintervals(eps,subint,ierr)
820:    double precision subint(npart+1) output
821: .ve

823:    Level: advanced

825: .seealso: EPSKrylovSchurSetSubintervals(), EPSKrylovSchurGetPartitions(), EPSSetInterval()
826: @*/
827: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal **subint)
828: {
831:   PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
832:   return 0;
833: }

835: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
836: {
837:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
838:   PetscInt        i,numsh;
839:   EPS_SR          sr = ctx->sr;

843:   switch (eps->state) {
844:   case EPS_STATE_INITIAL:
845:     break;
846:   case EPS_STATE_SETUP:
847:     numsh = ctx->npart+1;
848:     if (n) *n = numsh;
849:     if (shifts) {
850:       PetscMalloc1(numsh,shifts);
851:       (*shifts)[0] = eps->inta;
852:       if (ctx->npart==1) (*shifts)[1] = eps->intb;
853:       else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
854:     }
855:     if (inertias) {
856:       PetscMalloc1(numsh,inertias);
857:       (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
858:       if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
859:       else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
860:     }
861:     break;
862:   case EPS_STATE_SOLVED:
863:   case EPS_STATE_EIGENVECTORS:
864:     numsh = ctx->nshifts;
865:     if (n) *n = numsh;
866:     if (shifts) {
867:       PetscMalloc1(numsh,shifts);
868:       for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
869:     }
870:     if (inertias) {
871:       PetscMalloc1(numsh,inertias);
872:       for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
873:     }
874:     break;
875:   }
876:   return 0;
877: }

879: /*@C
880:    EPSKrylovSchurGetInertias - Gets the values of the shifts and their
881:    corresponding inertias in case of doing spectrum slicing for a
882:    computational interval.

884:    Not Collective

886:    Input Parameter:
887: .  eps - the eigenproblem solver context

889:    Output Parameters:
890: +  n        - number of shifts, including the endpoints of the interval
891: .  shifts   - the values of the shifts used internally in the solver
892: -  inertias - the values of the inertia in each shift

894:    Notes:
895:    If called after EPSSolve(), all shifts used internally by the solver are
896:    returned (including both endpoints and any intermediate ones). If called
897:    before EPSSolve() and after EPSSetUp() then only the information of the
898:    endpoints of subintervals is available.

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

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

905:    Fortran Notes:
906:    The calling sequence from Fortran is
907: .vb
908:    EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
909:    integer n
910:    double precision shifts(*)
911:    integer inertias(*)
912: .ve
913:    The arrays should be at least of length n. The value of n can be determined
914:    by an initial call
915: .vb
916:    EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
917: .ve

919:    Level: advanced

921: .seealso: EPSSetInterval(), EPSKrylovSchurSetSubintervals()
922: @*/
923: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
924: {
927:   PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
928:   return 0;
929: }

931: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
932: {
933:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
934:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

938:   if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
939:   if (n) *n = sr->numEigs;
940:   if (v) BVCreateVec(sr->V,v);
941:   return 0;
942: }

944: /*@C
945:    EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
946:    doing spectrum slicing for a computational interval with multiple
947:    communicators.

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

951:    Input Parameter:
952: .  eps - the eigenproblem solver context

954:    Output Parameters:
955: +  k - index of the subinterval for the calling process
956: .  n - number of eigenvalues found in the k-th subinterval
957: -  v - a vector owned by processes in the subcommunicator with dimensions
958:        compatible for locally computed eigenvectors (or NULL)

960:    Notes:
961:    This function is only available for spectrum slicing runs.

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

965:    Level: advanced

967: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommPairs()
968: @*/
969: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
970: {
972:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
973:   return 0;
974: }

976: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
977: {
978:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
979:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

981:   EPSCheckSolved(eps,1);
984:   if (eig) *eig = sr->eigr[sr->perm[i]];
985:   if (v) BVCopyVec(sr->V,sr->perm[i],v);
986:   return 0;
987: }

989: /*@
990:    EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
991:    internally in the subcommunicator to which the calling process belongs.

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

995:    Input Parameters:
996: +  eps - the eigenproblem solver context
997: -  i   - index of the solution

999:    Output Parameters:
1000: +  eig - the eigenvalue
1001: -  v   - the eigenvector

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

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

1011:    Level: advanced

1013: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo(), EPSKrylovSchurGetSubcommMats()
1014: @*/
1015: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1016: {
1019:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1020:   return 0;
1021: }

1023: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1024: {
1025:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

1029:   EPSGetOperators(ctx->eps,A,B);
1030:   return 0;
1031: }

1033: /*@C
1034:    EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1035:    internally in the subcommunicator to which the calling process belongs.

1037:    Collective on the subcommunicator

1039:    Input Parameter:
1040: .  eps - the eigenproblem solver context

1042:    Output Parameters:
1043: +  A  - the matrix associated with the eigensystem
1044: -  B  - the second matrix in the case of generalized eigenproblems

1046:    Notes:
1047:    This is the analog of EPSGetOperators(), but returns the matrices distributed
1048:    differently (in the subcommunicator rather than in the parent communicator).

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

1052:    Level: advanced

1054: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo()
1055: @*/
1056: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1057: {
1059:   PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1060:   return 0;
1061: }

1063: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1064: {
1065:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1066:   Mat             A,B=NULL,Ag,Bg=NULL;
1067:   PetscBool       reuse=PETSC_TRUE;

1071:   EPSGetOperators(eps,&Ag,&Bg);
1072:   EPSGetOperators(ctx->eps,&A,&B);

1074:   MatScale(A,a);
1075:   if (Au) MatAXPY(A,ap,Au,str);
1076:   if (B) MatScale(B,b);
1077:   if (Bu) MatAXPY(B,bp,Bu,str);
1078:   EPSSetOperators(ctx->eps,A,B);

1080:   /* Update stored matrix state */
1081:   subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1082:   PetscObjectStateGet((PetscObject)A,&subctx->Astate);
1083:   if (B) PetscObjectStateGet((PetscObject)B,&subctx->Bstate);

1085:   /* Update matrices in the parent communicator if requested by user */
1086:   if (globalup) {
1087:     if (ctx->npart>1) {
1088:       if (!ctx->isrow) {
1089:         MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol);
1090:         reuse = PETSC_FALSE;
1091:       }
1092:       if (str==DIFFERENT_NONZERO_PATTERN || str==UNKNOWN_NONZERO_PATTERN) reuse = PETSC_FALSE;
1093:       if (ctx->submata && !reuse) MatDestroyMatrices(1,&ctx->submata);
1094:       MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata);
1095:       MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag);
1096:       if (B) {
1097:         if (ctx->submatb && !reuse) MatDestroyMatrices(1,&ctx->submatb);
1098:         MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb);
1099:         MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg);
1100:       }
1101:     }
1102:     PetscObjectStateGet((PetscObject)Ag,&ctx->Astate);
1103:     if (Bg) PetscObjectStateGet((PetscObject)Bg,&ctx->Bstate);
1104:   }
1105:   EPSSetOperators(eps,Ag,Bg);
1106:   return 0;
1107: }

1109: /*@
1110:    EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1111:    internally in the subcommunicator to which the calling process belongs.

1113:    Collective on eps

1115:    Input Parameters:
1116: +  eps - the eigenproblem solver context
1117: .  s   - scalar that multiplies the existing A matrix
1118: .  a   - scalar used in the axpy operation on A
1119: .  Au  - matrix used in the axpy operation on A
1120: .  t   - scalar that multiplies the existing B matrix
1121: .  b   - scalar used in the axpy operation on B
1122: .  Bu  - matrix used in the axpy operation on B
1123: .  str - structure flag
1124: -  globalup - flag indicating if global matrices must be updated

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

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

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

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

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

1143:    Level: advanced

1145: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommMats()
1146: @*/
1147: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1148: {
1158:   PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1159:   return 0;
1160: }

1162: PetscErrorCode EPSKrylovSchurGetChildEPS(EPS eps,EPS *childeps)
1163: {
1164:   EPS_KRYLOVSCHUR  *ctx=(EPS_KRYLOVSCHUR*)eps->data,*ctx_local;
1165:   Mat              A,B=NULL,Ar=NULL,Br=NULL;
1166:   PetscMPIInt      rank;
1167:   PetscObjectState Astate,Bstate=0;
1168:   PetscObjectId    Aid,Bid=0;
1169:   STType           sttype;
1170:   PetscInt         nmat;
1171:   const char       *prefix;
1172:   MPI_Comm         child;

1174:   EPSGetOperators(eps,&A,&B);
1175:   if (ctx->npart==1) {
1176:     if (!ctx->eps) EPSCreate(((PetscObject)eps)->comm,&ctx->eps);
1177:     EPSGetOptionsPrefix(eps,&prefix);
1178:     EPSSetOptionsPrefix(ctx->eps,prefix);
1179:     EPSSetOperators(ctx->eps,A,B);
1180:   } else {
1181:     PetscObjectStateGet((PetscObject)A,&Astate);
1182:     PetscObjectGetId((PetscObject)A,&Aid);
1183:     if (B) {
1184:       PetscObjectStateGet((PetscObject)B,&Bstate);
1185:       PetscObjectGetId((PetscObject)B,&Bid);
1186:     }
1187:     if (!ctx->subc) {
1188:       /* Create context for subcommunicators */
1189:       PetscSubcommCreate(PetscObjectComm((PetscObject)eps),&ctx->subc);
1190:       PetscSubcommSetNumber(ctx->subc,ctx->npart);
1191:       PetscSubcommSetType(ctx->subc,PETSC_SUBCOMM_CONTIGUOUS);
1192:       PetscSubcommGetChild(ctx->subc,&child);

1194:       /* Duplicate matrices */
1195:       MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar);
1196:       ctx->Astate = Astate;
1197:       ctx->Aid = Aid;
1198:       MatPropagateSymmetryOptions(A,Ar);
1199:       if (B) {
1200:         MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br);
1201:         ctx->Bstate = Bstate;
1202:         ctx->Bid = Bid;
1203:         MatPropagateSymmetryOptions(B,Br);
1204:       }
1205:     } else {
1206:       PetscSubcommGetChild(ctx->subc,&child);
1207:       if (ctx->Astate != Astate || (B && ctx->Bstate != Bstate) || ctx->Aid != Aid || (B && ctx->Bid != Bid)) {
1208:         STGetNumMatrices(ctx->eps->st,&nmat);
1209:         if (nmat) EPSGetOperators(ctx->eps,&Ar,&Br);
1210:         MatCreateRedundantMatrix(A,0,child,MAT_INITIAL_MATRIX,&Ar);
1211:         ctx->Astate = Astate;
1212:         ctx->Aid = Aid;
1213:         MatPropagateSymmetryOptions(A,Ar);
1214:         if (B) {
1215:           MatCreateRedundantMatrix(B,0,child,MAT_INITIAL_MATRIX,&Br);
1216:           ctx->Bstate = Bstate;
1217:           ctx->Bid = Bid;
1218:           MatPropagateSymmetryOptions(B,Br);
1219:         }
1220:         EPSSetOperators(ctx->eps,Ar,Br);
1221:         MatDestroy(&Ar);
1222:         MatDestroy(&Br);
1223:       }
1224:     }

1226:     /* Create auxiliary EPS */
1227:     if (!ctx->eps) {
1228:       EPSCreate(child,&ctx->eps);
1229:       EPSGetOptionsPrefix(eps,&prefix);
1230:       EPSSetOptionsPrefix(ctx->eps,prefix);
1231:       EPSSetOperators(ctx->eps,Ar,Br);
1232:       MatDestroy(&Ar);
1233:       MatDestroy(&Br);
1234:     }
1235:     /* Create subcommunicator grouping processes with same rank */
1236:     if (!ctx->commset) {
1237:       MPI_Comm_rank(child,&rank);
1238:       MPI_Comm_split(((PetscObject)eps)->comm,rank,ctx->subc->color,&ctx->commrank);
1239:       ctx->commset = PETSC_TRUE;
1240:     }
1241:   }
1242:   EPSSetType(ctx->eps,((PetscObject)eps)->type_name);
1243:   STGetType(eps->st,&sttype);
1244:   STSetType(ctx->eps->st,sttype);

1246:   ctx_local = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1247:   ctx_local->npart = ctx->npart;
1248:   ctx_local->global = PETSC_FALSE;
1249:   ctx_local->eps = eps;
1250:   ctx_local->subc = ctx->subc;
1251:   ctx_local->commrank = ctx->commrank;
1252:   *childeps = ctx->eps;
1253:   return 0;
1254: }

1256: static PetscErrorCode EPSKrylovSchurGetKSP_KrylovSchur(EPS eps,KSP *ksp)
1257: {
1258:   EPS_KRYLOVSCHUR *ctx=(EPS_KRYLOVSCHUR*)eps->data;
1259:   ST              st;
1260:   PetscBool       isfilt;

1262:   PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1264:   EPSKrylovSchurGetChildEPS(eps,&ctx->eps);
1265:   EPSGetST(ctx->eps,&st);
1266:   STGetOperator(st,NULL);
1267:   STGetKSP(st,ksp);
1268:   return 0;
1269: }

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

1275:    Collective on eps

1277:    Input Parameter:
1278: .  eps - the eigenproblem solver context

1280:    Output Parameter:
1281: .  ksp - the internal KSP object

1283:    Notes:
1284:    When invoked to compute all eigenvalues in an interval with spectrum
1285:    slicing, EPSKRYLOVSCHUR creates another EPS object internally that is
1286:    used to compute eigenvalues by chunks near selected shifts. This function
1287:    allows access to the KSP object associated to this internal EPS object.

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

1294:    Level: advanced

1296: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions()
1297: @*/
1298: PetscErrorCode EPSKrylovSchurGetKSP(EPS eps,KSP *ksp)
1299: {
1301:   PetscUseMethod(eps,"EPSKrylovSchurGetKSP_C",(EPS,KSP*),(eps,ksp));
1302:   return 0;
1303: }

1305: PetscErrorCode EPSSetFromOptions_KrylovSchur(EPS eps,PetscOptionItems *PetscOptionsObject)
1306: {
1307:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1308:   PetscBool       flg,lock,b,f1,f2,f3,isfilt;
1309:   PetscReal       keep;
1310:   PetscInt        i,j,k;
1311:   KSP             ksp;

1313:   PetscOptionsHeadBegin(PetscOptionsObject,"EPS Krylov-Schur Options");

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

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

1321:     i = ctx->npart;
1322:     PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg);
1323:     if (flg) EPSKrylovSchurSetPartitions(eps,i);

1325:     b = ctx->detect;
1326:     PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg);
1327:     if (flg) EPSKrylovSchurSetDetectZeros(eps,b);

1329:     i = 1;
1330:     j = k = PETSC_DECIDE;
1331:     PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1);
1332:     PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2);
1333:     PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3);
1334:     if (f1 || f2 || f3) EPSKrylovSchurSetDimensions(eps,i,j,k);

1336:   PetscOptionsHeadEnd();

1338:   /* set options of child KSP in spectrum slicing */
1339:   if (eps->which==EPS_ALL) {
1340:     if (!eps->st) EPSGetST(eps,&eps->st);
1341:     EPSSetDefaultST(eps);
1342:     STSetFromOptions(eps->st);  /* need to advance this to check ST type */
1343:     PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1344:     if (!isfilt) {
1345:       EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp);
1346:       KSPSetFromOptions(ksp);
1347:     }
1348:   }
1349:   return 0;
1350: }

1352: PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1353: {
1354:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1355:   PetscBool       isascii,isfilt;
1356:   KSP             ksp;
1357:   PetscViewer     sviewer;

1359:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1360:   if (isascii) {
1361:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
1362:     PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
1363:     if (eps->which==EPS_ALL) {
1364:       PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1365:       if (isfilt) PetscViewerASCIIPrintf(viewer,"  using filtering to extract all eigenvalues in an interval\n");
1366:       else {
1367:         PetscViewerASCIIPrintf(viewer,"  doing spectrum slicing with nev=%" PetscInt_FMT ", ncv=%" PetscInt_FMT ", mpd=%" PetscInt_FMT "\n",ctx->nev,ctx->ncv,ctx->mpd);
1368:         if (ctx->npart>1) {
1369:           PetscViewerASCIIPrintf(viewer,"  multi-communicator spectrum slicing with %" PetscInt_FMT " partitions\n",ctx->npart);
1370:           if (ctx->detect) PetscViewerASCIIPrintf(viewer,"  detecting zeros when factorizing at subinterval boundaries\n");
1371:         }
1372:         /* view child KSP */
1373:         EPSKrylovSchurGetKSP_KrylovSchur(eps,&ksp);
1374:         PetscViewerASCIIPushTab(viewer);
1375:         if (ctx->npart>1 && ctx->subc) {
1376:           PetscViewerGetSubViewer(viewer,ctx->subc->child,&sviewer);
1377:           if (!ctx->subc->color) KSPView(ksp,sviewer);
1378:           PetscViewerFlush(sviewer);
1379:           PetscViewerRestoreSubViewer(viewer,ctx->subc->child,&sviewer);
1380:           PetscViewerFlush(viewer);
1381:           /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
1382:           PetscViewerASCIIPopSynchronized(viewer);
1383:         } else KSPView(ksp,viewer);
1384:         PetscViewerASCIIPopTab(viewer);
1385:       }
1386:     }
1387:   }
1388:   return 0;
1389: }

1391: PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1392: {
1393:   PetscBool      isfilt;

1395:   PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1396:   if (eps->which==EPS_ALL && !isfilt) EPSDestroy_KrylovSchur_Slice(eps);
1397:   PetscFree(eps->data);
1398:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL);
1399:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL);
1400:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL);
1401:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL);
1402:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL);
1403:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL);
1404:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL);
1405:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL);
1406:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL);
1407:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL);
1408:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL);
1409:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL);
1410:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL);
1411:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL);
1412:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL);
1413:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL);
1414:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL);
1415:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",NULL);
1416:   return 0;
1417: }

1419: PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1420: {
1421:   PetscBool      isfilt;

1423:   PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1424:   if (eps->which==EPS_ALL && !isfilt) EPSReset_KrylovSchur_Slice(eps);
1425:   return 0;
1426: }

1428: PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1429: {
1430:   if (eps->which==EPS_ALL) {
1431:     if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSINVERT);
1432:   }
1433:   return 0;
1434: }

1436: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1437: {
1438:   EPS_KRYLOVSCHUR *ctx;

1440:   PetscNew(&ctx);
1441:   eps->data   = (void*)ctx;
1442:   ctx->lock   = PETSC_TRUE;
1443:   ctx->nev    = 1;
1444:   ctx->ncv    = PETSC_DEFAULT;
1445:   ctx->mpd    = PETSC_DEFAULT;
1446:   ctx->npart  = 1;
1447:   ctx->detect = PETSC_FALSE;
1448:   ctx->global = PETSC_TRUE;

1450:   eps->useds = PETSC_TRUE;

1452:   /* solve and computevectors determined at setup */
1453:   eps->ops->setup          = EPSSetUp_KrylovSchur;
1454:   eps->ops->setupsort      = EPSSetUpSort_KrylovSchur;
1455:   eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1456:   eps->ops->destroy        = EPSDestroy_KrylovSchur;
1457:   eps->ops->reset          = EPSReset_KrylovSchur;
1458:   eps->ops->view           = EPSView_KrylovSchur;
1459:   eps->ops->backtransform  = EPSBackTransform_Default;
1460:   eps->ops->setdefaultst   = EPSSetDefaultST_KrylovSchur;

1462:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur);
1463:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur);
1464:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur);
1465:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur);
1466:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur);
1467:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur);
1468:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur);
1469:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur);
1470:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur);
1471:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur);
1472:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur);
1473:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur);
1474:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur);
1475:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur);
1476:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur);
1477:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur);
1478:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur);
1479:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetKSP_C",EPSKrylovSchurGetKSP_KrylovSchur);
1480:   return 0;
1481: }