Actual source code: epssetup.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    EPS routines related to problem setup
 12: */

 14: #include <slepc/private/epsimpl.h>

 16: /*
 17:    Let the solver choose the ST type that should be used by default,
 18:    otherwise set it to SHIFT.
 19:    This is called at EPSSetFromOptions (before STSetFromOptions)
 20:    and also at EPSSetUp (in case EPSSetFromOptions was not called).
 21: */
 22: PetscErrorCode EPSSetDefaultST(EPS eps)
 23: {
 24:   if (eps->ops->setdefaultst) (*eps->ops->setdefaultst)(eps);
 25:   if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSHIFT);
 26:   PetscFunctionReturn(0);
 27: }

 29: /*
 30:    This is done by preconditioned eigensolvers that use the PC only.
 31:    It sets STPRECOND with KSPPREONLY.
 32: */
 33: PetscErrorCode EPSSetDefaultST_Precond(EPS eps)
 34: {
 35:   KSP            ksp;

 37:   if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STPRECOND);
 38:   STGetKSP(eps->st,&ksp);
 39:   if (!((PetscObject)ksp)->type_name) KSPSetType(ksp,KSPPREONLY);
 40:   PetscFunctionReturn(0);
 41: }

 43: /*
 44:    This is done by preconditioned eigensolvers that can also use the KSP.
 45:    It sets STPRECOND with the default KSP (GMRES) and maxit=5.
 46: */
 47: PetscErrorCode EPSSetDefaultST_GMRES(EPS eps)
 48: {
 49:   KSP            ksp;

 51:   if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STPRECOND);
 52:   STPrecondSetKSPHasMat(eps->st,PETSC_TRUE);
 53:   STGetKSP(eps->st,&ksp);
 54:   if (!((PetscObject)ksp)->type_name) {
 55:     KSPSetType(ksp,KSPGMRES);
 56:     KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,5);
 57:   }
 58:   PetscFunctionReturn(0);
 59: }

 61: #if defined(SLEPC_HAVE_SCALAPACK) || defined(SLEPC_HAVE_ELPA) || defined(SLEPC_HAVE_ELEMENTAL) || defined(SLEPC_HAVE_EVSL)
 62: /*
 63:    This is for direct eigensolvers that work with A and B directly, so
 64:    no need to factorize B.
 65: */
 66: PetscErrorCode EPSSetDefaultST_NoFactor(EPS eps)
 67: {
 68:   KSP            ksp;
 69:   PC             pc;

 71:   if (!((PetscObject)eps->st)->type_name) STSetType(eps->st,STSHIFT);
 72:   STGetKSP(eps->st,&ksp);
 73:   if (!((PetscObject)ksp)->type_name) KSPSetType(ksp,KSPPREONLY);
 74:   KSPGetPC(ksp,&pc);
 75:   if (!((PetscObject)pc)->type_name) PCSetType(pc,PCNONE);
 76:   PetscFunctionReturn(0);
 77: }
 78: #endif

 80: /*
 81:    Check that the ST selected by the user is compatible with the EPS solver and options
 82: */
 83: PetscErrorCode EPSCheckCompatibleST(EPS eps)
 84: {
 85:   PetscBool      precond,shift,sinvert,cayley,lyapii;
 86: #if defined(PETSC_USE_COMPLEX)
 87:   PetscScalar    sigma;
 88: #endif

 90:   PetscObjectTypeCompare((PetscObject)eps->st,STPRECOND,&precond);
 91:   PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&shift);
 92:   PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert);
 93:   PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&cayley);

 95:   /* preconditioned eigensolvers */

 99:   /* harmonic extraction */

102:   /* real shifts in Hermitian problems */
103: #if defined(PETSC_USE_COMPLEX)
104:   STGetShift(eps->st,&sigma);
106: #endif

108:   /* Cayley with PGNHEP */

111:   /* make sure that the user does not specify smallest magnitude with shift-and-invert */
112:   if ((cayley || sinvert) && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER)) {
113:     PetscObjectTypeCompare((PetscObject)eps,EPSLYAPII,&lyapii);
115:   }
116:   PetscFunctionReturn(0);
117: }

