Actual source code: krylovschur.c

slepc-3.14.0 2020-09-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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: {
 38:   PetscInt       i,newi,ld,n,l;
 39:   Vec            xr=eps->work[0],xi=eps->work[1];
 40:   PetscScalar    re,im,*Zr,*Zi,*X;

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

 62: static PetscErrorCode EstimateRange(Mat A,PetscReal *left,PetscReal *right)
 63: {
 65:   PetscInt       nconv;
 66:   PetscScalar    eig0;
 67:   EPS            eps;

 70:   *left = 0.0; *right = 0.0;
 71:   EPSCreate(PetscObjectComm((PetscObject)A),&eps);
 72:   EPSSetOptionsPrefix(eps,"eps_filter_");
 73:   EPSSetOperators(eps,A,NULL);
 74:   EPSSetProblemType(eps,EPS_HEP);
 75:   EPSSetTolerances(eps,1e-3,50);
 76:   EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);
 77:   EPSSolve(eps);
 78:   EPSGetConverged(eps,&nconv);
 79:   if (nconv>0) {
 80:     EPSGetEigenvalue(eps,0,&eig0,NULL);
 81:   } else eig0 = eps->eigr[0];
 82:   *left = PetscRealPart(eig0);
 83:   EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
 84:   EPSSolve(eps);
 85:   EPSGetConverged(eps,&nconv);
 86:   if (nconv>0) {
 87:     EPSGetEigenvalue(eps,0,&eig0,NULL);
 88:   } else eig0 = eps->eigr[0];
 89:   *right = PetscRealPart(eig0);
 90:   EPSDestroy(&eps);
 91:   return(0);
 92: }

 94: static PetscErrorCode EPSSetUp_KrylovSchur_Filter(EPS eps)
 95: {
 97:   PetscReal      rleft,rright;
 98:   Mat            A;

101:   EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
102:   EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
103:   if (eps->intb >= PETSC_MAX_REAL && eps->inta <= PETSC_MIN_REAL) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"The defined computational interval should have at least one of their sides bounded");
104:   EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
105:   if (eps->tol==PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL*1e-2;  /* use tighter tolerance */
106:   STFilterSetInterval(eps->st,eps->inta,eps->intb);
107:   STGetMatrix(eps->st,0,&A);
108:   STFilterGetRange(eps->st,&rleft,&rright);
109:   if (!rleft && !rright) {
110:     EstimateRange(A,&rleft,&rright);
111:     PetscInfo2(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright);
112:     STFilterSetRange(eps->st,rleft,rright);
113:   }
114:   if (eps->ncv==PETSC_DEFAULT && eps->nev==1) eps->nev = 40;  /* user did not provide nev estimation */
115:   EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
116:   if (eps->ncv>eps->nev+eps->mpd) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must not be larger than nev+mpd");
117:   if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
118:   return(0);
119: }

