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 : PEP routines related to problem setup
12 : */
13 :
14 : #include <slepc/private/pepimpl.h> /*I "slepcpep.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 PEPSetFromOptions (before STSetFromOptions)
20 : and also at PEPSetUp (in case PEPSetFromOptions was not called).
21 : */
22 365 : PetscErrorCode PEPSetDefaultST(PEP pep)
23 : {
24 365 : PetscFunctionBegin;
25 365 : PetscTryTypeMethod(pep,setdefaultst);
26 365 : if (!((PetscObject)pep->st)->type_name) PetscCall(STSetType(pep->st,STSHIFT));
27 365 : PetscFunctionReturn(PETSC_SUCCESS);
28 : }
29 :
30 : /*
31 : This is used in Q-Arnoldi and STOAR to set the transform flag by
32 : default, otherwise the user has to explicitly run with -st_transform
33 : */
34 32 : PetscErrorCode PEPSetDefaultST_Transform(PEP pep)
35 : {
36 32 : PetscFunctionBegin;
37 32 : PetscCall(STSetTransform(pep->st,PETSC_TRUE));
38 32 : PetscFunctionReturn(PETSC_SUCCESS);
39 : }
40 :
41 : /*@
42 : PEPSetDSType - Sets the type of the internal DS object based on the current
43 : settings of the polynomial eigensolver.
44 :
45 : Collective
46 :
47 : Input Parameter:
48 : . pep - polynomial eigensolver context
49 :
50 : Note:
51 : This function need not be called explicitly, since it will be called at
52 : both PEPSetFromOptions() and PEPSetUp().
53 :
54 : Level: developer
55 :
56 : .seealso: PEPSetFromOptions(), PEPSetUp()
57 : @*/
58 364 : PetscErrorCode PEPSetDSType(PEP pep)
59 : {
60 364 : PetscFunctionBegin;
61 364 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
62 364 : PetscTryTypeMethod(pep,setdstype);
63 364 : PetscFunctionReturn(PETSC_SUCCESS);
64 : }
65 :
66 : /*@
67 : PEPSetUp - Sets up all the internal data structures necessary for the
68 : execution of the PEP solver.
69 :
70 : Collective
71 :
72 : Input Parameter:
73 : . pep - solver context
74 :
75 : Notes:
76 : This function need not be called explicitly in most cases, since PEPSolve()
77 : calls it. It can be useful when one wants to measure the set-up time
78 : separately from the solve time.
79 :
80 : Level: developer
81 :
82 : .seealso: PEPCreate(), PEPSolve(), PEPDestroy()
83 : @*/
84 192 : PetscErrorCode PEPSetUp(PEP pep)
85 : {
86 192 : SlepcSC sc;
87 192 : PetscBool istrivial,flg;
88 192 : PetscInt k;
89 192 : KSP ksp;
90 192 : PC pc;
91 192 : PetscMPIInt size;
92 192 : MatSolverType stype;
93 :
94 192 : PetscFunctionBegin;
95 192 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
96 192 : if (pep->state) PetscFunctionReturn(PETSC_SUCCESS);
97 189 : PetscCall(PetscLogEventBegin(PEP_SetUp,pep,0,0,0));
98 :
99 : /* reset the convergence flag from the previous solves */
100 189 : pep->reason = PEP_CONVERGED_ITERATING;
101 :
102 : /* set default solver type (PEPSetFromOptions was not called) */
103 189 : if (!((PetscObject)pep)->type_name) PetscCall(PEPSetType(pep,PEPTOAR));
104 189 : if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
105 189 : PetscCall(PEPSetDefaultST(pep));
106 189 : if (!pep->ds) PetscCall(PEPGetDS(pep,&pep->ds));
107 189 : PetscCall(PEPSetDSType(pep));
108 189 : if (!pep->rg) PetscCall(PEPGetRG(pep,&pep->rg));
109 189 : if (!((PetscObject)pep->rg)->type_name) PetscCall(RGSetType(pep->rg,RGINTERVAL));
110 :
111 : /* check matrices, transfer them to ST */
112 189 : PetscCheck(pep->A,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"PEPSetOperators must be called first");
113 189 : PetscCall(STSetMatrices(pep->st,pep->nmat,pep->A));
114 :
115 : /* set problem dimensions */
116 189 : PetscCall(MatGetSize(pep->A[0],&pep->n,NULL));
117 189 : PetscCall(MatGetLocalSize(pep->A[0],&pep->nloc,NULL));
118 :
119 : /* set default problem type */
120 189 : if (!pep->problem_type) PetscCall(PEPSetProblemType(pep,PEP_GENERAL));
121 189 : if (pep->nev > (pep->nmat-1)*pep->n) pep->nev = (pep->nmat-1)*pep->n;
122 189 : if (pep->ncv > (pep->nmat-1)*pep->n) pep->ncv = (pep->nmat-1)*pep->n;
123 :
124 : /* check consistency of refinement options */
125 189 : if (pep->refine) {
126 29 : if (!pep->scheme) { /* set default scheme */
127 4 : PetscCall(PEPRefineGetKSP(pep,&ksp));
128 4 : PetscCall(KSPGetPC(ksp,&pc));
129 4 : PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg));
130 4 : if (flg) PetscCall(PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,""));
131 8 : pep->scheme = flg? PEP_REFINE_SCHEME_MBE: PEP_REFINE_SCHEME_SCHUR;
132 : }
133 29 : if (pep->scheme==PEP_REFINE_SCHEME_MBE) {
134 8 : PetscCall(PEPRefineGetKSP(pep,&ksp));
135 8 : PetscCall(KSPGetPC(ksp,&pc));
136 8 : PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg));
137 8 : if (flg) PetscCall(PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,""));
138 8 : PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The MBE scheme for refinement requires a direct solver in KSP");
139 8 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size));
140 8 : if (size>1) { /* currently selected PC is a factorization */
141 0 : PetscCall(PCFactorGetMatSolverType(pc,&stype));
142 0 : PetscCall(PetscStrcmp(stype,MATSOLVERPETSC,&flg));
143 0 : PetscCheck(!flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"For Newton refinement, you chose to solve linear systems with a factorization, but in parallel runs you need to select an external package");
144 : }
145 : }
146 29 : if (pep->scheme==PEP_REFINE_SCHEME_SCHUR) {
147 9 : PetscCheck(pep->npart==1,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The Schur scheme for refinement does not support subcommunicators");
148 : }
149 : }
150 : /* call specific solver setup */
151 189 : PetscUseTypeMethod(pep,setup);
152 :
153 : /* set tolerance if not yet set */
154 189 : if (pep->tol==(PetscReal)PETSC_DETERMINE) pep->tol = SLEPC_DEFAULT_TOL;
155 189 : if (pep->refine) {
156 57 : if (pep->rtol==(PetscReal)PETSC_DETERMINE) pep->rtol = PetscMax(pep->tol/1000,PETSC_MACHINE_EPSILON);
157 43 : if (pep->rits==PETSC_DETERMINE) pep->rits = (pep->refine==PEP_REFINE_SIMPLE)? 10: 1;
158 : }
159 :
160 : /* set default extraction */
161 189 : if (!pep->extract) {
162 162 : pep->extract = (pep->basis==PEP_BASIS_MONOMIAL)? PEP_EXTRACT_NORM: PEP_EXTRACT_NONE;
163 : }
164 :
165 : /* fill sorting criterion context */
166 189 : switch (pep->which) {
167 52 : case PEP_LARGEST_MAGNITUDE:
168 52 : pep->sc->comparison = SlepcCompareLargestMagnitude;
169 52 : pep->sc->comparisonctx = NULL;
170 52 : break;
171 0 : case PEP_SMALLEST_MAGNITUDE:
172 0 : pep->sc->comparison = SlepcCompareSmallestMagnitude;
173 0 : pep->sc->comparisonctx = NULL;
174 0 : break;
175 2 : case PEP_LARGEST_REAL:
176 2 : pep->sc->comparison = SlepcCompareLargestReal;
177 2 : pep->sc->comparisonctx = NULL;
178 2 : break;
179 0 : case PEP_SMALLEST_REAL:
180 0 : pep->sc->comparison = SlepcCompareSmallestReal;
181 0 : pep->sc->comparisonctx = NULL;
182 0 : break;
183 0 : case PEP_LARGEST_IMAGINARY:
184 0 : pep->sc->comparison = SlepcCompareLargestImaginary;
185 0 : pep->sc->comparisonctx = NULL;
186 0 : break;
187 0 : case PEP_SMALLEST_IMAGINARY:
188 0 : pep->sc->comparison = SlepcCompareSmallestImaginary;
189 0 : pep->sc->comparisonctx = NULL;
190 0 : break;
191 110 : case PEP_TARGET_MAGNITUDE:
192 110 : pep->sc->comparison = SlepcCompareTargetMagnitude;
193 110 : pep->sc->comparisonctx = &pep->target;
194 110 : break;
195 0 : case PEP_TARGET_REAL:
196 0 : pep->sc->comparison = SlepcCompareTargetReal;
197 0 : pep->sc->comparisonctx = &pep->target;
198 0 : break;
199 0 : case PEP_TARGET_IMAGINARY:
200 : #if defined(PETSC_USE_COMPLEX)
201 0 : pep->sc->comparison = SlepcCompareTargetImaginary;
202 0 : pep->sc->comparisonctx = &pep->target;
203 : #endif
204 0 : break;
205 21 : case PEP_ALL:
206 21 : pep->sc->comparison = SlepcCompareSmallestReal;
207 21 : pep->sc->comparisonctx = NULL;
208 21 : break;
209 : case PEP_WHICH_USER:
210 : break;
211 : }
212 189 : pep->sc->map = NULL;
213 189 : pep->sc->mapobj = NULL;
214 :
215 : /* fill sorting criterion for DS */
216 189 : if (pep->which!=PEP_ALL) {
217 168 : PetscCall(DSGetSlepcSC(pep->ds,&sc));
218 168 : PetscCall(RGIsTrivial(pep->rg,&istrivial));
219 168 : sc->rg = istrivial? NULL: pep->rg;
220 168 : sc->comparison = pep->sc->comparison;
221 168 : sc->comparisonctx = pep->sc->comparisonctx;
222 168 : sc->map = SlepcMap_ST;
223 168 : sc->mapobj = (PetscObject)pep->st;
224 : }
225 : /* setup ST */
226 189 : PetscCall(STSetUp(pep->st));
227 :
228 : /* compute matrix coefficients */
229 189 : PetscCall(STGetTransform(pep->st,&flg));
230 189 : if (!flg) {
231 158 : if (pep->which!=PEP_ALL && pep->solvematcoeffs) PetscCall(STMatSetUp(pep->st,1.0,pep->solvematcoeffs));
232 : } else {
233 31 : PetscCheck(pep->basis==PEP_BASIS_MONOMIAL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Cannot use ST-transform with non-monomial basis in PEP");
234 : }
235 :
236 : /* compute scale factor if no set by user */
237 189 : PetscCall(PEPComputeScaleFactor(pep));
238 :
239 : /* build balancing matrix if required */
240 189 : if (pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) {
241 18 : if (!pep->Dl) PetscCall(BVCreateVec(pep->V,&pep->Dl));
242 18 : if (!pep->Dr) PetscCall(BVCreateVec(pep->V,&pep->Dr));
243 18 : PetscCall(PEPBuildDiagonalScaling(pep));
244 : }
245 :
246 : /* process initial vectors */
247 189 : if (pep->nini<0) {
248 3 : k = -pep->nini;
249 3 : PetscCheck(k<=pep->ncv,PetscObjectComm((PetscObject)pep),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
250 3 : PetscCall(BVInsertVecs(pep->V,0,&k,pep->IS,PETSC_TRUE));
251 3 : PetscCall(SlepcBasisDestroy_Private(&pep->nini,&pep->IS));
252 3 : pep->nini = k;
253 : }
254 189 : PetscCall(PetscLogEventEnd(PEP_SetUp,pep,0,0,0));
255 189 : pep->state = PEP_STATE_SETUP;
256 189 : PetscFunctionReturn(PETSC_SUCCESS);
257 : }
258 :
259 : /*@
260 : PEPSetOperators - Sets the coefficient matrices associated with the polynomial
261 : eigenvalue problem.
262 :
263 : Collective
264 :
265 : Input Parameters:
266 : + pep - the eigenproblem solver context
267 : . nmat - number of matrices in array A
268 : - A - the array of matrices associated with the eigenproblem
269 :
270 : Notes:
271 : The polynomial eigenproblem is defined as P(l)*x=0, where l is
272 : the eigenvalue, x is the eigenvector, and P(l) is defined as
273 : P(l) = A_0 + l*A_1 + ... + l^d*A_d, with d=nmat-1 (the degree of P).
274 : For non-monomial bases, this expression is different.
275 :
276 : Level: beginner
277 :
278 : .seealso: PEPSolve(), PEPGetOperators(), PEPGetNumMatrices(), PEPSetBasis()
279 : @*/
280 190 : PetscErrorCode PEPSetOperators(PEP pep,PetscInt nmat,Mat A[])
281 : {
282 190 : PetscInt i,n=0,m,m0=0,mloc,nloc,mloc0=0;
283 :
284 190 : PetscFunctionBegin;
285 190 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
286 570 : PetscValidLogicalCollectiveInt(pep,nmat,2);
287 190 : PetscCheck(nmat>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Non-positive value of nmat: %" PetscInt_FMT,nmat);
288 190 : PetscCheck(nmat>2,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Cannot solve linear eigenproblems with PEP; use EPS instead");
289 190 : PetscAssertPointer(A,3);
290 :
291 915 : for (i=0;i<nmat;i++) {
292 725 : PetscValidHeaderSpecific(A[i],MAT_CLASSID,3);
293 725 : PetscCheckSameComm(pep,1,A[i],3);
294 725 : PetscCall(MatGetSize(A[i],&m,&n));
295 725 : PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
296 725 : PetscCheck(m==n,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,m,n);
297 725 : PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
298 725 : if (!i) { m0 = m; mloc0 = mloc; }
299 725 : PetscCheck(m==m0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_INCOMP,"Dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,m,m0);
300 725 : PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_INCOMP,"Local dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
301 725 : PetscCall(PetscObjectReference((PetscObject)A[i]));
302 : }
303 :
304 190 : if (pep->state && (n!=pep->n || nloc!=pep->nloc)) PetscCall(PEPReset(pep));
305 186 : else if (pep->nmat) {
306 1 : PetscCall(MatDestroyMatrices(pep->nmat,&pep->A));
307 1 : PetscCall(PetscFree2(pep->pbc,pep->nrma));
308 1 : PetscCall(PetscFree(pep->solvematcoeffs));
309 : }
310 :
311 190 : PetscCall(PetscMalloc1(nmat,&pep->A));
312 190 : PetscCall(PetscCalloc2(3*nmat,&pep->pbc,nmat,&pep->nrma));
313 915 : for (i=0;i<nmat;i++) {
314 725 : pep->A[i] = A[i];
315 725 : pep->pbc[i] = 1.0; /* default to monomial basis */
316 : }
317 190 : pep->nmat = nmat;
318 190 : pep->state = PEP_STATE_INITIAL;
319 190 : PetscFunctionReturn(PETSC_SUCCESS);
320 : }
321 :
322 : /*@
323 : PEPGetOperators - Gets the matrices associated with the polynomial eigensystem.
324 :
325 : Not Collective
326 :
327 : Input Parameters:
328 : + pep - the PEP context
329 : - k - the index of the requested matrix (starting in 0)
330 :
331 : Output Parameter:
332 : . A - the requested matrix
333 :
334 : Level: intermediate
335 :
336 : .seealso: PEPSolve(), PEPSetOperators(), PEPGetNumMatrices()
337 : @*/
338 58 : PetscErrorCode PEPGetOperators(PEP pep,PetscInt k,Mat *A)
339 : {
340 58 : PetscFunctionBegin;
341 58 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
342 58 : PetscAssertPointer(A,3);
343 58 : PetscCheck(k>=0 && k<pep->nmat,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,pep->nmat-1);
344 58 : *A = pep->A[k];
345 58 : PetscFunctionReturn(PETSC_SUCCESS);
346 : }
347 :
348 : /*@
349 : PEPGetNumMatrices - Returns the number of matrices stored in the PEP.
350 :
351 : Not Collective
352 :
353 : Input Parameter:
354 : . pep - the PEP context
355 :
356 : Output Parameters:
357 : . nmat - the number of matrices passed in PEPSetOperators()
358 :
359 : Level: intermediate
360 :
361 : .seealso: PEPSetOperators()
362 : @*/
363 1 : PetscErrorCode PEPGetNumMatrices(PEP pep,PetscInt *nmat)
364 : {
365 1 : PetscFunctionBegin;
366 1 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
367 1 : PetscAssertPointer(nmat,2);
368 1 : *nmat = pep->nmat;
369 1 : PetscFunctionReturn(PETSC_SUCCESS);
370 : }
371 :
372 : /*@
373 : PEPSetInitialSpace - Specify a basis of vectors that constitute the initial
374 : space, that is, the subspace from which the solver starts to iterate.
375 :
376 : Collective
377 :
378 : Input Parameters:
379 : + pep - the polynomial eigensolver context
380 : . n - number of vectors
381 : - is - set of basis vectors of the initial space
382 :
383 : Notes:
384 : Some solvers start to iterate on a single vector (initial vector). In that case,
385 : the other vectors are ignored.
386 :
387 : These vectors do not persist from one PEPSolve() call to the other, so the
388 : initial space should be set every time.
389 :
390 : The vectors do not need to be mutually orthonormal, since they are explicitly
391 : orthonormalized internally.
392 :
393 : Common usage of this function is when the user can provide a rough approximation
394 : of the wanted eigenspace. Then, convergence may be faster.
395 :
396 : Level: intermediate
397 :
398 : .seealso: PEPSetUp()
399 : @*/
400 5 : PetscErrorCode PEPSetInitialSpace(PEP pep,PetscInt n,Vec is[])
401 : {
402 5 : PetscFunctionBegin;
403 5 : PetscValidHeaderSpecific(pep,PEP_CLASSID,1);
404 15 : PetscValidLogicalCollectiveInt(pep,n,2);
405 5 : PetscCheck(n>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
406 5 : if (n>0) {
407 5 : PetscAssertPointer(is,3);
408 5 : PetscValidHeaderSpecific(*is,VEC_CLASSID,3);
409 : }
410 5 : PetscCall(SlepcBasisReference_Private(n,is,&pep->nini,&pep->IS));
411 5 : if (n>0) pep->state = PEP_STATE_INITIAL;
412 5 : PetscFunctionReturn(PETSC_SUCCESS);
413 : }
414 :
415 : /*
416 : PEPSetDimensions_Default - Set reasonable values for ncv, mpd if not set
417 : by the user. This is called at setup.
418 : */
419 140 : PetscErrorCode PEPSetDimensions_Default(PEP pep,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
420 : {
421 140 : PetscBool krylov;
422 140 : PetscInt dim;
423 :
424 140 : PetscFunctionBegin;
425 140 : PetscCall(PetscObjectTypeCompareAny((PetscObject)pep,&krylov,PEPTOAR,PEPSTOAR,PEPQARNOLDI,""));
426 140 : dim = (pep->nmat-1)*pep->n;
427 140 : if (*ncv!=PETSC_DETERMINE) { /* ncv set */
428 70 : if (krylov) {
429 59 : PetscCheck(*ncv>nev || (*ncv==nev && *ncv==dim),PetscObjectComm((PetscObject)pep),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
430 : } else {
431 11 : PetscCheck(*ncv>=nev,PetscObjectComm((PetscObject)pep),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
432 : }
433 70 : } else if (*mpd!=PETSC_DETERMINE) { /* mpd set */
434 0 : *ncv = PetscMin(dim,nev+(*mpd));
435 : } else { /* neither set: defaults depend on nev being small or large */
436 70 : if (nev<500) *ncv = PetscMin(dim,PetscMax(2*nev,nev+15));
437 : else {
438 0 : *mpd = 500;
439 0 : *ncv = PetscMin(dim,nev+(*mpd));
440 : }
441 : }
442 140 : if (*mpd==PETSC_DETERMINE) *mpd = *ncv;
443 140 : PetscFunctionReturn(PETSC_SUCCESS);
444 : }
445 :
446 : /*@
447 : PEPAllocateSolution - Allocate memory storage for common variables such
448 : as eigenvalues and eigenvectors.
449 :
450 : Collective
451 :
452 : Input Parameters:
453 : + pep - eigensolver context
454 : - extra - number of additional positions, used for methods that require a
455 : working basis slightly larger than ncv
456 :
457 : Developer Notes:
458 : This is SLEPC_EXTERN because it may be required by user plugin PEP
459 : implementations.
460 :
461 : Level: developer
462 :
463 : .seealso: PEPSetUp()
464 : @*/
465 211 : PetscErrorCode PEPAllocateSolution(PEP pep,PetscInt extra)
466 : {
467 211 : PetscInt oldsize,requested,requestedbv;
468 211 : Vec t;
469 :
470 211 : PetscFunctionBegin;
471 211 : requested = (pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1)) + extra;
472 211 : requestedbv = pep->ncv + extra;
473 :
474 : /* oldsize is zero if this is the first time setup is called */
475 211 : PetscCall(BVGetSizes(pep->V,NULL,NULL,&oldsize));
476 :
477 : /* allocate space for eigenvalues and friends */
478 211 : if (requested != oldsize || !pep->eigr) {
479 207 : PetscCall(PetscFree4(pep->eigr,pep->eigi,pep->errest,pep->perm));
480 207 : PetscCall(PetscMalloc4(requested,&pep->eigr,requested,&pep->eigi,requested,&pep->errest,requested,&pep->perm));
481 : }
482 :
483 : /* allocate V */
484 211 : if (!pep->V) PetscCall(PEPGetBV(pep,&pep->V));
485 211 : if (!oldsize) {
486 189 : if (!((PetscObject)pep->V)->type_name) PetscCall(BVSetType(pep->V,BVMAT));
487 189 : PetscCall(STMatCreateVecsEmpty(pep->st,&t,NULL));
488 189 : PetscCall(BVSetSizesFromVec(pep->V,t,requestedbv));
489 189 : PetscCall(VecDestroy(&t));
490 22 : } else PetscCall(BVResize(pep->V,requestedbv,PETSC_FALSE));
491 211 : PetscFunctionReturn(PETSC_SUCCESS);
492 : }
|