119: /*
120:    MatEstimateSpectralRange_EPS: estimate the spectral range [left,right] of a
121:    symmetric/Hermitian matrix A using an auxiliary EPS object
122: */
123: PetscErrorCode MatEstimateSpectralRange_EPS(Mat A,PetscReal *left,PetscReal *right)
124: {
125:   PetscInt       nconv;
126:   PetscScalar    eig0;
127:   PetscReal      tol=1e-3,errest=tol;
128:   EPS            eps;

130:   *left = 0.0; *right = 0.0;
131:   EPSCreate(PetscObjectComm((PetscObject)A),&eps);
132:   EPSSetOptionsPrefix(eps,"eps_filter_");
133:   EPSSetOperators(eps,A,NULL);
134:   EPSSetProblemType(eps,EPS_HEP);
135:   EPSSetTolerances(eps,tol,50);
136:   EPSSetConvergenceTest(eps,EPS_CONV_ABS);
137:   EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);
138:   EPSSolve(eps);
139:   EPSGetConverged(eps,&nconv);
140:   if (nconv>0) {
141:     EPSGetEigenvalue(eps,0,&eig0,NULL);
142:     EPSGetErrorEstimate(eps,0,&errest);
143:   } else eig0 = eps->eigr[0];
144:   *left = PetscRealPart(eig0)-errest;
145:   EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
146:   EPSSolve(eps);
147:   EPSGetConverged(eps,&nconv);
148:   if (nconv>0) {
149:     EPSGetEigenvalue(eps,0,&eig0,NULL);
150:     EPSGetErrorEstimate(eps,0,&errest);
151:   } else eig0 = eps->eigr[0];
152:   *right = PetscRealPart(eig0)+errest;
153:   EPSDestroy(&eps);
154:   PetscFunctionReturn(0);
155: }

157: /*
158:    EPSSetUpSort_Basic: configure the EPS sorting criterion according to 'which'
159: */
160: PetscErrorCode EPSSetUpSort_Basic(EPS eps)
161: {
162:   switch (eps->which) {
163:     case EPS_LARGEST_MAGNITUDE:
164:       eps->sc->comparison    = SlepcCompareLargestMagnitude;
165:       eps->sc->comparisonctx = NULL;
166:       break;
167:     case EPS_SMALLEST_MAGNITUDE:
168:       eps->sc->comparison    = SlepcCompareSmallestMagnitude;
169:       eps->sc->comparisonctx = NULL;
170:       break;
171:     case EPS_LARGEST_REAL:
172:       eps->sc->comparison    = SlepcCompareLargestReal;
173:       eps->sc->comparisonctx = NULL;
174:       break;
175:     case EPS_SMALLEST_REAL:
176:       eps->sc->comparison    = SlepcCompareSmallestReal;
177:       eps->sc->comparisonctx = NULL;
178:       break;
179:     case EPS_LARGEST_IMAGINARY:
180:       eps->sc->comparison    = SlepcCompareLargestImaginary;
181:       eps->sc->comparisonctx = NULL;
182:       break;
183:     case EPS_SMALLEST_IMAGINARY:
184:       eps->sc->comparison    = SlepcCompareSmallestImaginary;
185:       eps->sc->comparisonctx = NULL;
186:       break;
187:     case EPS_TARGET_MAGNITUDE:
188:       eps->sc->comparison    = SlepcCompareTargetMagnitude;
189:       eps->sc->comparisonctx = &eps->target;
190:       break;
191:     case EPS_TARGET_REAL:
192:       eps->sc->comparison    = SlepcCompareTargetReal;
193:       eps->sc->comparisonctx = &eps->target;
194:       break;
195:     case EPS_TARGET_IMAGINARY:
196: #if defined(PETSC_USE_COMPLEX)
197:       eps->sc->comparison    = SlepcCompareTargetImaginary;
198:       eps->sc->comparisonctx = &eps->target;
199: #endif
200:       break;
201:     case EPS_ALL:
202:       eps->sc->comparison    = SlepcCompareSmallestReal;
203:       eps->sc->comparisonctx = NULL;
204:       break;
205:     case EPS_WHICH_USER:
206:       break;
207:   }
208:   eps->sc->map    = NULL;
209:   eps->sc->mapobj = NULL;
210:   PetscFunctionReturn(0);
211: }