121: PetscErrorCode EPSSetUp_KrylovSchur(EPS eps)
122: {
123:   PetscErrorCode    ierr;
124:   PetscReal         eta;
125:   PetscBool         isfilt=PETSC_FALSE;
126:   BVOrthogType      otype;
127:   BVOrthogBlockType obtype;
128:   EPS_KRYLOVSCHUR   *ctx = (EPS_KRYLOVSCHUR*)eps->data;
129:   enum { EPS_KS_DEFAULT,EPS_KS_SYMM,EPS_KS_SLICE,EPS_KS_FILTER,EPS_KS_INDEF,EPS_KS_TWOSIDED } variant;

132:   if (eps->which==EPS_ALL) {  /* default values in case of spectrum slicing or polynomial filter  */
133:     PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
134:     if (isfilt) {
135:       EPSSetUp_KrylovSchur_Filter(eps);
136:     } else {
137:       EPSSetUp_KrylovSchur_Slice(eps);
138:     }
139:   } else {
140:     EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
141:     if (eps->ncv>eps->nev+eps->mpd) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must not be larger than nev+mpd");
142:     if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
143:     if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
144:   }
145:   if (!ctx->lock && eps->mpd<eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");

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

149:   if (eps->extraction!=EPS_RITZ && eps->extraction!=EPS_HARMONIC) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");

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

153:   EPSAllocateSolution(eps,1);
154:   EPS_SetInnerProduct(eps);
155:   if (eps->arbitrary) {
156:     EPSSetWorkVecs(eps,2);
157:   } else if (eps->ishermitian && !eps->ispositive){
158:     EPSSetWorkVecs(eps,1);
159:   }

161:   /* dispatch solve method */
162:   if (eps->ishermitian) {
163:     if (eps->which==EPS_ALL) {
164:       if (eps->isgeneralized && !eps->ispositive) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectrum slicing not implemented for indefinite problems");
165:       else variant = isfilt? EPS_KS_FILTER: EPS_KS_SLICE;
166:     } else if (eps->isgeneralized && !eps->ispositive) {
167:       variant = EPS_KS_INDEF;
168:     } else {
169:       switch (eps->extraction) {
170:         case EPS_RITZ:     variant = EPS_KS_SYMM; break;
171:         case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
172:         default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
173:       }
174:     }
175:   } else if (eps->twosided) {
176:     variant = EPS_KS_TWOSIDED;
177:   } else {
178:     switch (eps->extraction) {
179:       case EPS_RITZ:     variant = EPS_KS_DEFAULT; break;
180:       case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
181:       default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
182:     }
183:   }
184:   switch (variant) {
185:     case EPS_KS_DEFAULT:
186:       eps->ops->solve = EPSSolve_KrylovSchur_Default;
187:       eps->ops->computevectors = EPSComputeVectors_Schur;
188:       DSSetType(eps->ds,DSNHEP);
189:       DSAllocate(eps->ds,eps->ncv+1);
190:       break;
191:     case EPS_KS_SYMM:
192:     case EPS_KS_FILTER:
193:       eps->ops->solve = EPSSolve_KrylovSchur_Symm;
194:       eps->ops->computevectors = EPSComputeVectors_Hermitian;
195:       DSSetType(eps->ds,DSHEP);
196:       DSSetCompact(eps->ds,PETSC_TRUE);
197:       DSSetExtraRow(eps->ds,PETSC_TRUE);
198:       DSAllocate(eps->ds,eps->ncv+1);
199:       break;
200:     case EPS_KS_SLICE:
201:       eps->ops->solve = EPSSolve_KrylovSchur_Slice;
202:       eps->ops->computevectors = EPSComputeVectors_Slice;
203:       break;
204:     case EPS_KS_INDEF:
205:       eps->ops->solve = EPSSolve_KrylovSchur_Indefinite;
206:       eps->ops->computevectors = EPSComputeVectors_Indefinite;
207:       DSSetType(eps->ds,DSGHIEP);
208:       DSSetCompact(eps->ds,PETSC_TRUE);
209:       DSAllocate(eps->ds,eps->ncv+1);
210:       /* force reorthogonalization for pseudo-Lanczos */
211:       BVGetOrthogonalization(eps->V,&otype,NULL,&eta,&obtype);
212:       BVSetOrthogonalization(eps->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype);
213:       break;
214:     case EPS_KS_TWOSIDED:
215:       eps->ops->solve = EPSSolve_KrylovSchur_TwoSided;
216:       eps->ops->computevectors = EPSComputeVectors_Schur;
217:       DSSetType(eps->ds,DSNHEP);
218:       DSAllocate(eps->ds,eps->ncv+1);
219:       DSSetType(eps->dsts,DSNHEP);
220:       DSAllocate(eps->dsts,eps->ncv+1);
221:       break;
222:     default: SETERRQ(PetscObjectComm((PetscObject)eps),1,"Unexpected error");
223:   }
224:   return(0);
225: }

227: PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)
228: {
229:   PetscErrorCode  ierr;
230:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
231:   SlepcSC         sc;
232:   PetscBool       isfilt;

235:   EPSSetUpSort_Default(eps);
236:   if (eps->which==EPS_ALL) {
237:     PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
238:     if (isfilt) {
239:       DSGetSlepcSC(eps->ds,&sc);
240:       sc->rg            = NULL;
241:       sc->comparison    = SlepcCompareLargestReal;
242:       sc->comparisonctx = NULL;
243:       sc->map           = NULL;
244:       sc->mapobj        = NULL;
245:     } else {
246:       if (!ctx->global && ctx->sr->numEigs>0) {
247:         DSGetSlepcSC(eps->ds,&sc);
248:         sc->rg            = NULL;
249:         sc->comparison    = SlepcCompareLargestMagnitude;
250:         sc->comparisonctx = NULL;
251:         sc->map           = NULL;
252:         sc->mapobj        = NULL;
253:       }
254:     }
255:   }
256:   return(0);
257: }

