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: }