213: /*
214:    EPSSetUpSort_Default: configure both EPS and DS sorting criterion
215: */
216: PetscErrorCode EPSSetUpSort_Default(EPS eps)
217: {
218:   SlepcSC        sc;
219:   PetscBool      istrivial;

221:   /* fill sorting criterion context */
222:   EPSSetUpSort_Basic(eps);
223:   /* fill sorting criterion for DS */
224:   DSGetSlepcSC(eps->ds,&sc);
225:   RGIsTrivial(eps->rg,&istrivial);
226:   sc->rg            = istrivial? NULL: eps->rg;
227:   sc->comparison    = eps->sc->comparison;
228:   sc->comparisonctx = eps->sc->comparisonctx;
229:   sc->map           = SlepcMap_ST;
230:   sc->mapobj        = (PetscObject)eps->st;
231:   PetscFunctionReturn(0);
232: }

234: /*@
235:    EPSSetUp - Sets up all the internal data structures necessary for the
236:    execution of the eigensolver. Then calls STSetUp() for any set-up
237:    operations associated to the ST object.

239:    Collective on eps

241:    Input Parameter:
242: .  eps   - eigenproblem solver context

244:    Notes:
245:    This function need not be called explicitly in most cases, since EPSSolve()
246:    calls it. It can be useful when one wants to measure the set-up time
247:    separately from the solve time.

249:    Level: developer

251: .seealso: EPSCreate(), EPSSolve(), EPSDestroy(), STSetUp(), EPSSetInitialSpace()
252: @*/
253: PetscErrorCode EPSSetUp(EPS eps)
254: {
255:   Mat            A,B;
256:   PetscInt       k,nmat;
257:   PetscBool      flg;

260:   if (eps->state) PetscFunctionReturn(0);
261:   PetscLogEventBegin(EPS_SetUp,eps,0,0,0);

263:   /* reset the convergence flag from the previous solves */
264:   eps->reason = EPS_CONVERGED_ITERATING;

266:   /* Set default solver type (EPSSetFromOptions was not called) */
267:   if (!((PetscObject)eps)->type_name) EPSSetType(eps,EPSKRYLOVSCHUR);
268:   if (!eps->st) EPSGetST(eps,&eps->st);
269:   EPSSetDefaultST(eps);

271:   STSetTransform(eps->st,PETSC_TRUE);
272:   if (eps->useds && !eps->ds) EPSGetDS(eps,&eps->ds);
273:   if (eps->twosided) {
275:   }
276:   if (!eps->rg) EPSGetRG(eps,&eps->rg);
277:   if (!((PetscObject)eps->rg)->type_name) RGSetType(eps->rg,RGINTERVAL);

279:   /* Set problem dimensions */
280:   STGetNumMatrices(eps->st,&nmat);
282:   STMatGetSize(eps->st,&eps->n,NULL);
283:   STMatGetLocalSize(eps->st,&eps->nloc,NULL);

285:   /* Set default problem type */
286:   if (!eps->problem_type) {
287:     if (nmat==1) EPSSetProblemType(eps,EPS_NHEP);
288:     else EPSSetProblemType(eps,EPS_GNHEP);
289:   } else if (nmat==1 && eps->isgeneralized) {
290:     PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n");
291:     eps->isgeneralized = PETSC_FALSE;
292:     eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;

295:   if (eps->nev > eps->n) eps->nev = eps->n;
296:   if (eps->ncv > eps->n) eps->ncv = eps->n;

298:   /* check some combinations of eps->which */

301:   /* initialization of matrix norms */
302:   if (eps->conv==EPS_CONV_NORM) {
303:     if (!eps->nrma) {
304:       STGetMatrix(eps->st,0,&A);
305:       MatNorm(A,NORM_INFINITY,&eps->nrma);
306:     }
307:     if (nmat>1 && !eps->nrmb) {
308:       STGetMatrix(eps->st,1,&B);
309:       MatNorm(B,NORM_INFINITY,&eps->nrmb);
310:     }
311:   }

313:   /* call specific solver setup */
314:   (*eps->ops->setup)(eps);

316:   /* if purification is set, check that it really makes sense */
317:   if (eps->purify) {
318:     if (eps->categ==EPS_CATEGORY_PRECOND || eps->categ==EPS_CATEGORY_CONTOUR) eps->purify = PETSC_FALSE;
319:     else {
320:       if (!eps->isgeneralized) eps->purify = PETSC_FALSE;
321:       else if (!eps->ishermitian && !eps->ispositive) eps->purify = PETSC_FALSE;
322:       else {
323:         PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg);
324:         if (flg) eps->purify = PETSC_FALSE;
325:       }
326:     }
327:   }

329:   /* set tolerance if not yet set */
330:   if (eps->tol==PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL;

332:   /* set up sorting criterion */
333:   if (eps->ops->setupsort) (*eps->ops->setupsort)(eps);

335:   /* Build balancing matrix if required */
336:   if (eps->balance!=EPS_BALANCE_USER) {
337:     STSetBalanceMatrix(eps->st,NULL);
338:     if (!eps->ishermitian && (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE)) {
339:       if (!eps->D) {
340:         BVCreateVec(eps->V,&eps->D);
341:         PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->D);
342:       }
343:       EPSBuildBalance_Krylov(eps);
344:       STSetBalanceMatrix(eps->st,eps->D);
345:     }
346:   }

348:   /* Setup ST */
349:   STSetUp(eps->st);
350:   EPSCheckCompatibleST(eps);

352:   /* process deflation and initial vectors */
353:   if (eps->nds<0) {
354:     k = -eps->nds;
355:     BVInsertConstraints(eps->V,&k,eps->defl);
356:     SlepcBasisDestroy_Private(&eps->nds,&eps->defl);
357:     eps->nds = k;
358:     STCheckNullSpace(eps->st,eps->V);
359:   }
360:   if (eps->nini<0) {
361:     k = -eps->nini;
363:     BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE);
364:     SlepcBasisDestroy_Private(&eps->nini,&eps->IS);
365:     eps->nini = k;
366:   }
367:   if (eps->twosided && eps->ninil<0) {
368:     k = -eps->ninil;
370:     BVInsertVecs(eps->W,0,&k,eps->ISL,PETSC_TRUE);
371:     SlepcBasisDestroy_Private(&eps->ninil,&eps->ISL);
372:     eps->ninil = k;
373:   }