259: PetscErrorCode EPSSolve_KrylovSchur_Default(EPS eps)
260: {
261:   PetscErrorCode  ierr;
262:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
263:   PetscInt        i,j,*pj,k,l,nv,ld,nconv;
264:   Mat             U,Op;
265:   PetscScalar     *S,*Q,*g;
266:   PetscReal       beta,gamma=1.0;
267:   PetscBool       breakdown,harmonic;

270:   DSGetLeadingDimension(eps->ds,&ld);
271:   harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
272:   if (harmonic) { PetscMalloc1(ld,&g); }
273:   if (eps->arbitrary) pj = &j;
274:   else pj = NULL;

276:   /* Get the starting Arnoldi vector */
277:   EPSGetStartVector(eps,0,NULL);
278:   l = 0;

280:   /* Restart loop */
281:   while (eps->reason == EPS_CONVERGED_ITERATING) {
282:     eps->its++;

284:     /* Compute an nv-step Arnoldi factorization */
285:     nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
286:     STGetOperator(eps->st,&Op);
287:     DSGetArray(eps->ds,DS_MAT_A,&S);
288:     BVMatArnoldi(eps->V,Op,S,ld,eps->nconv+l,&nv,&beta,&breakdown);
289:     DSRestoreArray(eps->ds,DS_MAT_A,&S);
290:     STRestoreOperator(eps->st,&Op);
291:     DSSetDimensions(eps->ds,nv,0,eps->nconv,eps->nconv+l);
292:     if (l==0) {
293:       DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
294:     } else {
295:       DSSetState(eps->ds,DS_STATE_RAW);
296:     }
297:     BVSetActiveColumns(eps->V,eps->nconv,nv);

299:     /* Compute translation of Krylov decomposition if harmonic extraction used */
300:     if (harmonic) {
301:       DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,g,&gamma);
302:     }

304:     /* Solve projected problem */
305:     DSSolve(eps->ds,eps->eigr,eps->eigi);
306:     if (eps->arbitrary) {
307:       EPSGetArbitraryValues(eps,eps->rr,eps->ri);
308:       j=1;
309:     }
310:     DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj);
311:     DSSynchronize(eps->ds,eps->eigr,eps->eigi);

313:     /* Check convergence */
314:     EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k);
315:     (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
316:     nconv = k;

318:     /* Update l */
319:     if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
320:     else {
321:       l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
322: #if !defined(PETSC_USE_COMPLEX)
323:       DSGetArray(eps->ds,DS_MAT_A,&S);
324:       if (S[k+l+(k+l-1)*ld] != 0.0) {
325:         if (k+l<nv-1) l = l+1;
326:         else l = l-1;
327:       }
328:       DSRestoreArray(eps->ds,DS_MAT_A,&S);
329: #endif
330:     }
331:     if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
332:     if (l) { PetscInfo1(eps,"Preparing to restart keeping l=%D vectors\n",l); }

334:     if (eps->reason == EPS_CONVERGED_ITERATING) {
335:       if (breakdown || k==nv) {
336:         /* Start a new Arnoldi factorization */
337:         PetscInfo2(eps,"Breakdown in Krylov-Schur method (it=%D norm=%g)\n",eps->its,(double)beta);
338:         if (k<eps->nev) {
339:           EPSGetStartVector(eps,k,&breakdown);
340:           if (breakdown) {
341:             eps->reason = EPS_DIVERGED_BREAKDOWN;
342:             PetscInfo(eps,"Unable to generate more start vectors\n");
343:           }
344:         }
345:       } else {
346:         /* Undo translation of Krylov decomposition */
347:         if (harmonic) {
348:           DSSetDimensions(eps->ds,nv,0,k,l);
349:           DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma);
350:           /* gamma u^ = u - U*g~ */
351:           BVSetActiveColumns(eps->V,0,nv);
352:           BVMultColumn(eps->V,-1.0,1.0,nv,g);
353:           BVScaleColumn(eps->V,nv,1.0/gamma);
354:           BVSetActiveColumns(eps->V,eps->nconv,nv);
355:         }
356:         /* Prepare the Rayleigh quotient for restart */
357:         DSGetArray(eps->ds,DS_MAT_A,&S);
358:         DSGetArray(eps->ds,DS_MAT_Q,&Q);
359:         for (i=k;i<k+l;i++) {
360:           S[k+l+i*ld] = Q[nv-1+i*ld]*beta*gamma;
361:         }
362:         DSRestoreArray(eps->ds,DS_MAT_A,&S);
363:         DSRestoreArray(eps->ds,DS_MAT_Q,&Q);
364:       }
365:     }
366:     /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
367:     DSGetMat(eps->ds,DS_MAT_Q,&U);
368:     BVMultInPlace(eps->V,U,eps->nconv,k+l);
369:     MatDestroy(&U);

371:     if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
372:       BVCopyColumn(eps->V,nv,k+l);
373:     }
374:     eps->nconv = k;
375:     EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv);
376:   }

378:   if (harmonic) { PetscFree(g); }
379:   /* truncate Schur decomposition and change the state to raw so that
380:      DSVectors() computes eigenvectors from scratch */
381:   DSSetDimensions(eps->ds,eps->nconv,0,0,0);
382:   DSSetState(eps->ds,DS_STATE_RAW);
383:   return(0);
384: }

386: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)
387: {
388:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

391:   if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
392:   else {
393:     if (keep<0.1 || keep>0.9) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument %g must be in the range [0.1,0.9]",keep);
394:     ctx->keep = keep;
395:   }
396:   return(0);
397: }

