Actual source code: epssetup.c
slepc-3.18.1 2022-11-02
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: EPS routines related to problem setup
12: */
14: #include <slepc/private/epsimpl.h>
16: /*
17: Let the solver choose the ST type that should be used by default,
18: otherwise set it to SHIFT.
19: This is called at EPSSetFromOptions (before STSetFromOptions)
20: and also at EPSSetUp (in case EPSSetFromOptions was not called).
21: */
22: PetscErrorCode EPSSetDefaultST(EPS eps)
23: {
24: PetscTryTypeMethod(eps,setdefaultst);
25: if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSHIFT);
26: return 0;
27: }
29: /*
30: This is done by preconditioned eigensolvers that use the PC only.
31: It sets STPRECOND with KSPPREONLY.
32: */
33: PetscErrorCode EPSSetDefaultST_Precond(EPS eps)
34: {
35: KSP ksp;
37: if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STPRECOND);
38: STGetKSP(eps->st,&ksp);
39: if (!((PetscObject)ksp)->type_name) KSPSetType(ksp,KSPPREONLY);
40: return 0;
41: }
43: /*
44: This is done by preconditioned eigensolvers that can also use the KSP.
45: It sets STPRECOND with the default KSP (GMRES) and maxit=5.
46: */
47: PetscErrorCode EPSSetDefaultST_GMRES(EPS eps)
48: {
49: KSP ksp;
51: if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STPRECOND);
52: STPrecondSetKSPHasMat(eps->st,PETSC_TRUE);
53: STGetKSP(eps->st,&ksp);
54: if (!((PetscObject)ksp)->type_name) {
55: KSPSetType(ksp,KSPGMRES);
56: KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,5);
57: }
58: return 0;
59: }
61: #if defined(SLEPC_HAVE_SCALAPACK) || defined(SLEPC_HAVE_ELPA) || defined(SLEPC_HAVE_ELEMENTAL) || defined(SLEPC_HAVE_EVSL)
62: /*
63: This is for direct eigensolvers that work with A and B directly, so
64: no need to factorize B.
65: */
66: PetscErrorCode EPSSetDefaultST_NoFactor(EPS eps)
67: {
68: KSP ksp;
69: PC pc;
71: if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSHIFT);
72: STGetKSP(eps->st,&ksp);
73: if (!((PetscObject)ksp)->type_name) KSPSetType(ksp,KSPPREONLY);
74: KSPGetPC(ksp,&pc);
75: if (!((PetscObject)pc)->type_name) PCSetType(pc,PCNONE);
76: return 0;
77: }
78: #endif
80: /*
81: Check that the ST selected by the user is compatible with the EPS solver and options
82: */
83: PetscErrorCode EPSCheckCompatibleST(EPS eps)
84: {
85: PetscBool precond,shift,sinvert,cayley,lyapii;
86: #if defined(PETSC_USE_COMPLEX)
87: PetscScalar sigma;
88: #endif
90: PetscObjectTypeCompare((PetscObject)eps->st,STPRECOND,&precond);
91: PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&shift);
92: PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert);
93: PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&cayley);
95: /* preconditioned eigensolvers */
99: /* harmonic extraction */
102: /* real shifts in Hermitian problems */
103: #if defined(PETSC_USE_COMPLEX)
104: STGetShift(eps->st,&sigma);
106: #endif
108: /* Cayley with PGNHEP */
111: /* make sure that the user does not specify smallest magnitude with shift-and-invert */
112: if ((cayley || sinvert) && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER)) {
113: PetscObjectTypeCompare((PetscObject)eps,EPSLYAPII,&lyapii);
115: }
116: return 0;
117: }
119: /*
120: MatEstimateSpectralRange_EPS: estimate the spectral range [left,right] of a
121: symmetric/Hermitian matrix A using an auxiliary EPS object
122: */
123: PetscErrorCode MatEstimateSpectralRange_EPS(Mat A,PetscReal *left,PetscReal *right)
124: {
125: PetscInt nconv;
126: PetscScalar eig0;
127: PetscReal tol=1e-3,errest=tol;
128: EPS eps;
130: *left = 0.0; *right = 0.0;
131: EPSCreate(PetscObjectComm((PetscObject)A),&eps);
132: EPSSetOptionsPrefix(eps,"eps_filter_");
133: EPSSetOperators(eps,A,NULL);
134: EPSSetProblemType(eps,EPS_HEP);
135: EPSSetTolerances(eps,tol,50);
136: EPSSetConvergenceTest(eps,EPS_CONV_ABS);
137: EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);
138: EPSSolve(eps);
139: EPSGetConverged(eps,&nconv);
140: if (nconv>0) {
141: EPSGetEigenvalue(eps,0,&eig0,NULL);
142: EPSGetErrorEstimate(eps,0,&errest);
143: } else eig0 = eps->eigr[0];
144: *left = PetscRealPart(eig0)-errest;
145: EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
146: EPSSolve(eps);
147: EPSGetConverged(eps,&nconv);
148: if (nconv>0) {
149: EPSGetEigenvalue(eps,0,&eig0,NULL);
150: EPSGetErrorEstimate(eps,0,&errest);
151: } else eig0 = eps->eigr[0];
152: *right = PetscRealPart(eig0)+errest;
153: EPSDestroy(&eps);
154: return 0;
155: }
157: /*
158: EPSSetUpSort_Basic: configure the EPS sorting criterion according to 'which'
159: */
160: PetscErrorCode EPSSetUpSort_Basic(EPS eps)
161: {
162: switch (eps->which) {
163: case EPS_LARGEST_MAGNITUDE:
164: eps->sc->comparison = SlepcCompareLargestMagnitude;
165: eps->sc->comparisonctx = NULL;
166: break;
167: case EPS_SMALLEST_MAGNITUDE:
168: eps->sc->comparison = SlepcCompareSmallestMagnitude;
169: eps->sc->comparisonctx = NULL;
170: break;
171: case EPS_LARGEST_REAL:
172: eps->sc->comparison = SlepcCompareLargestReal;
173: eps->sc->comparisonctx = NULL;
174: break;
175: case EPS_SMALLEST_REAL:
176: eps->sc->comparison = SlepcCompareSmallestReal;
177: eps->sc->comparisonctx = NULL;
178: break;
179: case EPS_LARGEST_IMAGINARY:
180: eps->sc->comparison = SlepcCompareLargestImaginary;
181: eps->sc->comparisonctx = NULL;
182: break;
183: case EPS_SMALLEST_IMAGINARY:
184: eps->sc->comparison = SlepcCompareSmallestImaginary;
185: eps->sc->comparisonctx = NULL;
186: break;
187: case EPS_TARGET_MAGNITUDE:
188: eps->sc->comparison = SlepcCompareTargetMagnitude;
189: eps->sc->comparisonctx = &eps->target;
190: break;
191: case EPS_TARGET_REAL:
192: eps->sc->comparison = SlepcCompareTargetReal;
193: eps->sc->comparisonctx = &eps->target;
194: break;
195: case EPS_TARGET_IMAGINARY:
196: #if defined(PETSC_USE_COMPLEX)
197: eps->sc->comparison = SlepcCompareTargetImaginary;
198: eps->sc->comparisonctx = &eps->target;
199: #endif
200: break;
201: case EPS_ALL:
202: eps->sc->comparison = SlepcCompareSmallestReal;
203: eps->sc->comparisonctx = NULL;
204: break;
205: case EPS_WHICH_USER:
206: break;
207: }
208: eps->sc->map = NULL;
209: eps->sc->mapobj = NULL;
210: return 0;
211: }
213: /*
214: EPSSetUpSort_Default: configure both EPS and DS sorting criterion
215: */
216: PetscErrorCode EPSSetUpSort_Default(EPS eps)
217: {
218: SlepcSC sc;
219: PetscBool istrivial;
221: /* fill sorting criterion context */
222: EPSSetUpSort_Basic(eps);
223: /* fill sorting criterion for DS */
224: DSGetSlepcSC(eps->ds,&sc);
225: RGIsTrivial(eps->rg,&istrivial);
226: sc->rg = istrivial? NULL: eps->rg;
227: sc->comparison = eps->sc->comparison;
228: sc->comparisonctx = eps->sc->comparisonctx;
229: sc->map = SlepcMap_ST;
230: sc->mapobj = (PetscObject)eps->st;
231: return 0;
232: }
234: /*@
235: EPSSetUp - Sets up all the internal data structures necessary for the
236: execution of the eigensolver. Then calls STSetUp() for any set-up
237: operations associated to the ST object.
239: Collective on eps
241: Input Parameter:
242: . eps - eigenproblem solver context
244: Notes:
245: This function need not be called explicitly in most cases, since EPSSolve()
246: calls it. It can be useful when one wants to measure the set-up time
247: separately from the solve time.
249: Level: developer
251: .seealso: EPSCreate(), EPSSolve(), EPSDestroy(), STSetUp(), EPSSetInitialSpace()
252: @*/
253: PetscErrorCode EPSSetUp(EPS eps)
254: {
255: Mat A,B;
256: PetscInt k,nmat;
257: PetscBool flg;
260: if (eps->state) return 0;
261: PetscLogEventBegin(EPS_SetUp,eps,0,0,0);
263: /* reset the convergence flag from the previous solves */
264: eps->reason = EPS_CONVERGED_ITERATING;
266: /* Set default solver type (EPSSetFromOptions was not called) */
267: if (!((PetscObject)eps)->type_name) EPSSetType(eps,EPSKRYLOVSCHUR);
268: if (!eps->st) EPSGetST(eps,&eps->st);
269: EPSSetDefaultST(eps);
271: STSetTransform(eps->st,PETSC_TRUE);
272: if (eps->useds && !eps->ds) EPSGetDS(eps,&eps->ds);
273: if (eps->twosided) {
275: }
276: if (!eps->rg) EPSGetRG(eps,&eps->rg);
277: if (!((PetscObject)eps->rg)->type_name) RGSetType(eps->rg,RGINTERVAL);
279: /* Set problem dimensions */
280: STGetNumMatrices(eps->st,&nmat);
282: STMatGetSize(eps->st,&eps->n,NULL);
283: STMatGetLocalSize(eps->st,&eps->nloc,NULL);
285: /* Set default problem type */
286: if (!eps->problem_type) {
287: if (nmat==1) EPSSetProblemType(eps,EPS_NHEP);
288: else EPSSetProblemType(eps,EPS_GNHEP);
289: } else if (nmat==1 && eps->isgeneralized) {
290: PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n");
291: eps->isgeneralized = PETSC_FALSE;
292: eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;
295: if (eps->nev > eps->n) eps->nev = eps->n;
296: if (eps->ncv > eps->n) eps->ncv = eps->n;
298: /* check some combinations of eps->which */
301: /* initialization of matrix norms */
302: if (eps->conv==EPS_CONV_NORM) {
303: if (!eps->nrma) {
304: STGetMatrix(eps->st,0,&A);
305: MatNorm(A,NORM_INFINITY,&eps->nrma);
306: }
307: if (nmat>1 && !eps->nrmb) {
308: STGetMatrix(eps->st,1,&B);
309: MatNorm(B,NORM_INFINITY,&eps->nrmb);
310: }
311: }
313: /* call specific solver setup */
314: PetscUseTypeMethod(eps,setup);
316: /* if purification is set, check that it really makes sense */
317: if (eps->purify) {
318: if (eps->categ==EPS_CATEGORY_PRECOND || eps->categ==EPS_CATEGORY_CONTOUR) eps->purify = PETSC_FALSE;
319: else {
320: if (!eps->isgeneralized) eps->purify = PETSC_FALSE;
321: else if (!eps->ishermitian && !eps->ispositive) eps->purify = PETSC_FALSE;
322: else {
323: PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg);
324: if (flg) eps->purify = PETSC_FALSE;
325: }
326: }
327: }
329: /* set tolerance if not yet set */
330: if (eps->tol==PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL;
332: /* set up sorting criterion */
333: PetscTryTypeMethod(eps,setupsort);
335: /* Build balancing matrix if required */
336: if (eps->balance!=EPS_BALANCE_USER) {
337: STSetBalanceMatrix(eps->st,NULL);
338: if (!eps->ishermitian && (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE)) {
339: if (!eps->D) BVCreateVec(eps->V,&eps->D);
340: EPSBuildBalance_Krylov(eps);
341: STSetBalanceMatrix(eps->st,eps->D);
342: }
343: }
345: /* Setup ST */
346: STSetUp(eps->st);
347: EPSCheckCompatibleST(eps);
349: /* process deflation and initial vectors */
350: if (eps->nds<0) {
351: k = -eps->nds;
352: BVInsertConstraints(eps->V,&k,eps->defl);
353: SlepcBasisDestroy_Private(&eps->nds,&eps->defl);
354: eps->nds = k;
355: STCheckNullSpace(eps->st,eps->V);
356: }
357: if (eps->nini<0) {
358: k = -eps->nini;
360: BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE);
361: SlepcBasisDestroy_Private(&eps->nini,&eps->IS);
362: eps->nini = k;
363: }
364: if (eps->twosided && eps->ninil<0) {
365: k = -eps->ninil;
367: BVInsertVecs(eps->W,0,&k,eps->ISL,PETSC_TRUE);
368: SlepcBasisDestroy_Private(&eps->ninil,&eps->ISL);
369: eps->ninil = k;
370: }
372: PetscLogEventEnd(EPS_SetUp,eps,0,0,0);
373: eps->state = EPS_STATE_SETUP;
374: return 0;
375: }
377: /*@
378: EPSSetOperators - Sets the matrices associated with the eigenvalue problem.
380: Collective on eps
382: Input Parameters:
383: + eps - the eigenproblem solver context
384: . A - the matrix associated with the eigensystem
385: - B - the second matrix in the case of generalized eigenproblems
387: Notes:
388: To specify a standard eigenproblem, use NULL for parameter B.
390: It must be called before EPSSetUp(). If it is called again after EPSSetUp() and
391: the matrix sizes have changed then the EPS object is reset.
393: Level: beginner
395: .seealso: EPSSolve(), EPSSetUp(), EPSReset(), EPSGetST(), STGetMatrix()
396: @*/
397: PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)
398: {
399: PetscInt m,n,m0,mloc,nloc,mloc0,nmat;
400: Mat mat[2];
408: /* Check matrix sizes */
409: MatGetSize(A,&m,&n);
410: MatGetLocalSize(A,&mloc,&nloc);
413: if (B) {
414: MatGetSize(B,&m0,&n);
415: MatGetLocalSize(B,&mloc0,&nloc);
420: }
421: if (eps->state && (n!=eps->n || nloc!=eps->nloc)) EPSReset(eps);
422: eps->nrma = 0.0;
423: eps->nrmb = 0.0;
424: if (!eps->st) EPSGetST(eps,&eps->st);
425: mat[0] = A;
426: if (B) {
427: mat[1] = B;
428: nmat = 2;
429: } else nmat = 1;
430: STSetMatrices(eps->st,nmat,mat);
431: eps->state = EPS_STATE_INITIAL;
432: return 0;
433: }
435: /*@
436: EPSGetOperators - Gets the matrices associated with the eigensystem.
438: Collective on eps
440: Input Parameter:
441: . eps - the EPS context
443: Output Parameters:
444: + A - the matrix associated with the eigensystem
445: - B - the second matrix in the case of generalized eigenproblems
447: Note:
448: Does not increase the reference count of the matrices, so you should not destroy them.
450: Level: intermediate
452: .seealso: EPSSolve(), EPSGetST(), STGetMatrix(), STSetMatrices()
453: @*/
454: PetscErrorCode EPSGetOperators(EPS eps,Mat *A,Mat *B)
455: {
456: ST st;
457: PetscInt k;
460: EPSGetST(eps,&st);
461: STGetNumMatrices(st,&k);
462: if (A) {
463: if (k<1) *A = NULL;
464: else STGetMatrix(st,0,A);
465: }
466: if (B) {
467: if (k<2) *B = NULL;
468: else STGetMatrix(st,1,B);
469: }
470: return 0;
471: }
473: /*@C
474: EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
475: space.
477: Collective on eps
479: Input Parameters:
480: + eps - the eigenproblem solver context
481: . n - number of vectors
482: - v - set of basis vectors of the deflation space
484: Notes:
485: When a deflation space is given, the eigensolver seeks the eigensolution
486: in the restriction of the problem to the orthogonal complement of this
487: space. This can be used for instance in the case that an invariant
488: subspace is known beforehand (such as the nullspace of the matrix).
490: These vectors do not persist from one EPSSolve() call to the other, so the
491: deflation space should be set every time.
493: The vectors do not need to be mutually orthonormal, since they are explicitly
494: orthonormalized internally.
496: Level: intermediate
498: .seealso: EPSSetInitialSpace()
499: @*/
500: PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec v[])
501: {
505: if (n>0) {
508: }
509: SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl);
510: if (n>0) eps->state = EPS_STATE_INITIAL;
511: return 0;
512: }
514: /*@C
515: EPSSetInitialSpace - Specify a basis of vectors that constitute the initial
516: space, that is, the subspace from which the solver starts to iterate.
518: Collective on eps
520: Input Parameters:
521: + eps - the eigenproblem solver context
522: . n - number of vectors
523: - is - set of basis vectors of the initial space
525: Notes:
526: Some solvers start to iterate on a single vector (initial vector). In that case,
527: the other vectors are ignored.
529: These vectors do not persist from one EPSSolve() call to the other, so the
530: initial space should be set every time.
532: The vectors do not need to be mutually orthonormal, since they are explicitly
533: orthonormalized internally.
535: Common usage of this function is when the user can provide a rough approximation
536: of the wanted eigenspace. Then, convergence may be faster.
538: Level: intermediate
540: .seealso: EPSSetLeftInitialSpace(), EPSSetDeflationSpace()
541: @*/
542: PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[])
543: {
547: if (n>0) {
550: }
551: SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS);
552: if (n>0) eps->state = EPS_STATE_INITIAL;
553: return 0;
554: }
556: /*@C
557: EPSSetLeftInitialSpace - Specify a basis of vectors that constitute the left
558: initial space, used by two-sided solvers to start the left subspace.
560: Collective on eps
562: Input Parameters:
563: + eps - the eigenproblem solver context
564: . n - number of vectors
565: - isl - set of basis vectors of the left initial space
567: Notes:
568: Left initial vectors are used to initiate the left search space in two-sided
569: eigensolvers. Users should pass here an approximation of the left eigenspace,
570: if available.
572: The same comments in EPSSetInitialSpace() are applicable here.
574: Level: intermediate
576: .seealso: EPSSetInitialSpace(), EPSSetTwoSided()
577: @*/
578: PetscErrorCode EPSSetLeftInitialSpace(EPS eps,PetscInt n,Vec isl[])
579: {
583: if (n>0) {
586: }
587: SlepcBasisReference_Private(n,isl,&eps->ninil,&eps->ISL);
588: if (n>0) eps->state = EPS_STATE_INITIAL;
589: return 0;
590: }
592: /*
593: EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
594: by the user. This is called at setup.
595: */
596: PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
597: {
598: PetscBool krylov;
600: if (*ncv!=PETSC_DEFAULT) { /* ncv set */
601: PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,"");
602: if (krylov) {
604: } else {
606: }
607: } else if (*mpd!=PETSC_DEFAULT) { /* mpd set */
608: *ncv = PetscMin(eps->n,nev+(*mpd));
609: } else { /* neither set: defaults depend on nev being small or large */
610: if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
611: else {
612: *mpd = 500;
613: *ncv = PetscMin(eps->n,nev+(*mpd));
614: }
615: }
616: if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
617: return 0;
618: }
620: /*@
621: EPSAllocateSolution - Allocate memory storage for common variables such
622: as eigenvalues and eigenvectors.
624: Collective on eps
626: Input Parameters:
627: + eps - eigensolver context
628: - extra - number of additional positions, used for methods that require a
629: working basis slightly larger than ncv
631: Developer Notes:
632: This is SLEPC_EXTERN because it may be required by user plugin EPS
633: implementations.
635: Level: developer
637: .seealso: EPSSetUp()
638: @*/
639: PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
640: {
641: PetscInt oldsize,requested;
642: PetscRandom rand;
643: Vec t;
645: requested = eps->ncv + extra;
647: /* oldsize is zero if this is the first time setup is called */
648: BVGetSizes(eps->V,NULL,NULL,&oldsize);
650: /* allocate space for eigenvalues and friends */
651: if (requested != oldsize || !eps->eigr) {
652: PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm);
653: PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm);
654: }
656: /* workspace for the case of arbitrary selection */
657: if (eps->arbitrary) {
658: if (eps->rr) PetscFree2(eps->rr,eps->ri);
659: PetscMalloc2(requested,&eps->rr,requested,&eps->ri);
660: }
662: /* allocate V */
663: if (!eps->V) EPSGetBV(eps,&eps->V);
664: if (!oldsize) {
665: if (!((PetscObject)(eps->V))->type_name) BVSetType(eps->V,BVSVEC);
666: STMatCreateVecsEmpty(eps->st,&t,NULL);
667: BVSetSizesFromVec(eps->V,t,requested);
668: VecDestroy(&t);
669: } else BVResize(eps->V,requested,PETSC_FALSE);
671: /* allocate W */
672: if (eps->twosided) {
673: BVGetRandomContext(eps->V,&rand); /* make sure the random context is available when duplicating */
674: BVDestroy(&eps->W);
675: BVDuplicate(eps->V,&eps->W);
676: }
677: return 0;
678: }