Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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 : */
13 :
14 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
15 :
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 1495 : PetscErrorCode EPSSetDefaultST(EPS eps)
23 : {
24 1495 : PetscFunctionBegin;
25 1495 : PetscTryTypeMethod(eps,setdefaultst);
26 1495 : if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSHIFT));
27 1495 : PetscFunctionReturn(PETSC_SUCCESS);
28 : }
29 :
30 : /*
31 : This is done by preconditioned eigensolvers that use the PC only.
32 : It sets STPRECOND with KSPPREONLY.
33 : */
34 103 : PetscErrorCode EPSSetDefaultST_Precond(EPS eps)
35 : {
36 103 : KSP ksp;
37 :
38 103 : PetscFunctionBegin;
39 103 : if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STPRECOND));
40 103 : PetscCall(STGetKSP(eps->st,&ksp));
41 103 : if (!((PetscObject)ksp)->type_name) PetscCall(KSPSetType(ksp,KSPPREONLY));
42 103 : PetscFunctionReturn(PETSC_SUCCESS);
43 : }
44 :
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 107 : PetscErrorCode EPSSetDefaultST_GMRES(EPS eps)
50 : {
51 107 : KSP ksp;
52 :
53 107 : PetscFunctionBegin;
54 107 : if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STPRECOND));
55 107 : PetscCall(STPrecondSetKSPHasMat(eps->st,PETSC_TRUE));
56 107 : PetscCall(STGetKSP(eps->st,&ksp));
57 107 : if (!((PetscObject)ksp)->type_name) {
58 43 : PetscCall(KSPSetType(ksp,KSPGMRES));
59 43 : PetscCall(KSPSetTolerances(ksp,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT,5));
60 : }
61 107 : PetscFunctionReturn(PETSC_SUCCESS);
62 : }
63 :
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 24 : PetscErrorCode EPSSetDefaultST_NoFactor(EPS eps)
70 : {
71 24 : KSP ksp;
72 24 : PC pc;
73 :
74 24 : PetscFunctionBegin;
75 24 : if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSHIFT));
76 24 : PetscCall(STGetKSP(eps->st,&ksp));
77 24 : if (!((PetscObject)ksp)->type_name) PetscCall(KSPSetType(ksp,KSPPREONLY));
78 24 : PetscCall(KSPGetPC(ksp,&pc));
79 24 : if (!((PetscObject)pc)->type_name) PetscCall(PCSetType(pc,PCNONE));
80 24 : PetscFunctionReturn(PETSC_SUCCESS);
81 : }
82 : #endif
83 :
84 : /*
85 : Check that the ST selected by the user is compatible with the EPS solver and options
86 : */
87 887 : static PetscErrorCode EPSCheckCompatibleST(EPS eps)
88 : {
89 887 : PetscBool precond,shift,sinvert,cayley,lyapii;
90 : #if defined(PETSC_USE_COMPLEX)
91 : PetscScalar sigma;
92 : #endif
93 :
94 887 : PetscFunctionBegin;
95 887 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STPRECOND,&precond));
96 887 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&shift));
97 887 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert));
98 887 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&cayley));
99 :
100 : /* preconditioned eigensolvers */
101 887 : PetscCheck(eps->categ!=EPS_CATEGORY_PRECOND || precond,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires ST=PRECOND");
102 887 : PetscCheck(eps->categ==EPS_CATEGORY_PRECOND || !precond,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"STPRECOND is intended for preconditioned eigensolvers only");
103 :
104 : /* harmonic extraction */
105 887 : PetscCheck(precond || shift || !eps->extraction || eps->extraction==EPS_RITZ,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot use a spectral transformation combined with harmonic extraction");
106 :
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
112 :
113 : /* Cayley with PGNHEP */
114 887 : PetscCheck(!cayley || eps->problem_type!=EPS_PGNHEP,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cayley spectral transformation is not compatible with PGNHEP");
115 :
116 : /* make sure that the user does not specify smallest magnitude with shift-and-invert */
117 887 : if ((cayley || sinvert) && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER)) {
118 122 : PetscCall(PetscObjectTypeCompare((PetscObject)eps,EPSLYAPII,&lyapii));
119 122 : 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 887 : PetscFunctionReturn(PETSC_SUCCESS);
122 : }
123 :
124 : /*
125 : MatEstimateSpectralRange_EPS: estimate the spectral range [left,right] of a
126 : symmetric/Hermitian matrix A using an auxiliary EPS object
127 : */
128 7 : PetscErrorCode MatEstimateSpectralRange_EPS(Mat A,PetscReal *left,PetscReal *right)
129 : {
130 7 : PetscInt nconv;
131 7 : PetscScalar eig0;
132 7 : PetscReal tol=1e-3,errest=tol;
133 7 : EPS eps;
134 :
135 7 : PetscFunctionBegin;
136 7 : *left = 0.0; *right = 0.0;
137 7 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)A),&eps));
138 7 : PetscCall(EPSSetOptionsPrefix(eps,"eps_filter_"));
139 7 : PetscCall(EPSSetOperators(eps,A,NULL));
140 7 : PetscCall(EPSSetProblemType(eps,EPS_HEP));
141 7 : PetscCall(EPSSetTolerances(eps,tol,50));
142 7 : PetscCall(EPSSetConvergenceTest(eps,EPS_CONV_ABS));
143 7 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
144 7 : PetscCall(EPSSolve(eps));
145 7 : PetscCall(EPSGetConverged(eps,&nconv));
146 7 : if (nconv>0) {
147 7 : PetscCall(EPSGetEigenvalue(eps,0,&eig0,NULL));
148 7 : PetscCall(EPSGetErrorEstimate(eps,0,&errest));
149 0 : } else eig0 = eps->eigr[0];
150 7 : *left = PetscRealPart(eig0)-errest;
151 7 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
152 7 : PetscCall(EPSSolve(eps));
153 7 : PetscCall(EPSGetConverged(eps,&nconv));
154 7 : if (nconv>0) {
155 7 : PetscCall(EPSGetEigenvalue(eps,0,&eig0,NULL));
156 7 : PetscCall(EPSGetErrorEstimate(eps,0,&errest));
157 0 : } else eig0 = eps->eigr[0];
158 7 : *right = PetscRealPart(eig0)+errest;
159 7 : PetscCall(EPSDestroy(&eps));
160 7 : PetscFunctionReturn(PETSC_SUCCESS);
161 : }
162 :
163 : /*
164 : EPSSetUpSort_Basic: configure the EPS sorting criterion according to 'which'
165 : */
166 826 : PetscErrorCode EPSSetUpSort_Basic(EPS eps)
167 : {
168 826 : PetscFunctionBegin;
169 826 : switch (eps->which) {
170 299 : case EPS_LARGEST_MAGNITUDE:
171 299 : eps->sc->comparison = SlepcCompareLargestMagnitude;
172 299 : eps->sc->comparisonctx = NULL;
173 299 : break;
174 14 : case EPS_SMALLEST_MAGNITUDE:
175 14 : eps->sc->comparison = SlepcCompareSmallestMagnitude;
176 14 : eps->sc->comparisonctx = NULL;
177 14 : break;
178 148 : case EPS_LARGEST_REAL:
179 148 : eps->sc->comparison = SlepcCompareLargestReal;
180 148 : eps->sc->comparisonctx = NULL;
181 148 : break;
182 159 : case EPS_SMALLEST_REAL:
183 159 : eps->sc->comparison = SlepcCompareSmallestReal;
184 159 : eps->sc->comparisonctx = NULL;
185 159 : break;
186 1 : case EPS_LARGEST_IMAGINARY:
187 1 : eps->sc->comparison = SlepcCompareLargestImaginary;
188 1 : eps->sc->comparisonctx = NULL;
189 1 : break;
190 1 : case EPS_SMALLEST_IMAGINARY:
191 1 : eps->sc->comparison = SlepcCompareSmallestImaginary;
192 1 : eps->sc->comparisonctx = NULL;
193 1 : break;
194 117 : case EPS_TARGET_MAGNITUDE:
195 117 : eps->sc->comparison = SlepcCompareTargetMagnitude;
196 117 : eps->sc->comparisonctx = &eps->target;
197 117 : break;
198 6 : case EPS_TARGET_REAL:
199 6 : eps->sc->comparison = SlepcCompareTargetReal;
200 6 : eps->sc->comparisonctx = &eps->target;
201 6 : break;
202 : case EPS_TARGET_IMAGINARY:
203 : #if defined(PETSC_USE_COMPLEX)
204 : eps->sc->comparison = SlepcCompareTargetImaginary;
205 : eps->sc->comparisonctx = &eps->target;
206 : #endif
207 : break;
208 60 : case EPS_ALL:
209 60 : eps->sc->comparison = SlepcCompareSmallestReal;
210 60 : eps->sc->comparisonctx = NULL;
211 60 : break;
212 : case EPS_WHICH_USER:
213 : break;
214 : }
215 826 : eps->sc->map = NULL;
216 826 : eps->sc->mapobj = NULL;
217 826 : PetscFunctionReturn(PETSC_SUCCESS);
218 : }
219 :
220 : /*
221 : EPSSetUpSort_Default: configure both EPS and DS sorting criterion
222 : */
223 779 : PetscErrorCode EPSSetUpSort_Default(EPS eps)
224 : {
225 779 : SlepcSC sc;
226 779 : PetscBool istrivial;
227 :
228 779 : PetscFunctionBegin;
229 : /* fill sorting criterion context */
230 779 : PetscCall(EPSSetUpSort_Basic(eps));
231 : /* fill sorting criterion for DS */
232 779 : PetscCall(DSGetSlepcSC(eps->ds,&sc));
233 779 : PetscCall(RGIsTrivial(eps->rg,&istrivial));
234 779 : sc->rg = istrivial? NULL: eps->rg;
235 779 : sc->comparison = eps->sc->comparison;
236 779 : sc->comparisonctx = eps->sc->comparisonctx;
237 779 : sc->map = SlepcMap_ST;
238 779 : sc->mapobj = (PetscObject)eps->st;
239 779 : PetscFunctionReturn(PETSC_SUCCESS);
240 : }
241 :
242 : /*@
243 : EPSSetDSType - Sets the type of the internal DS object based on the current
244 : settings of the eigenvalue solver.
245 :
246 : Collective
247 :
248 : Input Parameter:
249 : . eps - eigenproblem solver context
250 :
251 : Note:
252 : This function need not be called explicitly, since it will be called at
253 : both EPSSetFromOptions() and EPSSetUp().
254 :
255 : Level: developer
256 :
257 : .seealso: EPSSetFromOptions(), EPSSetUp()
258 : @*/
259 1328 : PetscErrorCode EPSSetDSType(EPS eps)
260 : {
261 1328 : PetscFunctionBegin;
262 1328 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
263 1328 : PetscTryTypeMethod(eps,setdstype);
264 1328 : PetscFunctionReturn(PETSC_SUCCESS);
265 : }
266 :
267 : /*@
268 : EPSSetUp - Sets up all the internal data structures necessary for the
269 : execution of the eigensolver. Then calls STSetUp() for any set-up
270 : operations associated to the ST object.
271 :
272 : Collective
273 :
274 : Input Parameter:
275 : . eps - eigenproblem solver context
276 :
277 : Notes:
278 : This function need not be called explicitly in most cases, since EPSSolve()
279 : calls it. It can be useful when one wants to measure the set-up time
280 : separately from the solve time.
281 :
282 : Level: developer
283 :
284 : .seealso: EPSCreate(), EPSSolve(), EPSDestroy(), STSetUp(), EPSSetInitialSpace()
285 : @*/
286 1062 : PetscErrorCode EPSSetUp(EPS eps)
287 : {
288 1062 : Mat A,B;
289 1062 : PetscInt k,nmat;
290 1062 : PetscBool flg;
291 :
292 1062 : PetscFunctionBegin;
293 1062 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
294 1062 : if (eps->state) PetscFunctionReturn(PETSC_SUCCESS);
295 887 : PetscCall(PetscLogEventBegin(EPS_SetUp,eps,0,0,0));
296 :
297 : /* reset the convergence flag from the previous solves */
298 887 : eps->reason = EPS_CONVERGED_ITERATING;
299 :
300 : /* Set default solver type (EPSSetFromOptions was not called) */
301 887 : if (!((PetscObject)eps)->type_name) PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
302 887 : if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
303 887 : PetscCall(EPSSetDefaultST(eps));
304 :
305 887 : PetscCall(STSetTransform(eps->st,PETSC_TRUE));
306 887 : PetscCall(STSetStructured(eps->st,PETSC_FALSE));
307 887 : if (eps->useds && !eps->ds) PetscCall(EPSGetDS(eps,&eps->ds));
308 887 : if (eps->useds) PetscCall(EPSSetDSType(eps));
309 887 : if (eps->twosided) {
310 19 : 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);
311 : }
312 887 : if (!eps->rg) PetscCall(EPSGetRG(eps,&eps->rg));
313 887 : if (!((PetscObject)eps->rg)->type_name) PetscCall(RGSetType(eps->rg,RGINTERVAL));
314 :
315 : /* Set problem dimensions */
316 887 : PetscCall(STGetNumMatrices(eps->st,&nmat));
317 887 : PetscCheck(nmat,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSSetOperators must be called first");
318 887 : PetscCall(STMatGetSize(eps->st,&eps->n,NULL));
319 887 : PetscCall(STMatGetLocalSize(eps->st,&eps->nloc,NULL));
320 :
321 : /* Set default problem type */
322 887 : if (!eps->problem_type) {
323 41 : if (nmat==1) PetscCall(EPSSetProblemType(eps,EPS_NHEP));
324 5 : else PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
325 846 : } else if (nmat==1 && eps->isgeneralized) {
326 0 : PetscCall(PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n"));
327 0 : eps->isgeneralized = PETSC_FALSE;
328 0 : eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;
329 846 : } 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 :
331 887 : if (eps->isstructured) {
332 : /* make sure the user has set the appropriate matrix */
333 14 : PetscCall(STGetMatrix(eps->st,0,&A));
334 14 : if (eps->problem_type==EPS_BSE) PetscCall(SlepcCheckMatStruct(A,SLEPC_MAT_STRUCT_BSE,NULL));
335 : }
336 :
337 887 : if (eps->nev > eps->n) eps->nev = eps->n;
338 887 : if (eps->ncv > eps->n) eps->ncv = eps->n;
339 :
340 : /* check some combinations of eps->which */
341 887 : 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");
342 :
343 : /* initialization of matrix norms */
344 887 : if (eps->conv==EPS_CONV_NORM) {
345 43 : if (!eps->nrma) {
346 43 : PetscCall(STGetMatrix(eps->st,0,&A));
347 43 : PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
348 : }
349 43 : if (nmat>1 && !eps->nrmb) {
350 42 : PetscCall(STGetMatrix(eps->st,1,&B));
351 42 : PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
352 : }
353 : }
354 :
355 : /* call specific solver setup */
356 887 : PetscUseTypeMethod(eps,setup);
357 :
358 : /* if purification is set, check that it really makes sense */
359 887 : if (eps->purify) {
360 641 : if (eps->categ==EPS_CATEGORY_PRECOND || eps->categ==EPS_CATEGORY_CONTOUR) eps->purify = PETSC_FALSE;
361 : else {
362 470 : if (!eps->isgeneralized) eps->purify = PETSC_FALSE;
363 147 : else if (!eps->ishermitian && !eps->ispositive) eps->purify = PETSC_FALSE;
364 : else {
365 120 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg));
366 120 : if (flg) eps->purify = PETSC_FALSE;
367 : }
368 : }
369 : }
370 :
371 : /* set tolerance if not yet set */
372 887 : if (eps->tol==(PetscReal)PETSC_DETERMINE) eps->tol = SLEPC_DEFAULT_TOL;
373 :
374 : /* set up sorting criterion */
375 887 : PetscTryTypeMethod(eps,setupsort);
376 :
377 : /* Build balancing matrix if required */
378 887 : if (eps->balance!=EPS_BALANCE_USER) {
379 887 : PetscCall(STSetBalanceMatrix(eps->st,NULL));
380 887 : if (!eps->ishermitian && (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE)) {
381 13 : if (!eps->D) PetscCall(BVCreateVec(eps->V,&eps->D));
382 13 : PetscCall(EPSBuildBalance_Krylov(eps));
383 13 : PetscCall(STSetBalanceMatrix(eps->st,eps->D));
384 : }
385 : }
386 :
387 : /* Setup ST */
388 887 : PetscCall(STSetUp(eps->st));
389 887 : PetscCall(EPSCheckCompatibleST(eps));
390 :
391 : /* process deflation and initial vectors */
392 887 : if (eps->nds<0) {
393 14 : PetscCheck(!eps->isstructured,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Deflation space is not supported in structured eigenproblems");
394 14 : k = -eps->nds;
395 14 : PetscCall(BVInsertConstraints(eps->V,&k,eps->defl));
396 14 : PetscCall(SlepcBasisDestroy_Private(&eps->nds,&eps->defl));
397 14 : eps->nds = k;
398 14 : PetscCall(STCheckNullSpace(eps->st,eps->V));
399 : }
400 887 : if (eps->nini<0) {
401 224 : k = -eps->nini;
402 224 : PetscCheck(k<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
403 224 : PetscCall(BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE));
404 224 : PetscCall(SlepcBasisDestroy_Private(&eps->nini,&eps->IS));
405 224 : eps->nini = k;
406 : }
407 887 : if (eps->twosided && eps->ninil<0) {
408 4 : k = -eps->ninil;
409 4 : PetscCheck(k<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
410 4 : PetscCall(BVInsertVecs(eps->W,0,&k,eps->ISL,PETSC_TRUE));
411 4 : PetscCall(SlepcBasisDestroy_Private(&eps->ninil,&eps->ISL));
412 4 : eps->ninil = k;
413 : }
414 :
415 887 : PetscCall(PetscLogEventEnd(EPS_SetUp,eps,0,0,0));
416 887 : eps->state = EPS_STATE_SETUP;
417 887 : PetscFunctionReturn(PETSC_SUCCESS);
418 : }
419 :
420 : /*@
421 : EPSSetOperators - Sets the matrices associated with the eigenvalue problem.
422 :
423 : Collective
424 :
425 : Input Parameters:
426 : + eps - the eigenproblem solver context
427 : . A - the matrix associated with the eigensystem
428 : - B - the second matrix in the case of generalized eigenproblems
429 :
430 : Notes:
431 : To specify a standard eigenproblem, use NULL for parameter B.
432 :
433 : It must be called before EPSSetUp(). If it is called again after EPSSetUp() and
434 : the matrix sizes have changed then the EPS object is reset.
435 :
436 : For structured eigenproblem types such as EPS_BSE (see EPSSetProblemType()), the
437 : provided matrices must have been created with the corresponding helper function,
438 : i.e., MatCreateBSE().
439 :
440 : Level: beginner
441 :
442 : .seealso: EPSSolve(), EPSSetUp(), EPSReset(), EPSGetST(), STGetMatrix(), EPSSetProblemType()
443 : @*/
444 750 : PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)
445 : {
446 750 : PetscInt m,n,m0,mloc,nloc,mloc0,nmat;
447 750 : Mat mat[2];
448 :
449 750 : PetscFunctionBegin;
450 750 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
451 750 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
452 750 : if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3);
453 750 : PetscCheckSameComm(eps,1,A,2);
454 750 : if (B) PetscCheckSameComm(eps,1,B,3);
455 :
456 : /* Check matrix sizes */
457 750 : PetscCall(MatGetSize(A,&m,&n));
458 750 : PetscCall(MatGetLocalSize(A,&mloc,&nloc));
459 750 : PetscCheck(m==n,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
460 750 : 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);
461 750 : if (B) {
462 231 : PetscCall(MatGetSize(B,&m0,&n));
463 231 : PetscCall(MatGetLocalSize(B,&mloc0,&nloc));
464 231 : PetscCheck(m0==n,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m0,n);
465 231 : 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);
466 231 : PetscCheck(m==m0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Dimensions of A and B do not match (%" PetscInt_FMT ", %" PetscInt_FMT ")",m,m0);
467 231 : 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);
468 : }
469 750 : if (eps->state && (n!=eps->n || nloc!=eps->nloc)) PetscCall(EPSReset(eps));
470 750 : eps->nrma = 0.0;
471 750 : eps->nrmb = 0.0;
472 750 : if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
473 750 : mat[0] = A;
474 750 : if (B) {
475 231 : mat[1] = B;
476 231 : nmat = 2;
477 : } else nmat = 1;
478 750 : PetscCall(STSetMatrices(eps->st,nmat,mat));
479 750 : eps->state = EPS_STATE_INITIAL;
480 750 : PetscFunctionReturn(PETSC_SUCCESS);
481 : }
482 :
483 : /*@
484 : EPSGetOperators - Gets the matrices associated with the eigensystem.
485 :
486 : Collective
487 :
488 : Input Parameter:
489 : . eps - the EPS context
490 :
491 : Output Parameters:
492 : + A - the matrix associated with the eigensystem
493 : - B - the second matrix in the case of generalized eigenproblems
494 :
495 : Note:
496 : Does not increase the reference count of the matrices, so you should not destroy them.
497 :
498 : Level: intermediate
499 :
500 : .seealso: EPSSolve(), EPSGetST(), STGetMatrix(), STSetMatrices()
501 : @*/
502 887 : PetscErrorCode EPSGetOperators(EPS eps,Mat *A,Mat *B)
503 : {
504 887 : ST st;
505 887 : PetscInt k;
506 :
507 887 : PetscFunctionBegin;
508 887 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
509 887 : PetscCall(EPSGetST(eps,&st));
510 887 : PetscCall(STGetNumMatrices(st,&k));
511 887 : if (A) {
512 450 : if (k<1) *A = NULL;
513 450 : else PetscCall(STGetMatrix(st,0,A));
514 : }
515 887 : if (B) {
516 594 : if (k<2) *B = NULL;
517 565 : else PetscCall(STGetMatrix(st,1,B));
518 : }
519 887 : PetscFunctionReturn(PETSC_SUCCESS);
520 : }
521 :
522 : /*@
523 : EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
524 : space.
525 :
526 : Collective
527 :
528 : Input Parameters:
529 : + eps - the eigenproblem solver context
530 : . n - number of vectors
531 : - v - set of basis vectors of the deflation space
532 :
533 : Notes:
534 : When a deflation space is given, the eigensolver seeks the eigensolution
535 : in the restriction of the problem to the orthogonal complement of this
536 : space. This can be used for instance in the case that an invariant
537 : subspace is known beforehand (such as the nullspace of the matrix).
538 :
539 : These vectors do not persist from one EPSSolve() call to the other, so the
540 : deflation space should be set every time.
541 :
542 : The vectors do not need to be mutually orthonormal, since they are explicitly
543 : orthonormalized internally.
544 :
545 : Level: intermediate
546 :
547 : .seealso: EPSSetInitialSpace()
548 : @*/
549 16 : PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec v[])
550 : {
551 16 : PetscFunctionBegin;
552 16 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
553 48 : PetscValidLogicalCollectiveInt(eps,n,2);
554 16 : PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
555 16 : if (n>0) {
556 14 : PetscAssertPointer(v,3);
557 14 : PetscValidHeaderSpecific(*v,VEC_CLASSID,3);
558 : }
559 16 : PetscCall(SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl));
560 16 : if (n>0) eps->state = EPS_STATE_INITIAL;
561 16 : PetscFunctionReturn(PETSC_SUCCESS);
562 : }
563 :
564 : /*@
565 : EPSSetInitialSpace - Specify a basis of vectors that constitute the initial
566 : space, that is, the subspace from which the solver starts to iterate.
567 :
568 : Collective
569 :
570 : Input Parameters:
571 : + eps - the eigenproblem solver context
572 : . n - number of vectors
573 : - is - set of basis vectors of the initial space
574 :
575 : Notes:
576 : Some solvers start to iterate on a single vector (initial vector). In that case,
577 : the other vectors are ignored.
578 :
579 : These vectors do not persist from one EPSSolve() call to the other, so the
580 : initial space should be set every time.
581 :
582 : The vectors do not need to be mutually orthonormal, since they are explicitly
583 : orthonormalized internally.
584 :
585 : Common usage of this function is when the user can provide a rough approximation
586 : of the wanted eigenspace. Then, convergence may be faster.
587 :
588 : Level: intermediate
589 :
590 : .seealso: EPSSetLeftInitialSpace(), EPSSetDeflationSpace()
591 : @*/
592 230 : PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[])
593 : {
594 230 : PetscFunctionBegin;
595 230 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
596 690 : PetscValidLogicalCollectiveInt(eps,n,2);
597 230 : PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
598 230 : if (n>0) {
599 224 : PetscAssertPointer(is,3);
600 224 : PetscValidHeaderSpecific(*is,VEC_CLASSID,3);
601 : }
602 230 : PetscCall(SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS));
603 230 : if (n>0) eps->state = EPS_STATE_INITIAL;
604 230 : PetscFunctionReturn(PETSC_SUCCESS);
605 : }
606 :
607 : /*@
608 : EPSSetLeftInitialSpace - Specify a basis of vectors that constitute the left
609 : initial space, used by two-sided solvers to start the left subspace.
610 :
611 : Collective
612 :
613 : Input Parameters:
614 : + eps - the eigenproblem solver context
615 : . n - number of vectors
616 : - isl - set of basis vectors of the left initial space
617 :
618 : Notes:
619 : Left initial vectors are used to initiate the left search space in two-sided
620 : eigensolvers. Users should pass here an approximation of the left eigenspace,
621 : if available.
622 :
623 : The same comments in EPSSetInitialSpace() are applicable here.
624 :
625 : Level: intermediate
626 :
627 : .seealso: EPSSetInitialSpace(), EPSSetTwoSided()
628 : @*/
629 4 : PetscErrorCode EPSSetLeftInitialSpace(EPS eps,PetscInt n,Vec isl[])
630 : {
631 4 : PetscFunctionBegin;
632 4 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
633 12 : PetscValidLogicalCollectiveInt(eps,n,2);
634 4 : PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
635 4 : if (n>0) {
636 4 : PetscAssertPointer(isl,3);
637 4 : PetscValidHeaderSpecific(*isl,VEC_CLASSID,3);
638 : }
639 4 : PetscCall(SlepcBasisReference_Private(n,isl,&eps->ninil,&eps->ISL));
640 4 : if (n>0) eps->state = EPS_STATE_INITIAL;
641 4 : PetscFunctionReturn(PETSC_SUCCESS);
642 : }
643 :
644 : /*
645 : EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
646 : by the user. This is called at setup.
647 : */
648 591 : PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
649 : {
650 591 : PetscBool krylov;
651 :
652 591 : PetscFunctionBegin;
653 591 : if (*ncv!=PETSC_DETERMINE) { /* ncv set */
654 290 : PetscCall(PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,""));
655 290 : if (krylov) {
656 274 : PetscCheck(*ncv>=nev+1 || (*ncv==nev && *ncv==eps->n),PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
657 : } else {
658 16 : PetscCheck(*ncv>=nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
659 : }
660 301 : } else if (*mpd!=PETSC_DETERMINE) { /* mpd set */
661 6 : *ncv = PetscMin(eps->n,nev+(*mpd));
662 : } else { /* neither set: defaults depend on nev being small or large */
663 295 : if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
664 : else {
665 0 : *mpd = 500;
666 0 : *ncv = PetscMin(eps->n,nev+(*mpd));
667 : }
668 : }
669 591 : if (*mpd==PETSC_DETERMINE) *mpd = *ncv;
670 591 : PetscFunctionReturn(PETSC_SUCCESS);
671 : }
672 :
673 : /*@
674 : EPSAllocateSolution - Allocate memory storage for common variables such
675 : as eigenvalues and eigenvectors.
676 :
677 : Collective
678 :
679 : Input Parameters:
680 : + eps - eigensolver context
681 : - extra - number of additional positions, used for methods that require a
682 : working basis slightly larger than ncv
683 :
684 : Developer Notes:
685 : This is SLEPC_EXTERN because it may be required by user plugin EPS
686 : implementations.
687 :
688 : Level: developer
689 :
690 : .seealso: EPSSetUp()
691 : @*/
692 887 : PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
693 : {
694 887 : PetscInt oldsize,requested;
695 887 : PetscRandom rand;
696 887 : Vec t;
697 :
698 887 : PetscFunctionBegin;
699 887 : requested = eps->ncv + extra;
700 :
701 : /* oldsize is zero if this is the first time setup is called */
702 887 : PetscCall(BVGetSizes(eps->V,NULL,NULL,&oldsize));
703 :
704 : /* allocate space for eigenvalues and friends */
705 887 : if (requested != oldsize || !eps->eigr) {
706 680 : PetscCall(PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm));
707 680 : PetscCall(PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm));
708 : }
709 :
710 : /* workspace for the case of arbitrary selection */
711 887 : if (eps->arbitrary) {
712 8 : if (eps->rr) PetscCall(PetscFree2(eps->rr,eps->ri));
713 8 : PetscCall(PetscMalloc2(requested,&eps->rr,requested,&eps->ri));
714 : }
715 :
716 : /* allocate V */
717 887 : if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
718 887 : if (!oldsize) {
719 667 : if (!((PetscObject)eps->V)->type_name) PetscCall(BVSetType(eps->V,BVMAT));
720 667 : PetscCall(STMatCreateVecsEmpty(eps->st,&t,NULL));
721 667 : PetscCall(BVSetSizesFromVec(eps->V,t,requested));
722 667 : PetscCall(VecDestroy(&t));
723 220 : } else PetscCall(BVResize(eps->V,requested,PETSC_FALSE));
724 :
725 : /* allocate W */
726 887 : if (eps->twosided) {
727 19 : PetscCall(BVGetRandomContext(eps->V,&rand)); /* make sure the random context is available when duplicating */
728 19 : PetscCall(BVDestroy(&eps->W));
729 19 : PetscCall(BVDuplicate(eps->V,&eps->W));
730 : }
731 887 : PetscFunctionReturn(PETSC_SUCCESS);
732 : }
|