399: /*@
400:    EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
401:    method, in particular the proportion of basis vectors that must be kept
402:    after restart.

404:    Logically Collective on eps

406:    Input Parameters:
407: +  eps - the eigenproblem solver context
408: -  keep - the number of vectors to be kept at restart

410:    Options Database Key:
411: .  -eps_krylovschur_restart - Sets the restart parameter

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

416:    Level: advanced

418: .seealso: EPSKrylovSchurGetRestart()
419: @*/
420: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)
421: {

427:   PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
428:   return(0);
429: }

431: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)
432: {
433:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

436:   *keep = ctx->keep;
437:   return(0);
438: }

440: /*@
441:    EPSKrylovSchurGetRestart - Gets the restart parameter used in the
442:    Krylov-Schur method.

444:    Not Collective

446:    Input Parameter:
447: .  eps - the eigenproblem solver context

449:    Output Parameter:
450: .  keep - the restart parameter

452:    Level: advanced

454: .seealso: EPSKrylovSchurSetRestart()
455: @*/
456: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)
457: {

463:   PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
464:   return(0);
465: }

467: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)
468: {
469:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

472:   ctx->lock = lock;
473:   return(0);
474: }

476: /*@
477:    EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
478:    the Krylov-Schur method.

480:    Logically Collective on eps

482:    Input Parameters:
483: +  eps  - the eigenproblem solver context
484: -  lock - true if the locking variant must be selected

486:    Options Database Key:
487: .  -eps_krylovschur_locking - Sets the locking flag

489:    Notes:
490:    The default is to lock converged eigenpairs when the method restarts.
491:    This behaviour can be changed so that all directions are kept in the
492:    working subspace even if already converged to working accuracy (the
493:    non-locking variant).

495:    Level: advanced

497: .seealso: EPSKrylovSchurGetLocking()
498: @*/
499: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)
500: {

506:   PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
507:   return(0);
508: }

510: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)
511: {
512:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

515:   *lock = ctx->lock;
516:   return(0);
517: }

519: /*@
520:    EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
521:    method.

523:    Not Collective

525:    Input Parameter:
526: .  eps - the eigenproblem solver context

528:    Output Parameter:
529: .  lock - the locking flag

531:    Level: advanced

533: .seealso: EPSKrylovSchurSetLocking()
534: @*/
535: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)
536: {

542:   PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
543:   return(0);
544: }

546: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)
547: {
548:   PetscErrorCode  ierr;
549:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
550:   PetscMPIInt     size;

553:   if (ctx->npart!=npart) {
554:     if (ctx->commset) { PetscSubcommDestroy(&ctx->subc); }
555:     EPSDestroy(&ctx->eps);
556:   }
557:   if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
558:     ctx->npart = 1;
559:   } else {
560:     MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
561:     if (npart<1 || npart>size) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
562:     ctx->npart = npart;
563:   }
564:   eps->state = EPS_STATE_INITIAL;
565:   return(0);
566: }

568: /*@
569:    EPSKrylovSchurSetPartitions - Sets the number of partitions for the
570:    case of doing spectrum slicing for a computational interval with the
571:    communicator split in several sub-communicators.

573:    Logically Collective on eps

575:    Input Parameters:
576: +  eps   - the eigenproblem solver context
577: -  npart - number of partitions

579:    Options Database Key:
580: .  -eps_krylovschur_partitions <npart> - Sets the number of partitions

582:    Notes:
583:    By default, npart=1 so all processes in the communicator participate in
584:    the processing of the whole interval. If npart>1 then the interval is
585:    divided into npart subintervals, each of them being processed by a
586:    subset of processes.

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

591:    Level: advanced

593: .seealso: EPSKrylovSchurSetSubintervals(), EPSSetInterval()
594: @*/
595: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)
596: {

602:   PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
603:   return(0);
604: }

606: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)
607: {
608:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

611:   *npart  = ctx->npart;
612:   return(0);
613: }

615: /*@
616:    EPSKrylovSchurGetPartitions - Gets the number of partitions of the
617:    communicator in case of spectrum slicing.

619:    Not Collective

621:    Input Parameter:
622: .  eps - the eigenproblem solver context

624:    Output Parameter:
625: .  npart - number of partitions

627:    Level: advanced

629: .seealso: EPSKrylovSchurSetPartitions()
630: @*/
631: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)
632: {

638:   PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
639:   return(0);
640: }

642: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)
643: {
644:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

647:   ctx->detect = detect;
648:   eps->state  = EPS_STATE_INITIAL;
649:   return(0);
650: }