375:   PetscLogEventEnd(EPS_SetUp,eps,0,0,0);
376:   eps->state = EPS_STATE_SETUP;
377:   PetscFunctionReturn(0);
378: }

380: /*@
381:    EPSSetOperators - Sets the matrices associated with the eigenvalue problem.

383:    Collective on eps

385:    Input Parameters:
386: +  eps - the eigenproblem solver context
387: .  A  - the matrix associated with the eigensystem
388: -  B  - the second matrix in the case of generalized eigenproblems

390:    Notes:
391:    To specify a standard eigenproblem, use NULL for parameter B.

393:    It must be called before EPSSetUp(). If it is called again after EPSSetUp() and
394:    the matrix sizes have changed then the EPS object is reset.

396:    Level: beginner

398: .seealso: EPSSolve(), EPSSetUp(), EPSReset(), EPSGetST(), STGetMatrix()
399: @*/
400: PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)
401: {
402:   PetscInt       m,n,m0,mloc,nloc,mloc0,nmat;
403:   Mat            mat[2];


411:   /* Check matrix sizes */
412:   MatGetSize(A,&m,&n);
413:   MatGetLocalSize(A,&mloc,&nloc);
416:   if (B) {
417:     MatGetSize(B,&m0,&n);
418:     MatGetLocalSize(B,&mloc0,&nloc);
423:   }
424:   if (eps->state && (n!=eps->n || nloc!=eps->nloc)) EPSReset(eps);
425:   eps->nrma = 0.0;
426:   eps->nrmb = 0.0;
427:   if (!eps->st) EPSGetST(eps,&eps->st);
428:   mat[0] = A;
429:   if (B) {
430:     mat[1] = B;
431:     nmat = 2;
432:   } else nmat = 1;
433:   STSetMatrices(eps->st,nmat,mat);
434:   eps->state = EPS_STATE_INITIAL;
435:   PetscFunctionReturn(0);
436: }

