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