652: /*@
653:    EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
654:    zeros during the factorizations throughout the spectrum slicing computation.

656:    Logically Collective on eps

658:    Input Parameters:
659: +  eps    - the eigenproblem solver context
660: -  detect - check for zeros

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

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

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

674:    Level: advanced

676: .seealso: EPSKrylovSchurSetPartitions(), EPSSetInterval()
677: @*/
678: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)
679: {

685:   PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
686:   return(0);
687: }

689: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)
690: {
691:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

694:   *detect = ctx->detect;
695:   return(0);
696: }

698: /*@
699:    EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
700:    in spectrum slicing.

702:    Not Collective

704:    Input Parameter:
705: .  eps - the eigenproblem solver context

707:    Output Parameter:
708: .  detect - whether zeros detection is enforced during factorizations

710:    Level: advanced

712: .seealso: EPSKrylovSchurSetDetectZeros()
713: @*/
714: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)
715: {

721:   PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
722:   return(0);
723: }

725: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
726: {
727:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

730:   if (nev<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
731:   ctx->nev = nev;
732:   if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
733:     ctx->ncv = PETSC_DEFAULT;
734:   } else {
735:     if (ncv<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
736:     ctx->ncv = ncv;
737:   }
738:   if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
739:     ctx->mpd = PETSC_DEFAULT;
740:   } else {
741:     if (mpd<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
742:     ctx->mpd = mpd;
743:   }
744:   eps->state = EPS_STATE_INITIAL;
745:   return(0);
746: }

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

753:    Logically Collective on eps

755:    Input Parameters:
756: +  eps - the eigenproblem solver context
757: .  nev - number of eigenvalues to compute
758: .  ncv - the maximum dimension of the subspace to be used by the subsolve
759: -  mpd - the maximum dimension allowed for the projected problem

761:    Options Database Key:
762: +  -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
763: .  -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
764: -  -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension

766:    Level: advanced

768: .seealso: EPSKrylovSchurGetDimensions(), EPSSetDimensions(), EPSSetInterval()
769: @*/
770: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)
771: {

779:   PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
780:   return(0);
781: }

783: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
784: {
785:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

788:   if (nev) *nev = ctx->nev;
789:   if (ncv) *ncv = ctx->ncv;
790:   if (mpd) *mpd = ctx->mpd;
791:   return(0);
792: }

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

798:    Not Collective

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

803:    Output Parameters:
804: +  nev - number of eigenvalues to compute
805: .  ncv - the maximum dimension of the subspace to be used by the subsolve
806: -  mpd - the maximum dimension allowed for the projected problem

808:    Level: advanced

810: .seealso: EPSKrylovSchurSetDimensions()
811: @*/
812: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
813: {

818:   PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
819:   return(0);
820: }

822: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)
823: {
824:   PetscErrorCode  ierr;
825:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
826:   PetscInt        i;

829:   if (subint[0]!=eps->inta || subint[ctx->npart]!=eps->intb) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"First and last values must match the endpoints of EPSSetInterval()");
830:   for (i=0;i<ctx->npart;i++) if (subint[i]>subint[i+1]) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Array must contain values in ascending order");
831:   if (ctx->subintervals) { PetscFree(ctx->subintervals); }
832:   PetscMalloc1(ctx->npart+1,&ctx->subintervals);
833:   for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
834:   ctx->subintset = PETSC_TRUE;
835:   eps->state = EPS_STATE_INITIAL;
836:   return(0);
837: }

839: /*@C
840:    EPSKrylovSchurSetSubintervals - Sets the points that delimit the
841:    subintervals to be used in spectrum slicing with several partitions.

843:    Logically Collective on eps

845:    Input Parameters:
846: +  eps    - the eigenproblem solver context
847: -  subint - array of real values specifying subintervals

849:    Notes:
850:    This function must be called after EPSKrylovSchurSetPartitions(). For npart
851:    partitions, the argument subint must contain npart+1 real values sorted in
852:    ascending order: subint_0, subint_1, ..., subint_npart, where the first
853:    and last values must coincide with the interval endpoints set with
854:    EPSSetInterval().

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

859:    Level: advanced

861: .seealso: EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubintervals(), EPSSetInterval()
862: @*/
863: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal *subint)
864: {

869:   PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
870:   return(0);
871: }

873: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal **subint)
874: {
875:   PetscErrorCode  ierr;
876:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
877:   PetscInt        i;

880:   if (!ctx->subintset) {
881:     if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
882:     if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
883:   }
884:   PetscMalloc1(ctx->npart+1,subint);
885:   for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
886:   return(0);
887: }

