Actual source code: epssetup.c
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: PetscFunctionBegin;
25: PetscTryTypeMethod(eps,setdefaultst);
26: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSHIFT));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: /*
31: This is done by preconditioned eigensolvers that use the PC only.
32: It sets STPRECOND with KSPPREONLY.
33: */
34: PetscErrorCode EPSSetDefaultST_Precond(EPS eps)
35: {
36: KSP ksp;
38: PetscFunctionBegin;
39: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STPRECOND));
40: PetscCall(STGetKSP(eps->st,&ksp));
41: if (!((PetscObject)ksp)->type_name) PetscCall(KSPSetType(ksp,KSPPREONLY));
42: PetscFunctionReturn(PETSC_SUCCESS);
43: }
45: /*
46: This is done by preconditioned eigensolvers that can also use the KSP.
47: It sets STPRECOND with the default KSP (GMRES) and maxit=5.
48: */
49: PetscErrorCode EPSSetDefaultST_GMRES(EPS eps)
50: {
51: KSP ksp;
53: PetscFunctionBegin;
54: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STPRECOND));
55: PetscCall(STPrecondSetKSPHasMat(eps->st,PETSC_TRUE));
56: PetscCall(STGetKSP(eps->st,&ksp));
57: if (!((PetscObject)ksp)->type_name) {
58: PetscCall(KSPSetType(ksp,KSPGMRES));
59: PetscCall(KSPSetTolerances(ksp,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT,5));
60: }
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
64: #if defined(SLEPC_HAVE_SCALAPACK) || defined(SLEPC_HAVE_ELPA) || defined(SLEPC_HAVE_ELEMENTAL) || defined(SLEPC_HAVE_EVSL)
65: /*
66: This is for direct eigensolvers that work with A and B directly, so
67: no need to factorize B.
68: */
69: PetscErrorCode EPSSetDefaultST_NoFactor(EPS eps)
70: {
71: KSP ksp;
72: PC pc;
74: PetscFunctionBegin;
75: if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSHIFT));
76: PetscCall(STGetKSP(eps->st,&ksp));
77: if (!((PetscObject)ksp)->type_name) PetscCall(KSPSetType(ksp,KSPPREONLY));
78: PetscCall(KSPGetPC(ksp,&pc));
79: if (!((PetscObject)pc)->type_name) PetscCall(PCSetType(pc,PCNONE));
80: PetscFunctionReturn(PETSC_SUCCESS);
81: }
82: #endif
84: /*
85: Check that the ST selected by the user is compatible with the EPS solver and options
86: */
87: static PetscErrorCode EPSCheckCompatibleST(EPS eps)
88: {
89: PetscBool precond,shift,sinvert,cayley,lyapii;
90: #if defined(PETSC_USE_COMPLEX)
91: PetscScalar sigma;
92: #endif
94: PetscFunctionBegin;
95: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STPRECOND,&precond));
96: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&shift));
97: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert));
98: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&cayley));
100: /* preconditioned eigensolvers */
101: PetscCheck(eps->categ!=EPS_CATEGORY_PRECOND || precond,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires ST=PRECOND");
102: PetscCheck(eps->categ==EPS_CATEGORY_PRECOND || !precond,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"STPRECOND is intended for preconditioned eigensolvers only");
104: /* harmonic extraction */
105: PetscCheck(precond || shift || !eps->extraction || eps->extraction==EPS_RITZ,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot use a spectral transformation combined with harmonic extraction");
107: /* real shifts in Hermitian problems */
108: #if defined(PETSC_USE_COMPLEX)
109: PetscCall(STGetShift(eps->st,&sigma));
110: PetscCheck(!eps->ishermitian || PetscImaginaryPart(sigma)==0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Hermitian problems are not compatible with complex shifts");
111: #endif
113: /* Cayley with PGNHEP */
114: PetscCheck(!cayley || eps->problem_type!=EPS_PGNHEP,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cayley spectral transformation is not compatible with PGNHEP");
116: /* make sure that the user does not specify smallest magnitude with shift-and-invert */
117: if ((cayley || sinvert) && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER)) {
118: PetscCall(PetscObjectTypeCompare((PetscObject)eps,EPSLYAPII,&lyapii));
119: PetscCheck(lyapii || eps->which==EPS_TARGET_MAGNITUDE || eps->which==EPS_TARGET_REAL || eps->which==EPS_TARGET_IMAGINARY || eps->which==EPS_ALL || eps->which==EPS_WHICH_USER,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"Shift-and-invert requires a target 'which' (see EPSSetWhichEigenpairs), for instance -st_type sinvert -eps_target 0 -eps_target_magnitude");
120: }
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: /*
125: MatEstimateSpectralRange_EPS: estimate the spectral range [left,right] of a
126: symmetric/Hermitian matrix A using an auxiliary EPS object
127: */
128: PetscErrorCode MatEstimateSpectralRange_EPS(Mat A,PetscReal *left,PetscReal *right)
129: {
130: PetscInt nconv;
131: PetscScalar eig0;
132: PetscReal tol=1e-3,errest=tol;
133: EPS eps;
135: PetscFunctionBegin;
136: *left = 0.0; *right = 0.0;
137: PetscCall(EPSCreate(PetscObjectComm((PetscObject)A),&eps));
138: PetscCall(EPSSetOptionsPrefix(eps,"eps_filter_"));
139: PetscCall(EPSSetOperators(eps,A,NULL));
140: PetscCall(EPSSetProblemType(eps,EPS_HEP));
141: PetscCall(EPSSetTolerances(eps,tol,50));
142: PetscCall(EPSSetConvergenceTest(eps,EPS_CONV_ABS));
143: PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
144: PetscCall(EPSSolve(eps));
145: PetscCall(EPSGetConverged(eps,&nconv));
146: if (nconv>0) {
147: PetscCall(EPSGetEigenvalue(eps,0,&eig0,NULL));
148: PetscCall(EPSGetErrorEstimate(eps,0,&errest));
149: } else eig0 = eps->eigr[0];
150: *left = PetscRealPart(eig0)-errest;
151: PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
152: PetscCall(EPSSolve(eps));
153: PetscCall(EPSGetConverged(eps,&nconv));
154: if (nconv>0) {
155: PetscCall(EPSGetEigenvalue(eps,0,&eig0,NULL));
156: PetscCall(EPSGetErrorEstimate(eps,0,&errest));
157: } else eig0 = eps->eigr[0];
158: *right = PetscRealPart(eig0)+errest;
159: PetscCall(EPSDestroy(&eps));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: /*
164: EPSSetUpSort_Basic: configure the EPS sorting criterion according to 'which'
165: */
166: PetscErrorCode EPSSetUpSort_Basic(EPS eps)
167: {
168: PetscFunctionBegin;
169: switch (eps->which) {
170: case EPS_LARGEST_MAGNITUDE:
171: eps->sc->comparison = SlepcCompareLargestMagnitude;
172: eps->sc->comparisonctx = NULL;
173: break;
174: case EPS_SMALLEST_MAGNITUDE:
175: eps->sc->comparison = SlepcCompareSmallestMagnitude;
176: eps->sc->comparisonctx = NULL;
177: break;
178: case EPS_LARGEST_REAL:
179: eps->sc->comparison = SlepcCompareLargestReal;
180: eps->sc->comparisonctx = NULL;
181: break;
182: case EPS_SMALLEST_REAL:
183: eps->sc->comparison = SlepcCompareSmallestReal;
184: eps->sc->comparisonctx = NULL;
185: break;
186: case EPS_LARGEST_IMAGINARY:
187: eps->sc->comparison = SlepcCompareLargestImaginary;
188: eps->sc->comparisonctx = NULL;
189: break;
190: case EPS_SMALLEST_IMAGINARY:
191: eps->sc->comparison = SlepcCompareSmallestImaginary;
192: eps->sc->comparisonctx = NULL;
193: break;
194: case EPS_TARGET_MAGNITUDE:
195: eps->sc->comparison = SlepcCompareTargetMagnitude;
196: eps->sc->comparisonctx = &eps->target;
197: break;
198: case EPS_TARGET_REAL:
199: eps->sc->comparison = SlepcCompareTargetReal;
200: eps->sc->comparisonctx = &eps->target;
201: break;
202: case EPS_TARGET_IMAGINARY:
203: eps->sc->comparison = SlepcCompareTargetImaginary;
204: eps->sc->comparisonctx = &eps->target;
205: break;
206: case EPS_ALL:
207: eps->sc->comparison = SlepcCompareSmallestReal;
208: eps->sc->comparisonctx = NULL;
209: break;
210: case EPS_WHICH_USER:
211: break;
212: }
213: eps->sc->map = NULL;
214: eps->sc->mapobj = NULL;
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*
219: EPSSetUpSort_Default: configure both EPS and DS sorting criterion
220: */
221: PetscErrorCode EPSSetUpSort_Default(EPS eps)
222: {
223: SlepcSC sc;
224: PetscBool istrivial;
226: PetscFunctionBegin;
227: /* fill sorting criterion context */
228: PetscCall(EPSSetUpSort_Basic(eps));
229: /* fill sorting criterion for DS */
230: PetscCall(DSGetSlepcSC(eps->ds,&sc));
231: PetscCall(RGIsTrivial(eps->rg,&istrivial));
232: sc->rg = istrivial? NULL: eps->rg;
233: sc->comparison = eps->sc->comparison;
234: sc->comparisonctx = eps->sc->comparisonctx;
235: sc->map = SlepcMap_ST;
236: sc->mapobj = (PetscObject)eps->st;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*@
241: EPSSetDSType - Sets the type of the internal `DS` object based on the current
242: settings of the eigenvalue solver.
244: Collective
246: Input Parameter:
247: . eps - the linear eigensolver context
249: Note:
250: This function need not be called explicitly, since it will be called at
251: both `EPSSetFromOptions()` and `EPSSetUp()`.
253: Level: developer
255: .seealso: [](ch:eps), `EPSSetFromOptions()`, `EPSSetUp()`
256: @*/
257: PetscErrorCode EPSSetDSType(EPS eps)
258: {
259: PetscFunctionBegin;
261: PetscTryTypeMethod(eps,setdstype);
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*@
266: EPSSetUp - Sets up all the internal data structures necessary for the
267: execution of the eigensolver. Then calls `STSetUp()` for any set-up
268: operations associated to the internal `ST` object.
270: Collective
272: Input Parameter:
273: . eps - the linear eigensolver context
275: Notes:
276: This function need not be called explicitly in most cases, since `EPSSolve()`
277: calls it. It can be useful when one wants to measure the set-up time
278: separately from the solve time.
280: Level: developer
282: .seealso: [](ch:eps), `EPSCreate()`, `EPSSolve()`, `EPSDestroy()`, `STSetUp()`, `EPSSetInitialSpace()`, `EPSSetDeflationSpace()`
283: @*/
284: PetscErrorCode EPSSetUp(EPS eps)
285: {
286: Mat A,B;
287: PetscInt k,nmat;
288: PetscBool flg,isfilter;
289: EPSStoppingCtx ctx;
291: PetscFunctionBegin;
293: if (eps->state) PetscFunctionReturn(PETSC_SUCCESS);
294: PetscCall(PetscLogEventBegin(EPS_SetUp,eps,0,0,0));
296: /* reset the convergence flag from the previous solves */
297: eps->reason = EPS_CONVERGED_ITERATING;
299: /* Set default solver type (EPSSetFromOptions was not called) */
300: if (!((PetscObject)eps)->type_name) PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
301: if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
302: PetscCall(EPSSetDefaultST(eps));
304: PetscCall(STSetTransform(eps->st,PETSC_TRUE));
305: PetscCall(STSetStructured(eps->st,PETSC_FALSE));
306: if (eps->useds && !eps->ds) PetscCall(EPSGetDS(eps,&eps->ds));
307: if (eps->useds) PetscCall(EPSSetDSType(eps));
308: if (eps->twosided) {
309: PetscCheck(!eps->ishermitian || (eps->isgeneralized && !eps->ispositive),PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided methods are not intended for %s problems",SLEPC_STRING_HERMITIAN);
310: }
311: if (!eps->rg) PetscCall(EPSGetRG(eps,&eps->rg));
312: if (!((PetscObject)eps->rg)->type_name) PetscCall(RGSetType(eps->rg,RGINTERVAL));
314: /* Set problem dimensions */
315: PetscCall(STGetNumMatrices(eps->st,&nmat));
316: PetscCheck(nmat,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSSetOperators must be called first");
317: PetscCall(STMatGetSize(eps->st,&eps->n,NULL));
318: PetscCall(STMatGetLocalSize(eps->st,&eps->nloc,NULL));
320: /* Set default problem type */
321: if (!eps->problem_type) {
322: if (nmat==1) PetscCall(EPSSetProblemType(eps,EPS_NHEP));
323: else PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
324: } else if (nmat==1 && eps->isgeneralized) {
325: PetscCall(PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n"));
326: eps->isgeneralized = PETSC_FALSE;
327: eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;
328: } else PetscCheck(nmat==1 || eps->isgeneralized,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Inconsistent EPS state: the problem type does not match the number of matrices");
330: if (eps->isstructured) {
331: /* make sure the user has set the appropriate matrix */
332: PetscCall(STGetMatrix(eps->st,0,&A));
333: if (eps->problem_type==EPS_BSE) PetscCall(SlepcCheckMatStruct(A,SLEPC_MAT_STRUCT_BSE,NULL));
334: if (eps->problem_type==EPS_HAMILT) PetscCall(SlepcCheckMatStruct(A,SLEPC_MAT_STRUCT_HAMILT,NULL));
335: if (eps->problem_type==EPS_LREP) PetscCall(SlepcCheckMatStruct(A,SLEPC_MAT_STRUCT_LREP,NULL));
336: }
338: /* safeguard for small problems */
339: if (eps->n == 0) PetscFunctionReturn(PETSC_SUCCESS);
340: if (eps->nev > eps->n) eps->nev = eps->n;
341: if (eps->problem_type == EPS_BSE || eps->problem_type == EPS_LREP) {
342: if (2*eps->ncv > eps->n) eps->ncv = eps->n/2;
343: } else {
344: if (eps->ncv > eps->n) eps->ncv = eps->n;
345: }
347: /* check some combinations of eps->which */
348: PetscCheck(!eps->ishermitian || (eps->isgeneralized && !eps->ispositive) || (eps->which!=EPS_LARGEST_IMAGINARY && eps->which!=EPS_SMALLEST_IMAGINARY && eps->which!=EPS_TARGET_IMAGINARY),PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Sorting the eigenvalues along the imaginary axis is not allowed when all eigenvalues are real");
350: /* initialization of matrix norms */
351: if (eps->conv==EPS_CONV_NORM) {
352: if (!eps->nrma) {
353: PetscCall(STGetMatrix(eps->st,0,&A));
354: PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
355: }
356: if (nmat>1 && !eps->nrmb) {
357: PetscCall(STGetMatrix(eps->st,1,&B));
358: PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
359: }
360: }
362: /* call specific solver setup */
363: PetscUseTypeMethod(eps,setup);
365: /* threshold stopping test */
366: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STFILTER,&isfilter));
367: if (eps->stop==EPS_STOP_THRESHOLD && !isfilter) {
368: PetscCheck(eps->thres!=PETSC_MIN_REAL,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"Must give a threshold value with EPSSetThreshold()");
369: PetscCheck(eps->which==EPS_LARGEST_MAGNITUDE || eps->which==EPS_SMALLEST_MAGNITUDE || eps->which==EPS_LARGEST_REAL || eps->which==EPS_SMALLEST_REAL || eps->which==EPS_TARGET_MAGNITUDE,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Threshold stopping test can only be used with largest/smallest/target magnitude or largest/smallest real selection of eigenvalues");
370: if (eps->which==EPS_LARGEST_REAL || eps->which==EPS_SMALLEST_REAL) PetscCheck(eps->problem_type==EPS_HEP || eps->problem_type==EPS_GHEP || eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Threshold stopping test with largest/smallest real can only be used in problems that have all eigenvalues real");
371: else PetscCheck(eps->thres>0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"In case of largest/smallest/target magnitude the threshold value must be positive");
372: PetscCheck(eps->which==EPS_LARGEST_MAGNITUDE || eps->which==EPS_TARGET_MAGNITUDE || !eps->threlative,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"Can only use a relative threshold with largest/target magnitude selection of eigenvalues");
373: PetscCall(PetscNew(&ctx));
374: ctx->thres = eps->thres;
375: ctx->threlative = eps->threlative;
376: ctx->which = eps->which;
377: PetscCall(EPSSetStoppingTestFunction(eps,EPSStoppingThreshold,ctx,PetscCtxDestroyDefault));
378: }
380: /* if purification is set, check that it really makes sense */
381: if (eps->purify) {
382: if (eps->categ==EPS_CATEGORY_PRECOND || eps->categ==EPS_CATEGORY_CONTOUR) eps->purify = PETSC_FALSE;
383: else {
384: if (!eps->isgeneralized) eps->purify = PETSC_FALSE;
385: else if (!eps->ishermitian && !eps->ispositive) eps->purify = PETSC_FALSE;
386: else {
387: PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg));
388: if (flg) eps->purify = PETSC_FALSE;
389: }
390: }
391: }
393: /* set tolerance if not yet set */
394: if (eps->tol==(PetscReal)PETSC_DETERMINE) eps->tol = SLEPC_DEFAULT_TOL;
396: /* set up sorting criterion */
397: PetscTryTypeMethod(eps,setupsort);
399: /* Build balancing matrix if required */
400: if (eps->balance!=EPS_BALANCE_USER) {
401: PetscCall(STSetBalanceMatrix(eps->st,NULL));
402: if (!eps->ishermitian && (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE)) {
403: if (!eps->D) PetscCall(BVCreateVec(eps->V,&eps->D));
404: PetscCall(EPSBuildBalance_Krylov(eps));
405: PetscCall(STSetBalanceMatrix(eps->st,eps->D));
406: }
407: }
409: /* Setup ST */
410: PetscCall(STSetUp(eps->st));
411: PetscCall(EPSCheckCompatibleST(eps));
413: /* threshold stopping test for STFILTER */
414: if (isfilter) {
415: PetscCall(PetscNew(&ctx));
416: PetscCall(STFilterGetThreshold(eps->st,&eps->thres));
417: eps->threlative = PETSC_FALSE;
418: ctx->thres = eps->thres;
419: ctx->threlative = eps->threlative;
420: ctx->which = EPS_LARGEST_MAGNITUDE;
421: ctx->its = -1;
422: PetscCall(EPSSetStoppingTestFunction(eps,EPSStoppingThreshold,ctx,PetscCtxDestroyDefault));
423: }
425: /* process deflation and initial vectors */
426: if (eps->nds<0) {
427: PetscCheck(!eps->isstructured,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Deflation space is not supported in structured eigenproblems");
428: k = -eps->nds;
429: PetscCall(BVInsertConstraints(eps->V,&k,eps->defl));
430: PetscCall(SlepcBasisDestroy_Private(&eps->nds,&eps->defl));
431: eps->nds = k;
432: PetscCall(STCheckNullSpace(eps->st,eps->V));
433: }
434: if (eps->nini<0) {
435: k = -eps->nini;
436: PetscCheck(k<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
437: PetscCall(BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE));
438: PetscCall(SlepcBasisDestroy_Private(&eps->nini,&eps->IS));
439: eps->nini = k;
440: }
441: if (eps->twosided && eps->ninil<0) {
442: k = -eps->ninil;
443: PetscCheck(k<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
444: PetscCall(BVInsertVecs(eps->W,0,&k,eps->ISL,PETSC_TRUE));
445: PetscCall(SlepcBasisDestroy_Private(&eps->ninil,&eps->ISL));
446: eps->ninil = k;
447: }
449: PetscCall(PetscLogEventEnd(EPS_SetUp,eps,0,0,0));
450: eps->state = EPS_STATE_SETUP;
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
454: /*@
455: EPSSetOperators - Sets the matrices associated with the eigenvalue problem.
457: Collective
459: Input Parameters:
460: + eps - the linear eigensolver context
461: . A - the matrix associated with the eigensystem
462: - B - the second matrix in the case of generalized eigenproblems
464: Notes:
465: To specify a standard eigenproblem, use `NULL` for parameter `B`.
467: It must be called before `EPSSetUp()`. If it is called again after `EPSSetUp()` and
468: the matrix sizes have changed then the `EPS` object is reset.
470: For structured eigenproblem types such as `EPS_BSE` (see `EPSSetProblemType()`), the
471: provided matrices must have been created with the corresponding helper function,
472: i.e., `MatCreateBSE()`.
474: Level: beginner
476: .seealso: [](ch:eps), `EPSSolve()`, `EPSSetUp()`, `EPSReset()`, `EPSSetProblemType()`
477: @*/
478: PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)
479: {
480: PetscInt m,n,m0,mloc,nloc,mloc0,nmat;
481: Mat mat[2];
482: VecType ta,tb;
483: PetscBool same;
485: PetscFunctionBegin;
489: PetscCheckSameComm(eps,1,A,2);
490: if (B) PetscCheckSameComm(eps,1,B,3);
492: /* Check matrix sizes */
493: PetscCall(MatGetSize(A,&m,&n));
494: PetscCall(MatGetLocalSize(A,&mloc,&nloc));
495: PetscCheck(m==n,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
496: PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A does not have equal row and column sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",mloc,nloc);
497: if (B) {
498: PetscCall(MatGetSize(B,&m0,&n));
499: PetscCall(MatGetLocalSize(B,&mloc0,&nloc));
500: PetscCheck(m0==n,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m0,n);
501: PetscCheck(mloc0==nloc,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",mloc0,nloc);
502: PetscCheck(m==m0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Dimensions of A and B do not match (%" PetscInt_FMT ", %" PetscInt_FMT ")",m,m0);
503: PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Local dimensions of A and B do not match (%" PetscInt_FMT ", %" PetscInt_FMT ")",mloc,mloc0);
504: /* make sure both matrices have compatible VecType */
505: PetscCall(MatGetVecType(A,&ta));
506: PetscCall(MatGetVecType(B,&tb));
507: PetscCall(PetscStrcmp(ta,tb,&same));
508: PetscCheck(same,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"The provided matrices have different vector types (%s vs %s), consider calling MatSetVecType() in the one that is not on GPU",ta,tb);
509: }
510: if (eps->state && (n!=eps->n || nloc!=eps->nloc)) PetscCall(EPSReset(eps));
511: eps->nrma = 0.0;
512: eps->nrmb = 0.0;
513: if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
514: mat[0] = A;
515: if (B) {
516: mat[1] = B;
517: nmat = 2;
518: } else nmat = 1;
519: PetscCall(STSetMatrices(eps->st,nmat,mat));
520: eps->state = EPS_STATE_INITIAL;
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /*@
525: EPSGetOperators - Gets the matrices associated with the eigensystem.
527: Collective
529: Input Parameter:
530: . eps - the linear eigensolver context
532: Output Parameters:
533: + A - the matrix associated with the eigensystem
534: - B - the second matrix in the case of generalized eigenproblems
536: Note:
537: Does not increase the reference count of the matrices, so you should not destroy them.
539: Level: intermediate
541: .seealso: [](ch:eps), `EPSSolve()`, `EPSGetST()`, `STGetMatrix()`, `STSetMatrices()`
542: @*/
543: PetscErrorCode EPSGetOperators(EPS eps,Mat *A,Mat *B)
544: {
545: ST st;
546: PetscInt k;
548: PetscFunctionBegin;
550: PetscCall(EPSGetST(eps,&st));
551: PetscCall(STGetNumMatrices(st,&k));
552: if (A) {
553: if (k<1) *A = NULL;
554: else PetscCall(STGetMatrix(st,0,A));
555: }
556: if (B) {
557: if (k<2) *B = NULL;
558: else PetscCall(STGetMatrix(st,1,B));
559: }
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: /*@
564: EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
565: space.
567: Collective
569: Input Parameters:
570: + eps - the linear eigensolver context
571: . n - number of vectors
572: - v - set of basis vectors of the deflation space
574: Notes:
575: When a deflation space is given, the eigensolver seeks the eigensolution
576: in the restriction of the problem to the orthogonal complement of this
577: space. This can be used for instance in the case that an invariant
578: subspace is known beforehand (such as the nullspace of the matrix).
580: These vectors do not persist from one `EPSSolve()` call to the other, so the
581: deflation space should be set every time.
583: The vectors do not need to be mutually orthonormal, since they are explicitly
584: orthonormalized internally.
586: Level: intermediate
588: .seealso: [](ch:eps), `EPSSetInitialSpace()`
589: @*/
590: PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec v[])
591: {
592: PetscFunctionBegin;
595: PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
596: if (n>0) {
597: PetscAssertPointer(v,3);
599: }
600: PetscCall(SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl));
601: if (n>0) eps->state = EPS_STATE_INITIAL;
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: /*@
606: EPSSetInitialSpace - Specify a basis of vectors that constitute the initial
607: space, that is, the subspace from which the solver starts to iterate.
609: Collective
611: Input Parameters:
612: + eps - the linear eigensolver context
613: . n - number of vectors
614: - is - set of basis vectors of the initial space
616: Notes:
617: Some solvers such as `EPSKRYLOVSCHUR` start to iterate on a single vector
618: (initial vector). In that case, only `is[0]` is taken into account and the
619: other vectors are ignored. But other solvers such as `EPSSUBSPACE` are
620: able to make use of the whole initial subspace as an initial guess.
622: These vectors do not persist from one `EPSSolve()` call to the other, so the
623: initial space should be set every time.
625: The vectors do not need to be mutually orthonormal, since they are explicitly
626: orthonormalized internally.
628: Common usage of this function is when the user can provide a rough approximation
629: of the wanted eigenspace. Then, convergence may be faster.
631: Level: intermediate
633: .seealso: [](ch:eps), `EPSSetLeftInitialSpace()`, `EPSSetDeflationSpace()`
634: @*/
635: PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[])
636: {
637: PetscFunctionBegin;
640: PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
641: if (n>0) {
642: PetscAssertPointer(is,3);
644: }
645: PetscCall(SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS));
646: if (n>0) eps->state = EPS_STATE_INITIAL;
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: /*@
651: EPSSetLeftInitialSpace - Specify a basis of vectors that constitute the left
652: initial space, used by two-sided solvers to start the left subspace.
654: Collective
656: Input Parameters:
657: + eps - the linear eigensolver context
658: . n - number of vectors
659: - isl - set of basis vectors of the left initial space
661: Notes:
662: Left initial vectors are used to initiate the left search space in two-sided
663: eigensolvers. Users should pass here an approximation of the left eigenspace,
664: if available.
666: The same comments in `EPSSetInitialSpace()` are applicable here.
668: Level: intermediate
670: .seealso: [](ch:eps), `EPSSetInitialSpace()`, `EPSSetTwoSided()`
671: @*/
672: PetscErrorCode EPSSetLeftInitialSpace(EPS eps,PetscInt n,Vec isl[])
673: {
674: PetscFunctionBegin;
677: PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
678: if (n>0) {
679: PetscAssertPointer(isl,3);
681: }
682: PetscCall(SlepcBasisReference_Private(n,isl,&eps->ninil,&eps->ISL));
683: if (n>0) eps->state = EPS_STATE_INITIAL;
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: /*
688: EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
689: by the user. This is called at setup.
690: */
691: PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt *nev,PetscInt *ncv,PetscInt *mpd)
692: {
693: PetscBool krylov;
694: PetscInt nev2, n = (eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP)? eps->n/2: eps->n;
696: PetscFunctionBegin;
697: if (*nev==0 && eps->stop!=EPS_STOP_THRESHOLD) *nev = 1;
698: nev2 = (eps->problem_type==EPS_BSE || eps->problem_type==EPS_LREP)? (*nev+1)/2: *nev;
699: if (*ncv!=PETSC_DETERMINE) { /* ncv set */
700: PetscCall(PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,""));
701: if (krylov) {
702: PetscCheck(*ncv>=nev2+1 || (*ncv==nev2 && *ncv==n),PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
703: } else {
704: PetscCheck(*ncv>=nev2,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
705: }
706: } else if (*mpd!=PETSC_DETERMINE) { /* mpd set */
707: *ncv = PetscMin(n,nev2+(*mpd));
708: } else { /* neither set: defaults depend on nev being small or large */
709: if (nev2<500) *ncv = PetscMin(n,PetscMax(2*(nev2),nev2+15));
710: else {
711: *mpd = 500;
712: *ncv = PetscMin(n,nev2+(*mpd));
713: }
714: }
715: if (*mpd==PETSC_DETERMINE) *mpd = *ncv;
716: PetscFunctionReturn(PETSC_SUCCESS);
717: }
719: /*@
720: EPSAllocateSolution - Allocate memory storage for common variables such
721: as eigenvalues and eigenvectors.
723: Collective
725: Input Parameters:
726: + eps - the linear eigensolver context
727: - extra - number of additional positions, used for methods that require a
728: working basis slightly larger than `ncv`
730: Developer Note:
731: This is `SLEPC_EXTERN` because it may be required by user plugin `EPS`
732: implementations.
734: Level: developer
736: .seealso: [](ch:eps), `EPSSetUp()`, `EPSSetDimensions()`
737: @*/
738: PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
739: {
740: PetscInt oldsize,requested;
741: PetscRandom rand;
742: Vec t;
744: PetscFunctionBegin;
745: requested = eps->ncv + extra;
747: /* oldsize is zero if this is the first time setup is called */
748: PetscCall(BVGetSizes(eps->V,NULL,NULL,&oldsize));
750: /* allocate space for eigenvalues and friends */
751: if (requested != oldsize || !eps->eigr) {
752: PetscCall(PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm));
753: PetscCall(PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm));
754: }
756: /* workspace for the case of arbitrary selection */
757: if (eps->arbitrary) {
758: if (eps->rr) PetscCall(PetscFree2(eps->rr,eps->ri));
759: PetscCall(PetscMalloc2(requested,&eps->rr,requested,&eps->ri));
760: }
762: /* allocate V */
763: if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
764: if (!oldsize) {
765: if (!((PetscObject)eps->V)->type_name) PetscCall(BVSetType(eps->V,BVMAT));
766: PetscCall(STMatCreateVecsEmpty(eps->st,&t,NULL));
767: PetscCall(BVSetSizesFromVec(eps->V,t,requested));
768: PetscCall(VecDestroy(&t));
769: } else PetscCall(BVResize(eps->V,requested,PETSC_FALSE));
771: /* allocate W */
772: if (eps->twosided) {
773: PetscCall(BVGetRandomContext(eps->V,&rand)); /* make sure the random context is available when duplicating */
774: PetscCall(BVDestroy(&eps->W));
775: PetscCall(BVDuplicate(eps->V,&eps->W));
776: }
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: /*@
781: EPSReallocateSolution - Reallocate memory storage for common variables such
782: as the eigenvalues and the basis vectors.
784: Collective
786: Input Parameters:
787: + eps - the linear eigensolver context
788: - newsize - new size
790: Developer Notes:
791: This is `SLEPC_EXTERN` because it may be required by user plugin `EPS`
792: implementations.
794: This is called during the iteration in case the threshold stopping test has
795: been selected.
797: Level: developer
799: .seealso: [](ch:eps), `EPSAllocateSolution()`, `EPSSetThreshold()`
800: @*/
801: PetscErrorCode EPSReallocateSolution(EPS eps,PetscInt newsize)
802: {
803: PetscInt oldsize,*nperm;
804: PetscReal *nerrest;
805: PetscScalar *neigr,*neigi;
807: PetscFunctionBegin;
808: PetscCall(BVGetSizes(eps->V,NULL,NULL,&oldsize));
809: if (oldsize>=newsize) PetscFunctionReturn(PETSC_SUCCESS);
810: PetscCall(PetscInfo(eps,"Reallocating basis vectors to %" PetscInt_FMT " columns\n",newsize));
812: /* reallocate eigenvalues */
813: PetscCall(PetscMalloc4(newsize,&neigr,newsize,&neigi,newsize,&nerrest,newsize,&nperm));
814: PetscCall(PetscArraycpy(neigr,eps->eigr,oldsize));
815: PetscCall(PetscArraycpy(neigi,eps->eigi,oldsize));
816: PetscCall(PetscArraycpy(nerrest,eps->errest,oldsize));
817: PetscCall(PetscArraycpy(nperm,eps->perm,oldsize));
818: PetscCall(PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm));
819: eps->eigr = neigr;
820: eps->eigi = neigi;
821: eps->errest = nerrest;
822: eps->perm = nperm;
823: /* reallocate V,W */
824: PetscCall(BVResize(eps->V,newsize,PETSC_TRUE));
825: if (eps->twosided) PetscCall(BVResize(eps->W,newsize,PETSC_TRUE));
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }