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: {
96: PetscErrorCode ierr;
97: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
98: PetscBool estimaterange=PETSC_TRUE;
99: PetscReal rleft,rright;
100: Mat A;
103: EPSCheckHermitianCondition(eps,PETSC_TRUE," with polynomial filter");
104: EPSCheckStandardCondition(eps,PETSC_TRUE," with polynomial filter");
105: 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");
106: EPSCheckUnsupportedCondition(eps,EPS_FEATURE_ARBITRARY | EPS_FEATURE_REGION | EPS_FEATURE_EXTRACTION,PETSC_TRUE," with polynomial filter");
107: if (eps->tol==PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL*1e-2; /* use tighter tolerance */
108: STFilterSetInterval(eps->st,eps->inta,eps->intb);
109: if (!ctx->estimatedrange) {
110: STFilterGetRange(eps->st,&rleft,&rright);
111: estimaterange = (!rleft && !rright)? PETSC_TRUE: PETSC_FALSE;
112: }
113: if (estimaterange) { /* user did not set a range */
114: STGetMatrix(eps->st,0,&A);
115: EstimateRange(A,&rleft,&rright);
116: PetscInfo2(eps,"Setting eigenvalue range to [%g,%g]\n",(double)rleft,(double)rright);
117: STFilterSetRange(eps->st,rleft,rright);
118: ctx->estimatedrange = PETSC_TRUE;
119: }
120: if (eps->ncv==PETSC_DEFAULT && eps->nev==1) eps->nev = 40; /* user did not provide nev estimation */
121: EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
122: if (eps->ncv>eps->nev+eps->mpd) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must not be larger than nev+mpd");
123: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
124: return(0);
125: }
127: PetscErrorCode EPSSetUp_KrylovSchur(EPS eps)128: {
129: PetscErrorCode ierr;
130: PetscReal eta;
131: PetscBool isfilt=PETSC_FALSE;
132: BVOrthogType otype;
133: BVOrthogBlockType obtype;
134: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
135: enum { EPS_KS_DEFAULT,EPS_KS_SYMM,EPS_KS_SLICE,EPS_KS_FILTER,EPS_KS_INDEF,EPS_KS_TWOSIDED } variant;
138: if (eps->which==EPS_ALL) { /* default values in case of spectrum slicing or polynomial filter */
139: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
140: if (isfilt) {
141: EPSSetUp_KrylovSchur_Filter(eps);
142: } else {
143: EPSSetUp_KrylovSchur_Slice(eps);
144: }
145: } else {
146: EPSSetDimensions_Default(eps,eps->nev,&eps->ncv,&eps->mpd);
147: if (eps->ncv>eps->nev+eps->mpd) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must not be larger than nev+mpd");
148: if (eps->max_it==PETSC_DEFAULT) eps->max_it = PetscMax(100,2*eps->n/eps->ncv);
149: if (!eps->which) { EPSSetWhichEigenpairs_Default(eps); }
150: }
151: if (!ctx->lock && eps->mpd<eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
153: EPSCheckDefiniteCondition(eps,eps->arbitrary," with arbitrary selection of eigenpairs");
155: if (eps->extraction!=EPS_RITZ && eps->extraction!=EPS_HARMONIC) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
157: if (!ctx->keep) ctx->keep = 0.5;
159: EPSAllocateSolution(eps,1);
160: EPS_SetInnerProduct(eps);
161: if (eps->arbitrary) {
162: EPSSetWorkVecs(eps,2);
163: } else if (eps->ishermitian && !eps->ispositive){
164: EPSSetWorkVecs(eps,1);
165: }
167: /* dispatch solve method */
168: if (eps->ishermitian) {
169: if (eps->which==EPS_ALL) {
170: if (eps->isgeneralized && !eps->ispositive) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Spectrum slicing not implemented for indefinite problems");
171: else variant = isfilt? EPS_KS_FILTER: EPS_KS_SLICE;
172: } else if (eps->isgeneralized && !eps->ispositive) {
173: variant = EPS_KS_INDEF;
174: } else {
175: switch (eps->extraction) {
176: case EPS_RITZ: variant = EPS_KS_SYMM; break;
177: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
178: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
179: }
180: }
181: } else if (eps->twosided) {
182: variant = EPS_KS_TWOSIDED;
183: } else {
184: switch (eps->extraction) {
185: case EPS_RITZ: variant = EPS_KS_DEFAULT; break;
186: case EPS_HARMONIC: variant = EPS_KS_DEFAULT; break;
187: default: SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Unsupported extraction type");
188: }
189: }
190: switch (variant) {
191: case EPS_KS_DEFAULT:
192: eps->ops->solve = EPSSolve_KrylovSchur_Default;
193: eps->ops->computevectors = EPSComputeVectors_Schur;
194: DSSetType(eps->ds,DSNHEP);
195: DSAllocate(eps->ds,eps->ncv+1);
196: break;
197: case EPS_KS_SYMM:
198: case EPS_KS_FILTER:
199: eps->ops->solve = EPSSolve_KrylovSchur_Symm;
200: eps->ops->computevectors = EPSComputeVectors_Hermitian;
201: DSSetType(eps->ds,DSHEP);
202: DSSetCompact(eps->ds,PETSC_TRUE);
203: DSSetExtraRow(eps->ds,PETSC_TRUE);
204: DSAllocate(eps->ds,eps->ncv+1);
205: break;
206: case EPS_KS_SLICE:
207: eps->ops->solve = EPSSolve_KrylovSchur_Slice;
208: eps->ops->computevectors = EPSComputeVectors_Slice;
209: break;
210: case EPS_KS_INDEF:
211: eps->ops->solve = EPSSolve_KrylovSchur_Indefinite;
212: eps->ops->computevectors = EPSComputeVectors_Indefinite;
213: DSSetType(eps->ds,DSGHIEP);
214: DSSetCompact(eps->ds,PETSC_TRUE);
215: DSAllocate(eps->ds,eps->ncv+1);
216: /* force reorthogonalization for pseudo-Lanczos */
217: BVGetOrthogonalization(eps->V,&otype,NULL,&eta,&obtype);
218: BVSetOrthogonalization(eps->V,otype,BV_ORTHOG_REFINE_ALWAYS,eta,obtype);
219: break;
220: case EPS_KS_TWOSIDED:
221: eps->ops->solve = EPSSolve_KrylovSchur_TwoSided;
222: eps->ops->computevectors = EPSComputeVectors_Schur;
223: DSSetType(eps->ds,DSNHEP);
224: DSAllocate(eps->ds,eps->ncv+1);
225: DSSetType(eps->dsts,DSNHEP);
226: DSAllocate(eps->dsts,eps->ncv+1);
227: break;
228: default: SETERRQ(PetscObjectComm((PetscObject)eps),1,"Unexpected error");
229: }
230: return(0);
231: }
233: PetscErrorCode EPSSetUpSort_KrylovSchur(EPS eps)234: {
235: PetscErrorCode ierr;
236: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
237: SlepcSC sc;
238: PetscBool isfilt;
241: EPSSetUpSort_Default(eps);
242: if (eps->which==EPS_ALL) {
243: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
244: if (isfilt) {
245: DSGetSlepcSC(eps->ds,&sc);
246: sc->rg = NULL;
247: sc->comparison = SlepcCompareLargestReal;
248: sc->comparisonctx = NULL;
249: sc->map = NULL;
250: sc->mapobj = NULL;
251: } else {
252: if (!ctx->global && ctx->sr->numEigs>0) {
253: DSGetSlepcSC(eps->ds,&sc);
254: sc->rg = NULL;
255: sc->comparison = SlepcCompareLargestMagnitude;
256: sc->comparisonctx = NULL;
257: sc->map = NULL;
258: sc->mapobj = NULL;
259: }
260: }
261: }
262: return(0);
263: }
265: PetscErrorCode EPSSolve_KrylovSchur_Default(EPS eps)266: {
267: PetscErrorCode ierr;
268: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
269: PetscInt i,j,*pj,k,l,nv,ld,nconv;
270: Mat U,Op;
271: PetscScalar *S,*Q,*g;
272: PetscReal beta,gamma=1.0;
273: PetscBool breakdown,harmonic;
276: DSGetLeadingDimension(eps->ds,&ld);
277: harmonic = (eps->extraction==EPS_HARMONIC || eps->extraction==EPS_REFINED_HARMONIC)?PETSC_TRUE:PETSC_FALSE;
278: if (harmonic) { PetscMalloc1(ld,&g); }
279: if (eps->arbitrary) pj = &j;
280: else pj = NULL;
282: /* Get the starting Arnoldi vector */
283: EPSGetStartVector(eps,0,NULL);
284: l = 0;
286: /* Restart loop */
287: while (eps->reason == EPS_CONVERGED_ITERATING) {
288: eps->its++;
290: /* Compute an nv-step Arnoldi factorization */
291: nv = PetscMin(eps->nconv+eps->mpd,eps->ncv);
292: STGetOperator(eps->st,&Op);
293: DSGetArray(eps->ds,DS_MAT_A,&S);
294: BVMatArnoldi(eps->V,Op,S,ld,eps->nconv+l,&nv,&beta,&breakdown);
295: DSRestoreArray(eps->ds,DS_MAT_A,&S);
296: STRestoreOperator(eps->st,&Op);
297: DSSetDimensions(eps->ds,nv,0,eps->nconv,eps->nconv+l);
298: if (l==0) {
299: DSSetState(eps->ds,DS_STATE_INTERMEDIATE);
300: } else {
301: DSSetState(eps->ds,DS_STATE_RAW);
302: }
303: BVSetActiveColumns(eps->V,eps->nconv,nv);
305: /* Compute translation of Krylov decomposition if harmonic extraction used */
306: if (harmonic) {
307: DSTranslateHarmonic(eps->ds,eps->target,beta,PETSC_FALSE,g,&gamma);
308: }
310: /* Solve projected problem */
311: DSSolve(eps->ds,eps->eigr,eps->eigi);
312: if (eps->arbitrary) {
313: EPSGetArbitraryValues(eps,eps->rr,eps->ri);
314: j=1;
315: }
316: DSSort(eps->ds,eps->eigr,eps->eigi,eps->rr,eps->ri,pj);
317: DSSynchronize(eps->ds,eps->eigr,eps->eigi);
319: /* Check convergence */
320: EPSKrylovConvergence(eps,PETSC_FALSE,eps->nconv,nv-eps->nconv,beta,0.0,gamma,&k);
321: (*eps->stopping)(eps,eps->its,eps->max_it,k,eps->nev,&eps->reason,eps->stoppingctx);
322: nconv = k;
324: /* Update l */
325: if (eps->reason != EPS_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
326: else {
327: l = PetscMax(1,(PetscInt)((nv-k)*ctx->keep));
328: #if !defined(PETSC_USE_COMPLEX)
329: DSGetArray(eps->ds,DS_MAT_A,&S);
330: if (S[k+l+(k+l-1)*ld] != 0.0) {
331: if (k+l<nv-1) l = l+1;
332: else l = l-1;
333: }
334: DSRestoreArray(eps->ds,DS_MAT_A,&S);
335: #endif
336: }
337: if (!ctx->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged pairs */
338: if (l) { PetscInfo1(eps,"Preparing to restart keeping l=%D vectors\n",l); }
340: if (eps->reason == EPS_CONVERGED_ITERATING) {
341: if (breakdown || k==nv) {
342: /* Start a new Arnoldi factorization */
343: PetscInfo2(eps,"Breakdown in Krylov-Schur method (it=%D norm=%g)\n",eps->its,(double)beta);
344: if (k<eps->nev) {
345: EPSGetStartVector(eps,k,&breakdown);
346: if (breakdown) {
347: eps->reason = EPS_DIVERGED_BREAKDOWN;
348: PetscInfo(eps,"Unable to generate more start vectors\n");
349: }
350: }
351: } else {
352: /* Undo translation of Krylov decomposition */
353: if (harmonic) {
354: DSSetDimensions(eps->ds,nv,0,k,l);
355: DSTranslateHarmonic(eps->ds,0.0,beta,PETSC_TRUE,g,&gamma);
356: /* gamma u^ = u - U*g~ */
357: BVSetActiveColumns(eps->V,0,nv);
358: BVMultColumn(eps->V,-1.0,1.0,nv,g);
359: BVScaleColumn(eps->V,nv,1.0/gamma);
360: BVSetActiveColumns(eps->V,eps->nconv,nv);
361: }
362: /* Prepare the Rayleigh quotient for restart */
363: DSGetArray(eps->ds,DS_MAT_A,&S);
364: DSGetArray(eps->ds,DS_MAT_Q,&Q);
365: for (i=k;i<k+l;i++) {
366: S[k+l+i*ld] = Q[nv-1+i*ld]*beta*gamma;
367: }
368: DSRestoreArray(eps->ds,DS_MAT_A,&S);
369: DSRestoreArray(eps->ds,DS_MAT_Q,&Q);
370: }
371: }
372: /* Update the corresponding vectors V(:,idx) = V*Q(:,idx) */
373: DSGetMat(eps->ds,DS_MAT_Q,&U);
374: BVMultInPlace(eps->V,U,eps->nconv,k+l);
375: MatDestroy(&U);
377: if (eps->reason == EPS_CONVERGED_ITERATING && !breakdown) {
378: BVCopyColumn(eps->V,nv,k+l);
379: }
380: eps->nconv = k;
381: EPSMonitor(eps,eps->its,nconv,eps->eigr,eps->eigi,eps->errest,nv);
382: }
384: if (harmonic) { PetscFree(g); }
385: /* truncate Schur decomposition and change the state to raw so that
386: DSVectors() computes eigenvectors from scratch */
387: DSSetDimensions(eps->ds,eps->nconv,0,0,0);
388: DSSetState(eps->ds,DS_STATE_RAW);
389: return(0);
390: }
392: static PetscErrorCode EPSKrylovSchurSetRestart_KrylovSchur(EPS eps,PetscReal keep)393: {
394: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
397: if (keep==PETSC_DEFAULT) ctx->keep = 0.5;
398: else {
399: 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);
400: ctx->keep = keep;
401: }
402: return(0);
403: }
405: /*@
406: EPSKrylovSchurSetRestart - Sets the restart parameter for the Krylov-Schur
407: method, in particular the proportion of basis vectors that must be kept
408: after restart.
410: Logically Collective on eps
412: Input Parameters:
413: + eps - the eigenproblem solver context
414: - keep - the number of vectors to be kept at restart
416: Options Database Key:
417: . -eps_krylovschur_restart - Sets the restart parameter
419: Notes:
420: Allowed values are in the range [0.1,0.9]. The default is 0.5.
422: Level: advanced
424: .seealso: EPSKrylovSchurGetRestart()
425: @*/
426: PetscErrorCode EPSKrylovSchurSetRestart(EPS eps,PetscReal keep)427: {
433: PetscTryMethod(eps,"EPSKrylovSchurSetRestart_C",(EPS,PetscReal),(eps,keep));
434: return(0);
435: }
437: static PetscErrorCode EPSKrylovSchurGetRestart_KrylovSchur(EPS eps,PetscReal *keep)438: {
439: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
442: *keep = ctx->keep;
443: return(0);
444: }
446: /*@
447: EPSKrylovSchurGetRestart - Gets the restart parameter used in the
448: Krylov-Schur method.
450: Not Collective
452: Input Parameter:
453: . eps - the eigenproblem solver context
455: Output Parameter:
456: . keep - the restart parameter
458: Level: advanced
460: .seealso: EPSKrylovSchurSetRestart()
461: @*/
462: PetscErrorCode EPSKrylovSchurGetRestart(EPS eps,PetscReal *keep)463: {
469: PetscUseMethod(eps,"EPSKrylovSchurGetRestart_C",(EPS,PetscReal*),(eps,keep));
470: return(0);
471: }
473: static PetscErrorCode EPSKrylovSchurSetLocking_KrylovSchur(EPS eps,PetscBool lock)474: {
475: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
478: ctx->lock = lock;
479: return(0);
480: }
482: /*@
483: EPSKrylovSchurSetLocking - Choose between locking and non-locking variants of
484: the Krylov-Schur method.
486: Logically Collective on eps
488: Input Parameters:
489: + eps - the eigenproblem solver context
490: - lock - true if the locking variant must be selected
492: Options Database Key:
493: . -eps_krylovschur_locking - Sets the locking flag
495: Notes:
496: The default is to lock converged eigenpairs when the method restarts.
497: This behaviour can be changed so that all directions are kept in the
498: working subspace even if already converged to working accuracy (the
499: non-locking variant).
501: Level: advanced
503: .seealso: EPSKrylovSchurGetLocking()
504: @*/
505: PetscErrorCode EPSKrylovSchurSetLocking(EPS eps,PetscBool lock)506: {
512: PetscTryMethod(eps,"EPSKrylovSchurSetLocking_C",(EPS,PetscBool),(eps,lock));
513: return(0);
514: }
516: static PetscErrorCode EPSKrylovSchurGetLocking_KrylovSchur(EPS eps,PetscBool *lock)517: {
518: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
521: *lock = ctx->lock;
522: return(0);
523: }
525: /*@
526: EPSKrylovSchurGetLocking - Gets the locking flag used in the Krylov-Schur
527: method.
529: Not Collective
531: Input Parameter:
532: . eps - the eigenproblem solver context
534: Output Parameter:
535: . lock - the locking flag
537: Level: advanced
539: .seealso: EPSKrylovSchurSetLocking()
540: @*/
541: PetscErrorCode EPSKrylovSchurGetLocking(EPS eps,PetscBool *lock)542: {
548: PetscUseMethod(eps,"EPSKrylovSchurGetLocking_C",(EPS,PetscBool*),(eps,lock));
549: return(0);
550: }
552: static PetscErrorCode EPSKrylovSchurSetPartitions_KrylovSchur(EPS eps,PetscInt npart)553: {
554: PetscErrorCode ierr;
555: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
556: PetscMPIInt size;
559: if (ctx->npart!=npart) {
560: if (ctx->commset) { PetscSubcommDestroy(&ctx->subc); }
561: EPSDestroy(&ctx->eps);
562: }
563: if (npart == PETSC_DEFAULT || npart == PETSC_DECIDE) {
564: ctx->npart = 1;
565: } else {
566: MPI_Comm_size(PetscObjectComm((PetscObject)eps),&size);
567: if (npart<1 || npart>size) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of npart");
568: ctx->npart = npart;
569: }
570: eps->state = EPS_STATE_INITIAL;
571: return(0);
572: }
574: /*@
575: EPSKrylovSchurSetPartitions - Sets the number of partitions for the
576: case of doing spectrum slicing for a computational interval with the
577: communicator split in several sub-communicators.
579: Logically Collective on eps
581: Input Parameters:
582: + eps - the eigenproblem solver context
583: - npart - number of partitions
585: Options Database Key:
586: . -eps_krylovschur_partitions <npart> - Sets the number of partitions
588: Notes:
589: By default, npart=1 so all processes in the communicator participate in
590: the processing of the whole interval. If npart>1 then the interval is
591: divided into npart subintervals, each of them being processed by a
592: subset of processes.
594: The interval is split proportionally unless the separation points are
595: specified with EPSKrylovSchurSetSubintervals().
597: Level: advanced
599: .seealso: EPSKrylovSchurSetSubintervals(), EPSSetInterval()
600: @*/
601: PetscErrorCode EPSKrylovSchurSetPartitions(EPS eps,PetscInt npart)602: {
608: PetscTryMethod(eps,"EPSKrylovSchurSetPartitions_C",(EPS,PetscInt),(eps,npart));
609: return(0);
610: }
612: static PetscErrorCode EPSKrylovSchurGetPartitions_KrylovSchur(EPS eps,PetscInt *npart)613: {
614: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
617: *npart = ctx->npart;
618: return(0);
619: }
621: /*@
622: EPSKrylovSchurGetPartitions - Gets the number of partitions of the
623: communicator in case of spectrum slicing.
625: Not Collective
627: Input Parameter:
628: . eps - the eigenproblem solver context
630: Output Parameter:
631: . npart - number of partitions
633: Level: advanced
635: .seealso: EPSKrylovSchurSetPartitions()
636: @*/
637: PetscErrorCode EPSKrylovSchurGetPartitions(EPS eps,PetscInt *npart)638: {
644: PetscUseMethod(eps,"EPSKrylovSchurGetPartitions_C",(EPS,PetscInt*),(eps,npart));
645: return(0);
646: }
648: static PetscErrorCode EPSKrylovSchurSetDetectZeros_KrylovSchur(EPS eps,PetscBool detect)649: {
650: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
653: ctx->detect = detect;
654: eps->state = EPS_STATE_INITIAL;
655: return(0);
656: }
658: /*@
659: EPSKrylovSchurSetDetectZeros - Sets a flag to enforce detection of
660: zeros during the factorizations throughout the spectrum slicing computation.
662: Logically Collective on eps
664: Input Parameters:
665: + eps - the eigenproblem solver context
666: - detect - check for zeros
668: Options Database Key:
669: . -eps_krylovschur_detect_zeros - Check for zeros; this takes an optional
670: bool value (0/1/no/yes/true/false)
672: Notes:
673: A zero in the factorization indicates that a shift coincides with an eigenvalue.
675: This flag is turned off by default, and may be necessary in some cases,
676: especially when several partitions are being used. This feature currently
677: requires an external package for factorizations with support for zero
678: detection, e.g. MUMPS.
680: Level: advanced
682: .seealso: EPSKrylovSchurSetPartitions(), EPSSetInterval()
683: @*/
684: PetscErrorCode EPSKrylovSchurSetDetectZeros(EPS eps,PetscBool detect)685: {
691: PetscTryMethod(eps,"EPSKrylovSchurSetDetectZeros_C",(EPS,PetscBool),(eps,detect));
692: return(0);
693: }
695: static PetscErrorCode EPSKrylovSchurGetDetectZeros_KrylovSchur(EPS eps,PetscBool *detect)696: {
697: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
700: *detect = ctx->detect;
701: return(0);
702: }
704: /*@
705: EPSKrylovSchurGetDetectZeros - Gets the flag that enforces zero detection
706: in spectrum slicing.
708: Not Collective
710: Input Parameter:
711: . eps - the eigenproblem solver context
713: Output Parameter:
714: . detect - whether zeros detection is enforced during factorizations
716: Level: advanced
718: .seealso: EPSKrylovSchurSetDetectZeros()
719: @*/
720: PetscErrorCode EPSKrylovSchurGetDetectZeros(EPS eps,PetscBool *detect)721: {
727: PetscUseMethod(eps,"EPSKrylovSchurGetDetectZeros_C",(EPS,PetscBool*),(eps,detect));
728: return(0);
729: }
731: static PetscErrorCode EPSKrylovSchurSetDimensions_KrylovSchur(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)732: {
733: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
736: if (nev<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of nev. Must be > 0");
737: ctx->nev = nev;
738: if (ncv == PETSC_DECIDE || ncv == PETSC_DEFAULT) {
739: ctx->ncv = PETSC_DEFAULT;
740: } else {
741: if (ncv<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ncv. Must be > 0");
742: ctx->ncv = ncv;
743: }
744: if (mpd == PETSC_DECIDE || mpd == PETSC_DEFAULT) {
745: ctx->mpd = PETSC_DEFAULT;
746: } else {
747: if (mpd<1) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of mpd. Must be > 0");
748: ctx->mpd = mpd;
749: }
750: eps->state = EPS_STATE_INITIAL;
751: return(0);
752: }
754: /*@
755: EPSKrylovSchurSetDimensions - Sets the dimensions used for each subsolve
756: step in case of doing spectrum slicing for a computational interval.
757: The meaning of the parameters is the same as in EPSSetDimensions().
759: Logically Collective on eps
761: Input Parameters:
762: + eps - the eigenproblem solver context
763: . nev - number of eigenvalues to compute
764: . ncv - the maximum dimension of the subspace to be used by the subsolve
765: - mpd - the maximum dimension allowed for the projected problem
767: Options Database Key:
768: + -eps_krylovschur_nev <nev> - Sets the number of eigenvalues
769: . -eps_krylovschur_ncv <ncv> - Sets the dimension of the subspace
770: - -eps_krylovschur_mpd <mpd> - Sets the maximum projected dimension
772: Level: advanced
774: .seealso: EPSKrylovSchurGetDimensions(), EPSSetDimensions(), EPSSetInterval()
775: @*/
776: PetscErrorCode EPSKrylovSchurSetDimensions(EPS eps,PetscInt nev,PetscInt ncv,PetscInt mpd)777: {
785: PetscTryMethod(eps,"EPSKrylovSchurSetDimensions_C",(EPS,PetscInt,PetscInt,PetscInt),(eps,nev,ncv,mpd));
786: return(0);
787: }
789: static PetscErrorCode EPSKrylovSchurGetDimensions_KrylovSchur(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)790: {
791: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
794: if (nev) *nev = ctx->nev;
795: if (ncv) *ncv = ctx->ncv;
796: if (mpd) *mpd = ctx->mpd;
797: return(0);
798: }
800: /*@
801: EPSKrylovSchurGetDimensions - Gets the dimensions used for each subsolve
802: step in case of doing spectrum slicing for a computational interval.
804: Not Collective
806: Input Parameter:
807: . eps - the eigenproblem solver context
809: Output Parameters:
810: + nev - number of eigenvalues to compute
811: . ncv - the maximum dimension of the subspace to be used by the subsolve
812: - mpd - the maximum dimension allowed for the projected problem
814: Level: advanced
816: .seealso: EPSKrylovSchurSetDimensions()
817: @*/
818: PetscErrorCode EPSKrylovSchurGetDimensions(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)819: {
824: PetscUseMethod(eps,"EPSKrylovSchurGetDimensions_C",(EPS,PetscInt*,PetscInt*,PetscInt*),(eps,nev,ncv,mpd));
825: return(0);
826: }
828: static PetscErrorCode EPSKrylovSchurSetSubintervals_KrylovSchur(EPS eps,PetscReal* subint)829: {
830: PetscErrorCode ierr;
831: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
832: PetscInt i;
835: 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()");
836: 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");
837: if (ctx->subintervals) { PetscFree(ctx->subintervals); }
838: PetscMalloc1(ctx->npart+1,&ctx->subintervals);
839: for (i=0;i<ctx->npart+1;i++) ctx->subintervals[i] = subint[i];
840: ctx->subintset = PETSC_TRUE;
841: eps->state = EPS_STATE_INITIAL;
842: return(0);
843: }
845: /*@C
846: EPSKrylovSchurSetSubintervals - Sets the points that delimit the
847: subintervals to be used in spectrum slicing with several partitions.
849: Logically Collective on eps
851: Input Parameters:
852: + eps - the eigenproblem solver context
853: - subint - array of real values specifying subintervals
855: Notes:
856: This function must be called after EPSKrylovSchurSetPartitions(). For npart
857: partitions, the argument subint must contain npart+1 real values sorted in
858: ascending order: subint_0, subint_1, ..., subint_npart, where the first
859: and last values must coincide with the interval endpoints set with
860: EPSSetInterval().
862: The subintervals are then defined by two consecutive points: [subint_0,subint_1],
863: [subint_1,subint_2], and so on.
865: Level: advanced
867: .seealso: EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubintervals(), EPSSetInterval()
868: @*/
869: PetscErrorCode EPSKrylovSchurSetSubintervals(EPS eps,PetscReal *subint)870: {
875: PetscTryMethod(eps,"EPSKrylovSchurSetSubintervals_C",(EPS,PetscReal*),(eps,subint));
876: return(0);
877: }
879: static PetscErrorCode EPSKrylovSchurGetSubintervals_KrylovSchur(EPS eps,PetscReal **subint)880: {
881: PetscErrorCode ierr;
882: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
883: PetscInt i;
886: if (!ctx->subintset) {
887: if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
888: if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
889: }
890: PetscMalloc1(ctx->npart+1,subint);
891: for (i=0;i<=ctx->npart;i++) (*subint)[i] = ctx->subintervals[i];
892: return(0);
893: }
895: /*@C
896: EPSKrylovSchurGetSubintervals - Returns the points that delimit the
897: subintervals used in spectrum slicing with several partitions.
899: Logically Collective on eps
901: Input Parameter:
902: . eps - the eigenproblem solver context
904: Output Parameter:
905: . subint - array of real values specifying subintervals
907: Notes:
908: If the user passed values with EPSKrylovSchurSetSubintervals(), then the
909: same values are returned. Otherwise, the values computed internally are
910: obtained.
912: This function is only available for spectrum slicing runs.
914: The returned array has length npart+1 (see EPSKrylovSchurGetPartitions())
915: and should be freed by the user.
917: Fortran Notes:
918: The calling sequence from Fortran is
919: .vb
920: EPSKrylovSchurGetSubintervals(eps,subint,ierr)
921: double precision subint(npart+1) output
922: .ve
924: Level: advanced
926: .seealso: EPSKrylovSchurSetSubintervals(), EPSKrylovSchurGetPartitions(), EPSSetInterval()
927: @*/
928: PetscErrorCode EPSKrylovSchurGetSubintervals(EPS eps,PetscReal **subint)929: {
935: PetscUseMethod(eps,"EPSKrylovSchurGetSubintervals_C",(EPS,PetscReal**),(eps,subint));
936: return(0);
937: }
939: static PetscErrorCode EPSKrylovSchurGetInertias_KrylovSchur(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)940: {
941: PetscErrorCode ierr;
942: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
943: PetscInt i,numsh;
944: EPS_SR sr = ctx->sr;
947: if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
948: if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
949: switch (eps->state) {
950: case EPS_STATE_INITIAL:
951: break;
952: case EPS_STATE_SETUP:
953: numsh = ctx->npart+1;
954: if (n) *n = numsh;
955: if (shifts) {
956: PetscMalloc1(numsh,shifts);
957: (*shifts)[0] = eps->inta;
958: if (ctx->npart==1) (*shifts)[1] = eps->intb;
959: else for (i=1;i<numsh;i++) (*shifts)[i] = ctx->subintervals[i];
960: }
961: if (inertias) {
962: PetscMalloc1(numsh,inertias);
963: (*inertias)[0] = (sr->dir==1)?sr->inertia0:sr->inertia1;
964: if (ctx->npart==1) (*inertias)[1] = (sr->dir==1)?sr->inertia1:sr->inertia0;
965: else for (i=1;i<numsh;i++) (*inertias)[i] = (*inertias)[i-1]+ctx->nconv_loc[i-1];
966: }
967: break;
968: case EPS_STATE_SOLVED:
969: case EPS_STATE_EIGENVECTORS:
970: numsh = ctx->nshifts;
971: if (n) *n = numsh;
972: if (shifts) {
973: PetscMalloc1(numsh,shifts);
974: for (i=0;i<numsh;i++) (*shifts)[i] = ctx->shifts[i];
975: }
976: if (inertias) {
977: PetscMalloc1(numsh,inertias);
978: for (i=0;i<numsh;i++) (*inertias)[i] = ctx->inertias[i];
979: }
980: break;
981: }
982: return(0);
983: }
985: /*@C
986: EPSKrylovSchurGetInertias - Gets the values of the shifts and their
987: corresponding inertias in case of doing spectrum slicing for a
988: computational interval.
990: Not Collective
992: Input Parameter:
993: . eps - the eigenproblem solver context
995: Output Parameters:
996: + n - number of shifts, including the endpoints of the interval
997: . shifts - the values of the shifts used internally in the solver
998: - inertias - the values of the inertia in each shift
1000: Notes:
1001: If called after EPSSolve(), all shifts used internally by the solver are
1002: returned (including both endpoints and any intermediate ones). If called
1003: before EPSSolve() and after EPSSetUp() then only the information of the
1004: endpoints of subintervals is available.
1006: This function is only available for spectrum slicing runs.
1008: The returned arrays should be freed by the user. Can pass NULL in any of
1009: the two arrays if not required.
1011: Fortran Notes:
1012: The calling sequence from Fortran is
1013: .vb
1014: EPSKrylovSchurGetInertias(eps,n,shifts,inertias,ierr)
1015: integer n
1016: double precision shifts(*)
1017: integer inertias(*)
1018: .ve
1019: The arrays should be at least of length n. The value of n can be determined
1020: by an initial call
1021: .vb
1022: EPSKrylovSchurGetInertias(eps,n,PETSC_NULL_REAL,PETSC_NULL_INTEGER,ierr)
1023: .ve
1025: Level: advanced
1027: .seealso: EPSSetInterval(), EPSKrylovSchurSetSubintervals()
1028: @*/
1029: PetscErrorCode EPSKrylovSchurGetInertias(EPS eps,PetscInt *n,PetscReal **shifts,PetscInt **inertias)1030: {
1036: PetscUseMethod(eps,"EPSKrylovSchurGetInertias_C",(EPS,PetscInt*,PetscReal**,PetscInt**),(eps,n,shifts,inertias));
1037: return(0);
1038: }
1040: static PetscErrorCode EPSKrylovSchurGetSubcommInfo_KrylovSchur(EPS eps,PetscInt *k,PetscInt *n,Vec *v)1041: {
1042: PetscErrorCode ierr;
1043: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1044: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1047: if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1048: if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1049: if (k) *k = (ctx->npart==1)? 0: ctx->subc->color;
1050: if (n) *n = sr->numEigs;
1051: if (v) {
1052: BVCreateVec(sr->V,v);
1053: }
1054: return(0);
1055: }
1057: /*@C
1058: EPSKrylovSchurGetSubcommInfo - Gets information related to the case of
1059: doing spectrum slicing for a computational interval with multiple
1060: communicators.
1062: Collective on the subcommunicator (if v is given)
1064: Input Parameter:
1065: . eps - the eigenproblem solver context
1067: Output Parameters:
1068: + k - index of the subinterval for the calling process
1069: . n - number of eigenvalues found in the k-th subinterval
1070: - v - a vector owned by processes in the subcommunicator with dimensions
1071: compatible for locally computed eigenvectors (or NULL)
1073: Notes:
1074: This function is only available for spectrum slicing runs.
1076: The returned Vec should be destroyed by the user.
1078: Level: advanced
1080: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommPairs()
1081: @*/
1082: PetscErrorCode EPSKrylovSchurGetSubcommInfo(EPS eps,PetscInt *k,PetscInt *n,Vec *v)1083: {
1088: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommInfo_C",(EPS,PetscInt*,PetscInt*,Vec*),(eps,k,n,v));
1089: return(0);
1090: }
1092: static PetscErrorCode EPSKrylovSchurGetSubcommPairs_KrylovSchur(EPS eps,PetscInt i,PetscScalar *eig,Vec v)1093: {
1094: PetscErrorCode ierr;
1095: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1096: EPS_SR sr = ((EPS_KRYLOVSCHUR*)ctx->eps->data)->sr;
1099: EPSCheckSolved(eps,1);
1100: if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1101: if (i<0 || i>=sr->numEigs) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
1102: if (eig) *eig = sr->eigr[sr->perm[i]];
1103: if (v) { BVCopyVec(sr->V,sr->perm[i],v); }
1104: return(0);
1105: }
1107: /*@
1108: EPSKrylovSchurGetSubcommPairs - Gets the i-th eigenpair stored
1109: internally in the subcommunicator to which the calling process belongs.
1111: Collective on the subcommunicator (if v is given)
1113: Input Parameter:
1114: + eps - the eigenproblem solver context
1115: - i - index of the solution
1117: Output Parameters:
1118: + eig - the eigenvalue
1119: - v - the eigenvector
1121: Notes:
1122: It is allowed to pass NULL for v if the eigenvector is not required.
1123: Otherwise, the caller must provide a valid Vec objects, i.e.,
1124: it must be created by the calling program with EPSKrylovSchurGetSubcommInfo().
1126: The index i should be a value between 0 and n-1, where n is the number of
1127: vectors in the local subinterval, see EPSKrylovSchurGetSubcommInfo().
1129: Level: advanced
1131: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo(), EPSKrylovSchurGetSubcommMats()
1132: @*/
1133: PetscErrorCode EPSKrylovSchurGetSubcommPairs(EPS eps,PetscInt i,PetscScalar *eig,Vec v)1134: {
1140: PetscUseMethod(eps,"EPSKrylovSchurGetSubcommPairs_C",(EPS,PetscInt,PetscScalar*,Vec),(eps,i,eig,v));
1141: return(0);
1142: }
1144: static PetscErrorCode EPSKrylovSchurGetSubcommMats_KrylovSchur(EPS eps,Mat *A,Mat *B)1145: {
1146: PetscErrorCode ierr;
1147: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1150: if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1151: if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1152: EPSGetOperators(ctx->eps,A,B);
1153: return(0);
1154: }
1156: /*@C
1157: EPSKrylovSchurGetSubcommMats - Gets the eigenproblem matrices stored
1158: internally in the subcommunicator to which the calling process belongs.
1160: Collective on the subcommunicator
1162: Input Parameter:
1163: . eps - the eigenproblem solver context
1165: Output Parameters:
1166: + A - the matrix associated with the eigensystem
1167: - B - the second matrix in the case of generalized eigenproblems
1169: Notes:
1170: This is the analog of EPSGetOperators(), but returns the matrices distributed
1171: differently (in the subcommunicator rather than in the parent communicator).
1173: These matrices should not be modified by the user.
1175: Level: advanced
1177: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommInfo()
1178: @*/
1179: PetscErrorCode EPSKrylovSchurGetSubcommMats(EPS eps,Mat *A,Mat *B)1180: {
1185: PetscTryMethod(eps,"EPSKrylovSchurGetSubcommMats_C",(EPS,Mat*,Mat*),(eps,A,B));
1186: return(0);
1187: }
1189: static PetscErrorCode EPSKrylovSchurUpdateSubcommMats_KrylovSchur(EPS eps,PetscScalar a,PetscScalar ap,Mat Au,PetscScalar b,PetscScalar bp, Mat Bu,MatStructure str,PetscBool globalup)1190: {
1191: PetscErrorCode ierr;
1192: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data,*subctx;
1193: Mat A,B=NULL,Ag,Bg=NULL;
1194: PetscBool reuse=PETSC_TRUE;
1197: if (!ctx->sr) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Only available in interval computations, see EPSSetInterval()");
1198: if (!eps->state) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"Must call EPSSetUp() first");
1199: EPSGetOperators(eps,&Ag,&Bg);
1200: EPSGetOperators(ctx->eps,&A,&B);
1202: MatScale(A,a);
1203: if (Au) {
1204: MatAXPY(A,ap,Au,str);
1205: }
1206: if (B) MatScale(B,b);
1207: if (Bu) {
1208: MatAXPY(B,bp,Bu,str);
1209: }
1210: EPSSetOperators(ctx->eps,A,B);
1212: /* Update stored matrix state */
1213: subctx = (EPS_KRYLOVSCHUR*)ctx->eps->data;
1214: PetscObjectStateGet((PetscObject)A,&subctx->Astate);
1215: if (B) { PetscObjectStateGet((PetscObject)B,&subctx->Bstate); }
1217: /* Update matrices in the parent communicator if requested by user */
1218: if (globalup) {
1219: if (ctx->npart>1) {
1220: if (!ctx->isrow) {
1221: MatGetOwnershipIS(Ag,&ctx->isrow,&ctx->iscol);
1222: reuse = PETSC_FALSE;
1223: }
1224: if (str==DIFFERENT_NONZERO_PATTERN) reuse = PETSC_FALSE;
1225: if (ctx->submata && !reuse) {
1226: MatDestroyMatrices(1,&ctx->submata);
1227: }
1228: MatCreateSubMatrices(A,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submata);
1229: MatCreateMPIMatConcatenateSeqMat(((PetscObject)Ag)->comm,ctx->submata[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Ag);
1230: if (B) {
1231: if (ctx->submatb && !reuse) {
1232: MatDestroyMatrices(1,&ctx->submatb);
1233: }
1234: MatCreateSubMatrices(B,1,&ctx->isrow,&ctx->iscol,(reuse)?MAT_REUSE_MATRIX:MAT_INITIAL_MATRIX,&ctx->submatb);
1235: MatCreateMPIMatConcatenateSeqMat(((PetscObject)Bg)->comm,ctx->submatb[0],PETSC_DECIDE,MAT_REUSE_MATRIX,&Bg);
1236: }
1237: }
1238: PetscObjectStateGet((PetscObject)Ag,&ctx->Astate);
1239: if (Bg) { PetscObjectStateGet((PetscObject)Bg,&ctx->Bstate); }
1240: }
1241: EPSSetOperators(eps,Ag,Bg);
1242: return(0);
1243: }
1245: /*@
1246: EPSKrylovSchurUpdateSubcommMats - Update the eigenproblem matrices stored
1247: internally in the subcommunicator to which the calling process belongs.
1249: Collective on eps
1251: Input Parameters:
1252: + eps - the eigenproblem solver context
1253: . s - scalar that multiplies the existing A matrix
1254: . a - scalar used in the axpy operation on A
1255: . Au - matrix used in the axpy operation on A
1256: . t - scalar that multiplies the existing B matrix
1257: . b - scalar used in the axpy operation on B
1258: . Bu - matrix used in the axpy operation on B
1259: . str - structure flag
1260: - globalup - flag indicating if global matrices must be updated
1262: Notes:
1263: This function modifies the eigenproblem matrices at the subcommunicator level,
1264: and optionally updates the global matrices in the parent communicator. The updates
1265: are expressed as A <-- s*A + a*Au, B <-- t*B + b*Bu.
1267: It is possible to update one of the matrices, or both.
1269: The matrices Au and Bu must be equal in all subcommunicators.
1271: The str flag is passed to the MatAXPY() operations to perform the updates.
1273: If globalup is true, communication is carried out to reconstruct the updated
1274: matrices in the parent communicator. The user must be warned that if global
1275: matrices are not in sync with subcommunicator matrices, the errors computed
1276: by EPSComputeError() will be wrong even if the computed solution is correct
1277: (the synchronization may be done only once at the end).
1279: Level: advanced
1281: .seealso: EPSSetInterval(), EPSKrylovSchurSetPartitions(), EPSKrylovSchurGetSubcommMats()
1282: @*/
1283: PetscErrorCode EPSKrylovSchurUpdateSubcommMats(EPS eps,PetscScalar s,PetscScalar a,Mat Au,PetscScalar t,PetscScalar b,Mat Bu,MatStructure str,PetscBool globalup)1284: {
1297: PetscTryMethod(eps,"EPSKrylovSchurUpdateSubcommMats_C",(EPS,PetscScalar,PetscScalar,Mat,PetscScalar,PetscScalar,Mat,MatStructure,PetscBool),(eps,s,a,Au,t,b,Bu,str,globalup));
1298: return(0);
1299: }
1301: PetscErrorCode EPSSetFromOptions_KrylovSchur(PetscOptionItems *PetscOptionsObject,EPS eps)1302: {
1303: PetscErrorCode ierr;
1304: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1305: PetscBool flg,lock,b,f1,f2,f3;
1306: PetscReal keep;
1307: PetscInt i,j,k;
1310: PetscOptionsHead(PetscOptionsObject,"EPS Krylov-Schur Options");
1312: PetscOptionsReal("-eps_krylovschur_restart","Proportion of vectors kept after restart","EPSKrylovSchurSetRestart",0.5,&keep,&flg);
1313: if (flg) { EPSKrylovSchurSetRestart(eps,keep); }
1315: PetscOptionsBool("-eps_krylovschur_locking","Choose between locking and non-locking variants","EPSKrylovSchurSetLocking",PETSC_TRUE,&lock,&flg);
1316: if (flg) { EPSKrylovSchurSetLocking(eps,lock); }
1318: i = ctx->npart;
1319: PetscOptionsInt("-eps_krylovschur_partitions","Number of partitions of the communicator for spectrum slicing","EPSKrylovSchurSetPartitions",ctx->npart,&i,&flg);
1320: if (flg) { EPSKrylovSchurSetPartitions(eps,i); }
1322: b = ctx->detect;
1323: PetscOptionsBool("-eps_krylovschur_detect_zeros","Check zeros during factorizations at subinterval boundaries","EPSKrylovSchurSetDetectZeros",ctx->detect,&b,&flg);
1324: if (flg) { EPSKrylovSchurSetDetectZeros(eps,b); }
1326: i = 1;
1327: j = k = PETSC_DECIDE;
1328: PetscOptionsInt("-eps_krylovschur_nev","Number of eigenvalues to compute in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",40,&i,&f1);
1329: PetscOptionsInt("-eps_krylovschur_ncv","Number of basis vectors in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&j,&f2);
1330: PetscOptionsInt("-eps_krylovschur_mpd","Maximum dimension of projected problem in each subsolve (only for spectrum slicing)","EPSKrylovSchurSetDimensions",80,&k,&f3);
1331: if (f1 || f2 || f3) { EPSKrylovSchurSetDimensions(eps,i,j,k); }
1333: PetscOptionsTail();
1334: return(0);
1335: }
1337: PetscErrorCode EPSView_KrylovSchur(EPS eps,PetscViewer viewer)1338: {
1339: PetscErrorCode ierr;
1340: EPS_KRYLOVSCHUR *ctx = (EPS_KRYLOVSCHUR*)eps->data;
1341: PetscBool isascii,isfilt;
1344: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
1345: if (isascii) {
1346: PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*ctx->keep));
1347: PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",ctx->lock?"":"non-");
1348: if (eps->which==EPS_ALL) {
1349: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1350: if (isfilt) {
1351: PetscViewerASCIIPrintf(viewer," using filtering to extract all eigenvalues in an interval\n");
1352: } else {
1353: PetscViewerASCIIPrintf(viewer," doing spectrum slicing with nev=%D, ncv=%D, mpd=%D\n",ctx->nev,ctx->ncv,ctx->mpd);
1354: if (ctx->npart>1) {
1355: PetscViewerASCIIPrintf(viewer," multi-communicator spectrum slicing with %D partitions\n",ctx->npart);
1356: if (ctx->detect) { PetscViewerASCIIPrintf(viewer," detecting zeros when factorizing at subinterval boundaries\n"); }
1357: }
1358: }
1359: }
1360: }
1361: return(0);
1362: }
1364: PetscErrorCode EPSDestroy_KrylovSchur(EPS eps)1365: {
1369: PetscFree(eps->data);
1370: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",NULL);
1371: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",NULL);
1372: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",NULL);
1373: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",NULL);
1374: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",NULL);
1375: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",NULL);
1376: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",NULL);
1377: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",NULL);
1378: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",NULL);
1379: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",NULL);
1380: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",NULL);
1381: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",NULL);
1382: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",NULL);
1383: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",NULL);
1384: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",NULL);
1385: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",NULL);
1386: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",NULL);
1387: return(0);
1388: }
1390: PetscErrorCode EPSReset_KrylovSchur(EPS eps)1391: {
1393: PetscBool isfilt;
1396: PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilt);
1397: if (eps->which==EPS_ALL && !isfilt) {
1398: EPSReset_KrylovSchur_Slice(eps);
1399: }
1400: return(0);
1401: }
1403: PetscErrorCode EPSSetDefaultST_KrylovSchur(EPS eps)1404: {
1408: if (eps->which==EPS_ALL) {
1409: if (!((PetscObject)eps->st)->type_name) {
1410: STSetType(eps->st,STSINVERT);
1411: }
1412: }
1413: return(0);
1414: }
1416: SLEPC_EXTERN PetscErrorCode EPSCreate_KrylovSchur(EPS eps)1417: {
1418: EPS_KRYLOVSCHUR *ctx;
1419: PetscErrorCode ierr;
1422: PetscNewLog(eps,&ctx);
1423: eps->data = (void*)ctx;
1424: ctx->lock = PETSC_TRUE;
1425: ctx->nev = 1;
1426: ctx->ncv = PETSC_DEFAULT;
1427: ctx->mpd = PETSC_DEFAULT;
1428: ctx->npart = 1;
1429: ctx->detect = PETSC_FALSE;
1430: ctx->global = PETSC_TRUE;
1432: eps->useds = PETSC_TRUE;
1434: /* solve and computevectors determined at setup */
1435: eps->ops->setup = EPSSetUp_KrylovSchur;
1436: eps->ops->setupsort = EPSSetUpSort_KrylovSchur;
1437: eps->ops->setfromoptions = EPSSetFromOptions_KrylovSchur;
1438: eps->ops->destroy = EPSDestroy_KrylovSchur;
1439: eps->ops->reset = EPSReset_KrylovSchur;
1440: eps->ops->view = EPSView_KrylovSchur;
1441: eps->ops->backtransform = EPSBackTransform_Default;
1442: eps->ops->setdefaultst = EPSSetDefaultST_KrylovSchur;
1444: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetRestart_C",EPSKrylovSchurSetRestart_KrylovSchur);
1445: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetRestart_C",EPSKrylovSchurGetRestart_KrylovSchur);
1446: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetLocking_C",EPSKrylovSchurSetLocking_KrylovSchur);
1447: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetLocking_C",EPSKrylovSchurGetLocking_KrylovSchur);
1448: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetPartitions_C",EPSKrylovSchurSetPartitions_KrylovSchur);
1449: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetPartitions_C",EPSKrylovSchurGetPartitions_KrylovSchur);
1450: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDetectZeros_C",EPSKrylovSchurSetDetectZeros_KrylovSchur);
1451: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDetectZeros_C",EPSKrylovSchurGetDetectZeros_KrylovSchur);
1452: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetDimensions_C",EPSKrylovSchurSetDimensions_KrylovSchur);
1453: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetDimensions_C",EPSKrylovSchurGetDimensions_KrylovSchur);
1454: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurSetSubintervals_C",EPSKrylovSchurSetSubintervals_KrylovSchur);
1455: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubintervals_C",EPSKrylovSchurGetSubintervals_KrylovSchur);
1456: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetInertias_C",EPSKrylovSchurGetInertias_KrylovSchur);
1457: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommInfo_C",EPSKrylovSchurGetSubcommInfo_KrylovSchur);
1458: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommPairs_C",EPSKrylovSchurGetSubcommPairs_KrylovSchur);
1459: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurGetSubcommMats_C",EPSKrylovSchurGetSubcommMats_KrylovSchur);
1460: PetscObjectComposeFunction((PetscObject)eps,"EPSKrylovSchurUpdateSubcommMats_C",EPSKrylovSchurUpdateSubcommMats_KrylovSchur);
1461: return(0);
1462: }