889: /*@C
890:    EPSKrylovSchurGetSubintervals - Returns the points that delimit the
891:    subintervals used in spectrum slicing with several partitions.

893:    Logically Collective on eps

895:    Input Parameter:
896: .  eps    - the eigenproblem solver context

898:    Output Parameter:
899: .  subint - array of real values specifying subintervals

901:    Notes:
902:    If the user passed values with EPSKrylovSchurSetSubintervals(), then the
903:    same values are returned. Otherwise, the values computed internally are
904:    obtained.

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

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

911:    Fortran Notes:
912:    The calling sequence from Fortran is
913: .vb
914:    EPSKrylovSchurGetSubintervals(eps,subint,ierr)
915:    double precision subint(npart+1) output
916: .ve

918:    Level: advanced

920: .seealso: EPSKrylovSchurSetSubintervals(), EPSKrylovSchurGetPartitions(), EPSSetInterval()
921: @*/
922: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal **subint)
923: {

929:   PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
930:   return(0);
931: }

933: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
934: {
935:   PetscErrorCode  ierr;
936:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
937:   PetscInt        i,numsh;
938:   EPS_SR          sr = ctx->sr;

941:   if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
942:   if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
943:   switch (eps->state) {
944:   case EPS_STATE_INITIAL:
945:     break;
946:   case EPS_STATE_SETUP:
947:     numsh = ctx->npart+1;
948:     if (n) *n = numsh;
949:     if (shifts) {
950:       PetscMalloc1(numsh,shifts);
951:       (*shifts)[0] = eps->inta;
952:       if (ctx->npart==1) (*shifts)[1] = eps->intb;
953:       else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
954:     }
955:     if (inertias) {
956:       PetscMalloc1(numsh,inertias);
957:       (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
958:       if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
959:       else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
960:     }
961:     break;
962:   case EPS_STATE_SOLVED:
963:   case EPS_STATE_EIGENVECTORS:
964:     numsh = ctx->nshifts;
965:     if (n) *n = numsh;
966:     if (shifts) {
967:       PetscMalloc1(numsh,shifts);
968:       for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
969:     }
970:     if (inertias) {
971:       PetscMalloc1(numsh,inertias);
972:       for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
973:     }
974:     break;
975:   }
976:   return(0);
977: }

979: /*@C
980:    EPSKrylovSchurGetInertias - Gets the values of the shifts and their
981:    corresponding inertias in case of doing spectrum slicing for a
982:    computational interval.

984:    Not Collective

986:    Input Parameter:
987: .  eps - the eigenproblem solver context

989:    Output Parameters:
990: +  n        - number of shifts, including the endpoints of the interval
991: .  shifts   - the values of the shifts used internally in the solver
992: -  inertias - the values of the inertia in each shift

994:    Notes:
995:    If called after EPSSolve(), all shifts used internally by the solver are
996:    returned (including both endpoints and any intermediate ones). If called
997:    before EPSSolve() and after EPSSetUp() then only the information of the
998:    endpoints of subintervals is available.

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

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

1005:    Fortran Notes:
1006:    The calling sequence from Fortran is
1007: .vb
1008:    EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
1009:    integer n
1010:    double precision shifts(*)
1011:    integer inertias(*)
1012: .ve
1013:    The arrays should be at least of length n. The value of n can be determined
1014:    by an initial call
1015: .vb
1016:    EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
1017: .ve

1019:    Level: advanced

1021: .seealso: EPSSetInterval(), EPSKrylovSchurSetSubintervals()
1022: @*/
1023: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)
1024: {

1030:   PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
1031:   return(0);
1032: }

1034: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1035: {
1036:   PetscErrorCode  ierr;
1037:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1038:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

1041:   if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1042:   if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1043:   if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
1044:   if (n) *n = sr->numEigs;
1045:   if (v) {
1046:     BVCreateVec(sr->V,v);
1047:   }
1048:   return(0);
1049: }

1051: /*@C
1052:    EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
1053:    doing spectrum slicing for a computational interval with multiple
1054:    communicators.

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

1058:    Input Parameter:
1059: .  eps - the eigenproblem solver context

1061:    Output Parameters:
1062: +  k - index of the subinterval for the calling process
1063: .  n - number of eigenvalues found in the k-th subinterval
1064: -  v - a vector owned by processes in the subcommunicator with dimensions
1065:        compatible for locally computed eigenvectors (or NULL)

1067:    Notes:
1068:    This function is only available for spectrum slicing runs.

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

1072:    Level: advanced

1074: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommPairs()
1075: @*/
1076: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)
1077: {

1082:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1083:   return(0);
1084: }

1086: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1087: {
1088:   PetscErrorCode  ierr;
1089:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1090:   EPS_SR          sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;

1093:   EPSCheckSolved(eps,1);
1094:   if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1095:   if (i<0 || i>=sr->numEigs) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1096:   if (eig) *eig = sr->eigr[sr->perm[i]];
1097:   if (v) { BVCopyVec(sr->V,sr->perm[i],v); }
1098:   return(0);
1099: }