438: /*@
439:    EPSGetOperators - Gets the matrices associated with the eigensystem.

441:    Collective on eps

443:    Input Parameter:
444: .  eps - the EPS context

446:    Output Parameters:
447: +  A  - the matrix associated with the eigensystem
448: -  B  - the second matrix in the case of generalized eigenproblems

450:    Note:
451:    Does not increase the reference count of the matrices, so you should not destroy them.

453:    Level: intermediate

455: .seealso: EPSSolve(), EPSGetST(), STGetMatrix(), STSetMatrices()
456: @*/
457: PetscErrorCode EPSGetOperators(EPS eps,Mat *A,Mat *B)
458: {
459:   ST             st;
460:   PetscInt       k;

463:   EPSGetST(eps,&st);
464:   STGetNumMatrices(st,&k);
465:   if (A) {
466:     if (k<1) *A = NULL;
467:     else STGetMatrix(st,0,A);
468:   }
469:   if (B) {
470:     if (k<2) *B = NULL;
471:     else STGetMatrix(st,1,B);
472:   }
473:   PetscFunctionReturn(0);
474: }

476: /*@C
477:    EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
478:    space.

480:    Collective on eps

482:    Input Parameters:
483: +  eps - the eigenproblem solver context
484: .  n   - number of vectors
485: -  v   - set of basis vectors of the deflation space

487:    Notes:
488:    When a deflation space is given, the eigensolver seeks the eigensolution
489:    in the restriction of the problem to the orthogonal complement of this
490:    space. This can be used for instance in the case that an invariant
491:    subspace is known beforehand (such as the nullspace of the matrix).

493:    These vectors do not persist from one EPSSolve() call to the other, so the
494:    deflation space should be set every time.

496:    The vectors do not need to be mutually orthonormal, since they are explicitly
497:    orthonormalized internally.

499:    Level: intermediate

501: .seealso: EPSSetInitialSpace()
502: @*/
503: PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec v[])
504: {
508:   if (n>0) {
511:   }
512:   SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl);
513:   if (n>0) eps->state = EPS_STATE_INITIAL;
514:   PetscFunctionReturn(0);
515: }

517: /*@C
518:    EPSSetInitialSpace - Specify a basis of vectors that constitute the initial
519:    space, that is, the subspace from which the solver starts to iterate.

521:    Collective on eps

523:    Input Parameters:
524: +  eps - the eigenproblem solver context
525: .  n   - number of vectors
526: -  is  - set of basis vectors of the initial space

528:    Notes:
529:    Some solvers start to iterate on a single vector (initial vector). In that case,
530:    the other vectors are ignored.

532:    These vectors do not persist from one EPSSolve() call to the other, so the
533:    initial space should be set every time.

535:    The vectors do not need to be mutually orthonormal, since they are explicitly
536:    orthonormalized internally.

538:    Common usage of this function is when the user can provide a rough approximation
539:    of the wanted eigenspace. Then, convergence may be faster.

541:    Level: intermediate

543: .seealso: EPSSetLeftInitialSpace(), EPSSetDeflationSpace()
544: @*/
545: PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[])
546: {
550:   if (n>0) {
553:   }
554:   SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS);
555:   if (n>0) eps->state = EPS_STATE_INITIAL;
556:   PetscFunctionReturn(0);
557: }

