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 1467 : PetscErrorCode EPSSetDefaultST(EPS eps)
23 : {
24 1467 : PetscFunctionBegin;
25 1467 : PetscTryTypeMethod(eps,setdefaultst);
26 1467 : if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSHIFT));
27 1467 : 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_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,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 873 : static PetscErrorCode EPSCheckCompatibleST(EPS eps)
88 : {
89 873 : PetscBool precond,shift,sinvert,cayley,lyapii;
90 : #if defined(PETSC_USE_COMPLEX)
91 : PetscScalar sigma;
92 : #endif
93 :
94 873 : PetscFunctionBegin;
95 873 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STPRECOND,&precond));
96 873 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&shift));
97 873 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert));
98 873 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&cayley));
99 :
100 : /* preconditioned eigensolvers */
101 873 : PetscCheck(eps->categ!=EPS_CATEGORY_PRECOND || precond,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires ST=PRECOND");
102 873 : 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 873 : 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 873 : 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 873 : if ((cayley || sinvert) && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER)) {
118 119 : PetscCall(PetscObjectTypeCompare((PetscObject)eps,EPSLYAPII,&lyapii));
119 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 873 : 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 812 : PetscErrorCode EPSSetUpSort_Basic(EPS eps)
167 : {
168 812 : PetscFunctionBegin;
169 812 : 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 3 : case EPS_SMALLEST_MAGNITUDE:
175 3 : eps->sc->comparison = SlepcCompareSmallestMagnitude;
176 3 : eps->sc->comparisonctx = NULL;
177 3 : 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 114 : case EPS_TARGET_MAGNITUDE:
195 114 : eps->sc->comparison = SlepcCompareTargetMagnitude;
196 114 : eps->sc->comparisonctx = &eps->target;
197 114 : 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 812 : eps->sc->map = NULL;
216 812 : eps->sc->mapobj = NULL;
217 812 : PetscFunctionReturn(PETSC_SUCCESS);
218 : }
219 :
220 : /*
221 : EPSSetUpSort_Default: configure both EPS and DS sorting criterion
222 : */
223 765 : PetscErrorCode EPSSetUpSort_Default(EPS eps)
224 : {
225 765 : SlepcSC sc;
226 765 : PetscBool istrivial;
227 :
228 765 : PetscFunctionBegin;
229 : /* fill sorting criterion context */
230 765 : PetscCall(EPSSetUpSort_Basic(eps));
231 : /* fill sorting criterion for DS */
232 765 : PetscCall(DSGetSlepcSC(eps->ds,&sc));
233 765 : PetscCall(RGIsTrivial(eps->rg,&istrivial));
234 765 : sc->rg = istrivial? NULL: eps->rg;
235 765 : sc->comparison = eps->sc->comparison;
236 765 : sc->comparisonctx = eps->sc->comparisonctx;
237 765 : sc->map = SlepcMap_ST;
238 765 : sc->mapobj = (PetscObject)eps->st;
239 765 : 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 1300 : PetscErrorCode EPSSetDSType(EPS eps)
260 : {
261 1300 : PetscFunctionBegin;
262 1300 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
263 1300 : PetscTryTypeMethod(eps,setdstype);
264 1300 : 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 1048 : PetscErrorCode EPSSetUp(EPS eps)
287 : {
288 1048 : Mat A,B;
289 1048 : PetscInt k,nmat;
290 1048 : PetscBool flg;
291 :
292 1048 : PetscFunctionBegin;
293 1048 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
294 1048 : if (eps->state) PetscFunctionReturn(PETSC_SUCCESS);
295 873 : PetscCall(PetscLogEventBegin(EPS_SetUp,eps,0,0,0));
296 :
297 : /* reset the convergence flag from the previous solves */
298 873 : eps->reason = EPS_CONVERGED_ITERATING;
299 :
300 : /* Set default solver type (EPSSetFromOptions was not called) */
301 873 : if (!((PetscObject)eps)->type_name) PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
302 873 : if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
303 873 : PetscCall(EPSSetDefaultST(eps));
304 :
305 873 : PetscCall(STSetTransform(eps->st,PETSC_TRUE));
306 873 : if (eps->useds && !eps->ds) PetscCall(EPSGetDS(eps,&eps->ds));
307 873 : if (eps->useds) PetscCall(EPSSetDSType(eps));
308 873 : if (eps->twosided) {
309 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);
310 : }
311 873 : if (!eps->rg) PetscCall(EPSGetRG(eps,&eps->rg));
312 873 : if (!((PetscObject)eps->rg)->type_name) PetscCall(RGSetType(eps->rg,RGINTERVAL));
313 :
314 : /* Set problem dimensions */
315 873 : PetscCall(STGetNumMatrices(eps->st,&nmat));
316 873 : PetscCheck(nmat,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSSetOperators must be called first");
317 873 : PetscCall(STMatGetSize(eps->st,&eps->n,NULL));
318 873 : PetscCall(STMatGetLocalSize(eps->st,&eps->nloc,NULL));
319 :
320 : /* Set default problem type */
321 873 : if (!eps->problem_type) {
322 41 : if (nmat==1) PetscCall(EPSSetProblemType(eps,EPS_NHEP));
323 5 : else PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
324 832 : } else if (nmat==1 && eps->isgeneralized) {
325 0 : PetscCall(PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n"));
326 0 : eps->isgeneralized = PETSC_FALSE;
327 0 : eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;
328 832 : } 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");
329 :
330 873 : if (eps->nev > eps->n) eps->nev = eps->n;
331 873 : if (eps->ncv > eps->n) eps->ncv = eps->n;
332 :
333 : /* check some combinations of eps->which */
334 873 : 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");
335 :
336 : /* initialization of matrix norms */
337 873 : if (eps->conv==EPS_CONV_NORM) {
338 43 : if (!eps->nrma) {
339 43 : PetscCall(STGetMatrix(eps->st,0,&A));
340 43 : PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
341 : }
342 43 : if (nmat>1 && !eps->nrmb) {
343 42 : PetscCall(STGetMatrix(eps->st,1,&B));
344 42 : PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
345 : }
346 : }
347 :
348 : /* call specific solver setup */
349 873 : PetscUseTypeMethod(eps,setup);
350 :
351 : /* if purification is set, check that it really makes sense */
352 873 : if (eps->purify) {
353 627 : if (eps->categ==EPS_CATEGORY_PRECOND || eps->categ==EPS_CATEGORY_CONTOUR) eps->purify = PETSC_FALSE;
354 : else {
355 456 : if (!eps->isgeneralized) eps->purify = PETSC_FALSE;
356 147 : else if (!eps->ishermitian && !eps->ispositive) eps->purify = PETSC_FALSE;
357 : else {
358 120 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg));
359 120 : if (flg) eps->purify = PETSC_FALSE;
360 : }
361 : }
362 : }
363 :
364 : /* set tolerance if not yet set */
365 873 : if (eps->tol==(PetscReal)PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL;
366 :
367 : /* set up sorting criterion */
368 873 : PetscTryTypeMethod(eps,setupsort);
369 :
370 : /* Build balancing matrix if required */
371 873 : if (eps->balance!=EPS_BALANCE_USER) {
372 873 : PetscCall(STSetBalanceMatrix(eps->st,NULL));
373 873 : if (!eps->ishermitian && (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE)) {
374 13 : if (!eps->D) PetscCall(BVCreateVec(eps->V,&eps->D));
375 13 : PetscCall(EPSBuildBalance_Krylov(eps));
376 13 : PetscCall(STSetBalanceMatrix(eps->st,eps->D));
377 : }
378 : }
379 :
380 : /* Setup ST */
381 873 : PetscCall(STSetUp(eps->st));
382 873 : PetscCall(EPSCheckCompatibleST(eps));
383 :
384 : /* process deflation and initial vectors */
385 873 : if (eps->nds<0) {
386 14 : k = -eps->nds;
387 14 : PetscCall(BVInsertConstraints(eps->V,&k,eps->defl));
388 14 : PetscCall(SlepcBasisDestroy_Private(&eps->nds,&eps->defl));
389 14 : eps->nds = k;
390 14 : PetscCall(STCheckNullSpace(eps->st,eps->V));
391 : }
392 873 : if (eps->nini<0) {
393 224 : k = -eps->nini;
394 224 : PetscCheck(k<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
395 224 : PetscCall(BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE));
396 224 : PetscCall(SlepcBasisDestroy_Private(&eps->nini,&eps->IS));
397 224 : eps->nini = k;
398 : }
399 873 : if (eps->twosided && eps->ninil<0) {
400 4 : k = -eps->ninil;
401 4 : PetscCheck(k<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
402 4 : PetscCall(BVInsertVecs(eps->W,0,&k,eps->ISL,PETSC_TRUE));
403 4 : PetscCall(SlepcBasisDestroy_Private(&eps->ninil,&eps->ISL));
404 4 : eps->ninil = k;
405 : }
406 :
407 873 : PetscCall(PetscLogEventEnd(EPS_SetUp,eps,0,0,0));
408 873 : eps->state = EPS_STATE_SETUP;
409 873 : PetscFunctionReturn(PETSC_SUCCESS);
410 : }
411 :
412 : /*@
413 : EPSSetOperators - Sets the matrices associated with the eigenvalue problem.
414 :
415 : Collective
416 :
417 : Input Parameters:
418 : + eps - the eigenproblem solver context
419 : . A - the matrix associated with the eigensystem
420 : - B - the second matrix in the case of generalized eigenproblems
421 :
422 : Notes:
423 : To specify a standard eigenproblem, use NULL for parameter B.
424 :
425 : It must be called before EPSSetUp(). If it is called again after EPSSetUp() and
426 : the matrix sizes have changed then the EPS object is reset.
427 :
428 : Level: beginner
429 :
430 : .seealso: EPSSolve(), EPSSetUp(), EPSReset(), EPSGetST(), STGetMatrix()
431 : @*/
432 736 : PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)
433 : {
434 736 : PetscInt m,n,m0,mloc,nloc,mloc0,nmat;
435 736 : Mat mat[2];
436 :
437 736 : PetscFunctionBegin;
438 736 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
439 736 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
440 736 : if (B) PetscValidHeaderSpecific(B,MAT_CLASSID,3);
441 736 : PetscCheckSameComm(eps,1,A,2);
442 736 : if (B) PetscCheckSameComm(eps,1,B,3);
443 :
444 : /* Check matrix sizes */
445 736 : PetscCall(MatGetSize(A,&m,&n));
446 736 : PetscCall(MatGetLocalSize(A,&mloc,&nloc));
447 736 : PetscCheck(m==n,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
448 736 : 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);
449 736 : if (B) {
450 231 : PetscCall(MatGetSize(B,&m0,&n));
451 231 : PetscCall(MatGetLocalSize(B,&mloc0,&nloc));
452 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);
453 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);
454 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);
455 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);
456 : }
457 736 : if (eps->state && (n!=eps->n || nloc!=eps->nloc)) PetscCall(EPSReset(eps));
458 736 : eps->nrma = 0.0;
459 736 : eps->nrmb = 0.0;
460 736 : if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
461 736 : mat[0] = A;
462 736 : if (B) {
463 231 : mat[1] = B;
464 231 : nmat = 2;
465 : } else nmat = 1;
466 736 : PetscCall(STSetMatrices(eps->st,nmat,mat));
467 736 : eps->state = EPS_STATE_INITIAL;
468 736 : PetscFunctionReturn(PETSC_SUCCESS);
469 : }
470 :
471 : /*@
472 : EPSGetOperators - Gets the matrices associated with the eigensystem.
473 :
474 : Collective
475 :
476 : Input Parameter:
477 : . eps - the EPS context
478 :
479 : Output Parameters:
480 : + A - the matrix associated with the eigensystem
481 : - B - the second matrix in the case of generalized eigenproblems
482 :
483 : Note:
484 : Does not increase the reference count of the matrices, so you should not destroy them.
485 :
486 : Level: intermediate
487 :
488 : .seealso: EPSSolve(), EPSGetST(), STGetMatrix(), STSetMatrices()
489 : @*/
490 878 : PetscErrorCode EPSGetOperators(EPS eps,Mat *A,Mat *B)
491 : {
492 878 : ST st;
493 878 : PetscInt k;
494 :
495 878 : PetscFunctionBegin;
496 878 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
497 878 : PetscCall(EPSGetST(eps,&st));
498 878 : PetscCall(STGetNumMatrices(st,&k));
499 878 : if (A) {
500 441 : if (k<1) *A = NULL;
501 441 : else PetscCall(STGetMatrix(st,0,A));
502 : }
503 878 : if (B) {
504 594 : if (k<2) *B = NULL;
505 565 : else PetscCall(STGetMatrix(st,1,B));
506 : }
507 878 : PetscFunctionReturn(PETSC_SUCCESS);
508 : }
509 :
510 : /*@C
511 : EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
512 : space.
513 :
514 : Collective
515 :
516 : Input Parameters:
517 : + eps - the eigenproblem solver context
518 : . n - number of vectors
519 : - v - set of basis vectors of the deflation space
520 :
521 : Notes:
522 : When a deflation space is given, the eigensolver seeks the eigensolution
523 : in the restriction of the problem to the orthogonal complement of this
524 : space. This can be used for instance in the case that an invariant
525 : subspace is known beforehand (such as the nullspace of the matrix).
526 :
527 : These vectors do not persist from one EPSSolve() call to the other, so the
528 : deflation space should be set every time.
529 :
530 : The vectors do not need to be mutually orthonormal, since they are explicitly
531 : orthonormalized internally.
532 :
533 : Level: intermediate
534 :
535 : .seealso: EPSSetInitialSpace()
536 : @*/
537 16 : PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec v[])
538 : {
539 16 : PetscFunctionBegin;
540 16 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
541 64 : PetscValidLogicalCollectiveInt(eps,n,2);
542 16 : PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
543 16 : if (n>0) {
544 14 : PetscAssertPointer(v,3);
545 14 : PetscValidHeaderSpecific(*v,VEC_CLASSID,3);
546 : }
547 16 : PetscCall(SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl));
548 16 : if (n>0) eps->state = EPS_STATE_INITIAL;
549 16 : PetscFunctionReturn(PETSC_SUCCESS);
550 : }
551 :
552 : /*@C
553 : EPSSetInitialSpace - Specify a basis of vectors that constitute the initial
554 : space, that is, the subspace from which the solver starts to iterate.
555 :
556 : Collective
557 :
558 : Input Parameters:
559 : + eps - the eigenproblem solver context
560 : . n - number of vectors
561 : - is - set of basis vectors of the initial space
562 :
563 : Notes:
564 : Some solvers start to iterate on a single vector (initial vector). In that case,
565 : the other vectors are ignored.
566 :
567 : These vectors do not persist from one EPSSolve() call to the other, so the
568 : initial space should be set every time.
569 :
570 : The vectors do not need to be mutually orthonormal, since they are explicitly
571 : orthonormalized internally.
572 :
573 : Common usage of this function is when the user can provide a rough approximation
574 : of the wanted eigenspace. Then, convergence may be faster.
575 :
576 : Level: intermediate
577 :
578 : .seealso: EPSSetLeftInitialSpace(), EPSSetDeflationSpace()
579 : @*/
580 230 : PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[])
581 : {
582 230 : PetscFunctionBegin;
583 230 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
584 920 : PetscValidLogicalCollectiveInt(eps,n,2);
585 230 : PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
586 230 : if (n>0) {
587 224 : PetscAssertPointer(is,3);
588 224 : PetscValidHeaderSpecific(*is,VEC_CLASSID,3);
589 : }
590 230 : PetscCall(SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS));
591 230 : if (n>0) eps->state = EPS_STATE_INITIAL;
592 230 : PetscFunctionReturn(PETSC_SUCCESS);
593 : }
594 :
595 : /*@C
596 : EPSSetLeftInitialSpace - Specify a basis of vectors that constitute the left
597 : initial space, used by two-sided solvers to start the left subspace.
598 :
599 : Collective
600 :
601 : Input Parameters:
602 : + eps - the eigenproblem solver context
603 : . n - number of vectors
604 : - isl - set of basis vectors of the left initial space
605 :
606 : Notes:
607 : Left initial vectors are used to initiate the left search space in two-sided
608 : eigensolvers. Users should pass here an approximation of the left eigenspace,
609 : if available.
610 :
611 : The same comments in EPSSetInitialSpace() are applicable here.
612 :
613 : Level: intermediate
614 :
615 : .seealso: EPSSetInitialSpace(), EPSSetTwoSided()
616 : @*/
617 4 : PetscErrorCode EPSSetLeftInitialSpace(EPS eps,PetscInt n,Vec isl[])
618 : {
619 4 : PetscFunctionBegin;
620 4 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
621 16 : PetscValidLogicalCollectiveInt(eps,n,2);
622 4 : PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
623 4 : if (n>0) {
624 4 : PetscAssertPointer(isl,3);
625 4 : PetscValidHeaderSpecific(*isl,VEC_CLASSID,3);
626 : }
627 4 : PetscCall(SlepcBasisReference_Private(n,isl,&eps->ninil,&eps->ISL));
628 4 : if (n>0) eps->state = EPS_STATE_INITIAL;
629 4 : PetscFunctionReturn(PETSC_SUCCESS);
630 : }
631 :
632 : /*
633 : EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
634 : by the user. This is called at setup.
635 : */
636 577 : PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
637 : {
638 577 : PetscBool krylov;
639 :
640 577 : PetscFunctionBegin;
641 577 : if (*ncv!=PETSC_DEFAULT) { /* ncv set */
642 276 : PetscCall(PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,""));
643 276 : if (krylov) {
644 260 : 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");
645 : } else {
646 16 : PetscCheck(*ncv>=nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
647 : }
648 301 : } else if (*mpd!=PETSC_DEFAULT) { /* mpd set */
649 6 : *ncv = PetscMin(eps->n,nev+(*mpd));
650 : } else { /* neither set: defaults depend on nev being small or large */
651 295 : if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
652 : else {
653 0 : *mpd = 500;
654 0 : *ncv = PetscMin(eps->n,nev+(*mpd));
655 : }
656 : }
657 577 : if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
658 577 : PetscFunctionReturn(PETSC_SUCCESS);
659 : }
660 :
661 : /*@
662 : EPSAllocateSolution - Allocate memory storage for common variables such
663 : as eigenvalues and eigenvectors.
664 :
665 : Collective
666 :
667 : Input Parameters:
668 : + eps - eigensolver context
669 : - extra - number of additional positions, used for methods that require a
670 : working basis slightly larger than ncv
671 :
672 : Developer Notes:
673 : This is SLEPC_EXTERN because it may be required by user plugin EPS
674 : implementations.
675 :
676 : Level: developer
677 :
678 : .seealso: EPSSetUp()
679 : @*/
680 873 : PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
681 : {
682 873 : PetscInt oldsize,requested;
683 873 : PetscRandom rand;
684 873 : Vec t;
685 :
686 873 : PetscFunctionBegin;
687 873 : requested = eps->ncv + extra;
688 :
689 : /* oldsize is zero if this is the first time setup is called */
690 873 : PetscCall(BVGetSizes(eps->V,NULL,NULL,&oldsize));
691 :
692 : /* allocate space for eigenvalues and friends */
693 873 : if (requested != oldsize || !eps->eigr) {
694 666 : PetscCall(PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm));
695 666 : PetscCall(PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm));
696 : }
697 :
698 : /* workspace for the case of arbitrary selection */
699 873 : if (eps->arbitrary) {
700 8 : if (eps->rr) PetscCall(PetscFree2(eps->rr,eps->ri));
701 8 : PetscCall(PetscMalloc2(requested,&eps->rr,requested,&eps->ri));
702 : }
703 :
704 : /* allocate V */
705 873 : if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
706 873 : if (!oldsize) {
707 653 : if (!((PetscObject)eps->V)->type_name) PetscCall(BVSetType(eps->V,BVMAT));
708 653 : PetscCall(STMatCreateVecsEmpty(eps->st,&t,NULL));
709 653 : PetscCall(BVSetSizesFromVec(eps->V,t,requested));
710 653 : PetscCall(VecDestroy(&t));
711 220 : } else PetscCall(BVResize(eps->V,requested,PETSC_FALSE));
712 :
713 : /* allocate W */
714 873 : if (eps->twosided) {
715 19 : PetscCall(BVGetRandomContext(eps->V,&rand)); /* make sure the random context is available when duplicating */
716 19 : PetscCall(BVDestroy(&eps->W));
717 19 : PetscCall(BVDuplicate(eps->V,&eps->W));
718 : }
719 873 : PetscFunctionReturn(PETSC_SUCCESS);
720 : }
|