1101: /*@
1102:    EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
1103:    internally in the subcommunicator to which the calling process belongs.

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

1107:    Input Parameter:
1108: +  eps - the eigenproblem solver context
1109: -  i   - index of the solution

1111:    Output Parameters:
1112: +  eig - the eigenvalue
1113: -  v   - the eigenvector

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

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

1123:    Level: advanced

1125: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo(), EPSKrylovSchurGetSubcommMats()
1126: @*/
1127: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)
1128: {

1134:   PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1135:   return(0);
1136: }

1138: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)
1139: {
1140:   PetscErrorCode  ierr;
1141:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;

1144:   if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1145:   if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1146:   EPSGetOperators(ctx->eps,A,B);
1147:   return(0);
1148: }

1150: /*@C
1151:    EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1152:    internally in the subcommunicator to which the calling process belongs.

1154:    Collective on the subcommunicator

1156:    Input Parameter:
1157: .  eps - the eigenproblem solver context

1159:    Output Parameters:
1160: +  A  - the matrix associated with the eigensystem
1161: -  B  - the second matrix in the case of generalized eigenproblems

1163:    Notes:
1164:    This is the analog of EPSGetOperators(), but returns the matrices distributed
1165:    differently (in the subcommunicator rather than in the parent communicator).

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

1169:    Level: advanced

1171: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo()
1172: @*/
1173: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)
1174: {

1179:   PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1180:   return(0);
1181: }

1183: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)
1184: {
1185:   PetscErrorCode  ierr;
1186:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1187:   Mat             A,B=NULL,Ag,Bg=NULL;
1188:   PetscBool       reuse=PETSC_TRUE;

1191:   if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1192:   if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1193:   EPSGetOperators(eps,&Ag,&Bg);
1194:   EPSGetOperators(ctx->eps,&A,&B);

1196:   MatScale(A,a);
1197:   if (Au) {
1198:     MatAXPY(A,ap,Au,str);
1199:   }
1200:   if (B) MatScale(B,b);
1201:   if (Bu) {
1202:     MatAXPY(B,bp,Bu,str);
1203:   }
1204:   EPSSetOperators(ctx->eps,A,B);

1206:   /* Update stored matrix state */
1207:   subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1208:   PetscObjectStateGet((PetscObject)A,&subctx->Astate);
1209:   if (B) { PetscObjectStateGet((PetscObject)B,&subctx->Bstate); }

1211:   /* Update matrices in the parent communicator if requested by user */
1212:   if (globalup) {
1213:     if (ctx->npart>1) {
1214:       if (!ctx->isrow) {
1215:         MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol);
1216:         reuse = PETSC_FALSE;
1217:       }
1218:       if (str==DIFFERENT_NONZERO_PATTERN) reuse = PETSC_FALSE;
1219:       if (ctx->submata && !reuse) {
1220:         MatDestroyMatrices(1,&ctx->submata);
1221:       }
1222:       MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata);
1223:       MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag);
1224:       if (B) {
1225:         if (ctx->submatb && !reuse) {
1226:           MatDestroyMatrices(1,&ctx->submatb);
1227:         }
1228:         MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb);
1229:         MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg);
1230:       }
1231:     }
1232:     PetscObjectStateGet((PetscObject)Ag,&ctx->Astate);
1233:     if (Bg) { PetscObjectStateGet((PetscObject)Bg,&ctx->Bstate); }
1234:   }
1235:   EPSSetOperators(eps,Ag,Bg);
1236:   return(0);
1237: }

1239: /*@
1240:    EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1241:    internally in the subcommunicator to which the calling process belongs.

1243:    Collective on eps

1245:    Input Parameters:
1246: +  eps - the eigenproblem solver context
1247: .  s   - scalar that multiplies the existing A matrix
1248: .  a   - scalar used in the axpy operation on A
1249: .  Au  - matrix used in the axpy operation on A
1250: .  t   - scalar that multiplies the existing B matrix
1251: .  b   - scalar used in the axpy operation on B
1252: .  Bu  - matrix used in the axpy operation on B
1253: .  str - structure flag
1254: -  globalup - flag indicating if global matrices must be updated

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

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

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

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

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

1273:    Level: advanced

1275: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommMats()
1276: @*/
1277: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)
1278: {

1291:   PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1292:   return(0);
1293: }

1295: PetscErrorCode EPSSetFromOptions_KrylovSchur(PetscOptionItems *PetscOptionsObject,EPS eps)
1296: {
1297:   PetscErrorCode  ierr;
1298:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1299:   PetscBool       flg,lock,b,f1,f2,f3;
1300:   PetscReal       keep;
1301:   PetscInt        i,j,k;

1304:   PetscOptionsHead(PetscOptionsObject,"EPS Krylov-Schur Options");

1306:     PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg);
1307:     if (flg) { EPSKrylovSchurSetRestart(eps,keep); }