559: /*@C
560:    EPSSetLeftInitialSpace - Specify a basis of vectors that constitute the left
561:    initial space, used by two-sided solvers to start the left subspace.

563:    Collective on eps

565:    Input Parameters:
566: +  eps - the eigenproblem solver context
567: .  n   - number of vectors
568: -  isl - set of basis vectors of the left initial space

570:    Notes:
571:    Left initial vectors are used to initiate the left search space in two-sided
572:    eigensolvers. Users should pass here an approximation of the left eigenspace,
573:    if available.

575:    The same comments in EPSSetInitialSpace() are applicable here.

577:    Level: intermediate

579: .seealso: EPSSetInitialSpace(), EPSSetTwoSided()
580: @*/
581: PetscErrorCode EPSSetLeftInitialSpace(EPS eps,PetscInt n,Vec isl[])
582: {
586:   if (n>0) {
589:   }
590:   SlepcBasisReference_Private(n,isl,&eps->ninil,&eps->ISL);
591:   if (n>0) eps->state = EPS_STATE_INITIAL;
592:   PetscFunctionReturn(0);
593: }

595: /*
596:   EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
597:   by the user. This is called at setup.
598:  */
599: PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
600: {
601:   PetscBool      krylov;

603:   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
604:     PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,"");
605:     if (krylov) {
607:     } else {
609:     }
610:   } else if (*mpd!=PETSC_DEFAULT) { /* mpd set */
611:     *ncv = PetscMin(eps->n,nev+(*mpd));
612:   } else { /* neither set: defaults depend on nev being small or large */
613:     if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
614:     else {
615:       *mpd = 500;
616:       *ncv = PetscMin(eps->n,nev+(*mpd));
617:     }
618:   }
619:   if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
620:   PetscFunctionReturn(0);
621: }

623: /*@
624:    EPSAllocateSolution - Allocate memory storage for common variables such
625:    as eigenvalues and eigenvectors.

627:    Collective on eps

629:    Input Parameters:
630: +  eps   - eigensolver context
631: -  extra - number of additional positions, used for methods that require a
632:            working basis slightly larger than ncv

634:    Developer Notes:
635:    This is SLEPC_EXTERN because it may be required by user plugin EPS
636:    implementations.

638:    Level: developer

640: .seealso: EPSSetUp()
641: @*/
642: PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
643: {
644:   PetscInt       oldsize,newc,requested;
645:   PetscLogDouble cnt;
646:   PetscRandom    rand;
647:   Vec            t;

649:   requested = eps->ncv + extra;

651:   /* oldsize is zero if this is the first time setup is called */
652:   BVGetSizes(eps->V,NULL,NULL,&oldsize);
653:   newc = PetscMax(0,requested-oldsize);

655:   /* allocate space for eigenvalues and friends */
656:   if (requested != oldsize || !eps->eigr) {
657:     PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm);
658:     PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm);
659:     cnt = 2*newc*sizeof(PetscScalar) + 2*newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
660:     PetscLogObjectMemory((PetscObject)eps,cnt);
661:   }

663:   /* workspace for the case of arbitrary selection */
664:   if (eps->arbitrary) {
665:     if (eps->rr) PetscFree2(eps->rr,eps->ri);
666:     PetscMalloc2(requested,&eps->rr,requested,&eps->ri);
667:     PetscLogObjectMemory((PetscObject)eps,2*newc*sizeof(PetscScalar));
668:   }

670:   /* allocate V */
671:   if (!eps->V) EPSGetBV(eps,&eps->V);
672:   if (!oldsize) {
673:     if (!((PetscObject)(eps->V))->type_name) BVSetType(eps->V,BVSVEC);
674:     STMatCreateVecsEmpty(eps->st,&t,NULL);
675:     BVSetSizesFromVec(eps->V,t,requested);
676:     VecDestroy(&t);
677:   } else BVResize(eps->V,requested,PETSC_FALSE);

679:   /* allocate W */
680:   if (eps->twosided) {
681:     BVGetRandomContext(eps->V,&rand);  /* make sure the random context is available when duplicating */
682:     BVDestroy(&eps->W);
683:     BVDuplicate(eps->V,&eps->W);
684:   }
685:   PetscFunctionReturn(0);
686: }