1309:     PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg);
1310:     if (flg) { EPSKrylovSchurSetLocking(eps,lock); }

1312:     i = ctx->npart;
1313:     PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg);
1314:     if (flg) { EPSKrylovSchurSetPartitions(eps,i); }

1316:     b = ctx->detect;
1317:     PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg);
1318:     if (flg) { EPSKrylovSchurSetDetectZeros(eps,b); }

1320:     i = 1;
1321:     j = k = PETSC_DECIDE;
1322:     PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1);
1323:     PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2);
1324:     PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3);
1325:     if (f1 || f2 || f3) { EPSKrylovSchurSetDimensions(eps,i,j,k); }

1327:   PetscOptionsTail();
1328:   return(0);
1329: }

1331: PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)
1332: {
1333:   PetscErrorCode  ierr;
1334:   EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1335:   PetscBool       isascii,isfilt;

1338:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1339:   if (isascii) {
1340:     PetscViewerASCIIPrintf(viewer,"  %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
1341:     PetscViewerASCIIPrintf(viewer,"  using the %slocking variant\n",ctx->lock?"":"non-");
1342:     if (eps->which==EPS_ALL) {
1343:       PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1344:       if (isfilt) {
1345:         PetscViewerASCIIPrintf(viewer,"  using filtering to extract all eigenvalues in an interval\n");
1346:       } else {
1347:         PetscViewerASCIIPrintf(viewer,"  doing spectrum slicing with nev=%D, ncv=%D, mpd=%D\n",ctx->nev,ctx->ncv,ctx->mpd);
1348:         if (ctx->npart>1) {
1349:           PetscViewerASCIIPrintf(viewer,"  multi-communicator spectrum slicing with %D partitions\n",ctx->npart);
1350:           if (ctx->detect) { PetscViewerASCIIPrintf(viewer,"  detecting zeros when factorizing at subinterval boundaries\n"); }
1351:         }
1352:       }
1353:     }
1354:   }
1355:   return(0);
1356: }

1358: PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)
1359: {

1363:   PetscFree(eps->data);
1364:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL);
1365:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL);
1366:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL);
1367:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL);
1368:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL);
1369:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL);
1370:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL);
1371:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL);
1372:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL);
1373:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL);
1374:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL);
1375:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL);
1376:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL);
1377:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL);
1378:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL);
1379:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL);
1380:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL);
1381:   return(0);
1382: }

1384: PetscErrorCode EPSReset_KrylovSchur(EPS eps)
1385: {
1387:   PetscBool      isfilt;

1390:   PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1391:   if (eps->which==EPS_ALL && !isfilt) {
1392:     EPSReset_KrylovSchur_Slice(eps);
1393:   }
1394:   return(0);
1395: }

1397: PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)
1398: {

1402:   if (eps->which==EPS_ALL) {
1403:     if (!((PetscObject)eps->st)->type_name) {
1404:       STSetType(eps->st,STSINVERT);
1405:     }
1406:   }
1407:   return(0);
1408: }

1410: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)
1411: {
1412:   EPS_KRYLOVSCHUR *ctx;
1413:   PetscErrorCode  ierr;

1416:   PetscNewLog(eps,&ctx);
1417:   eps->data   = (void*)ctx;
1418:   ctx->lock   = PETSC_TRUE;
1419:   ctx->nev    = 1;
1420:   ctx->ncv    = PETSC_DEFAULT;
1421:   ctx->mpd    = PETSC_DEFAULT;
1422:   ctx->npart  = 1;
1423:   ctx->detect = PETSC_FALSE;
1424:   ctx->global = PETSC_TRUE;

1426:   eps->useds = PETSC_TRUE;

1428:   /* solve and computevectors determined at setup */
1429:   eps->ops->setup          = EPSSetUp_KrylovSchur;
1430:   eps->ops->setupsort      = EPSSetUpSort_KrylovSchur;
1431:   eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1432:   eps->ops->destroy        = EPSDestroy_KrylovSchur;
1433:   eps->ops->reset          = EPSReset_KrylovSchur;
1434:   eps->ops->view           = EPSView_KrylovSchur;
1435:   eps->ops->backtransform  = EPSBackTransform_Default;
1436:   eps->ops->setdefaultst   = EPSSetDefaultST_KrylovSchur;

1438:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur);
1439:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur);
1440:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur);
1441:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur);
1442:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur);
1443:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur);
1444:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur);
1445:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur);
1446:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur);
1447:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur);
1448:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur);
1449:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur);
1450:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur);
1451:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur);
1452:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur);
1453:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur);
1454:   PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur);
1455:   return(0);
1456: }