Actual source code: epssetup.c

slepc-3.13.3 2020-06-14
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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: {

 27:   if (eps->ops->setdefaultst) { (*eps->ops->setdefaultst)(eps); }
 28:   if (!((PetscObject)eps->st)->type_name) {
 29:     STSetType(eps->st,STSHIFT);
 30:   }
 31:   return(0);
 32: }

 34: /*
 35:    This is done by preconditioned eigensolvers that use the PC only.
 36:    It sets STPRECOND with KSPPREONLY.
 37: */
 38: PetscErrorCode EPSSetDefaultST_Precond(EPS eps)
 39: {
 41:   KSP            ksp;

 44:   if (!((PetscObject)eps->st)->type_name) {
 45:     STSetType(eps->st,STPRECOND);
 46:   }
 47:   STGetKSP(eps->st,&ksp);
 48:   if (!((PetscObject)ksp)->type_name) {
 49:     KSPSetType(ksp,KSPPREONLY);
 50:   }
 51:   return(0);
 52: }

 54: /*
 55:    This is done by preconditioned eigensolvers that can also use the KSP.
 56:    It sets STPRECOND with the default KSP (GMRES) and maxit=5.
 57: */
 58: PetscErrorCode EPSSetDefaultST_GMRES(EPS eps)
 59: {
 61:   KSP            ksp;

 64:   if (!((PetscObject)eps->st)->type_name) {
 65:     STSetType(eps->st,STPRECOND);
 66:   }
 67:   STPrecondSetKSPHasMat(eps->st,PETSC_TRUE);
 68:   STGetKSP(eps->st,&ksp);
 69:   if (!((PetscObject)ksp)->type_name) {
 70:     KSPSetType(ksp,KSPGMRES);
 71:     KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,5);
 72:   }
 73:   return(0);
 74: }

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

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

 93:   /* preconditioned eigensolvers */
 94:   if (eps->categ==EPS_CATEGORY_PRECOND && !precond) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires ST=PRECOND");
 95:   if (eps->categ!=EPS_CATEGORY_PRECOND && precond) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"STPRECOND is intended for preconditioned eigensolvers only");

 97:   /* harmonic extraction */
 98:   if (!(precond || shift) && eps->extraction && eps->extraction!=EPS_RITZ) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot use a spectral transformation combined with harmonic extraction");

100:   /* real shifts in Hermitian problems */
101: #if defined(PETSC_USE_COMPLEX)
102:   STGetShift(eps->st,&sigma);
103:   if (eps->ishermitian && PetscImaginaryPart(sigma) != 0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Hermitian problems are not compatible with complex shifts");
104: #endif

106:   /* Cayley with PGNHEP */
107:   if (cayley && eps->problem_type == EPS_PGNHEP) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cayley spectral transformation is not compatible with PGNHEP");

109:   /* make sure that the user does not specify smallest magnitude with shift-and-invert */
110:   if ((cayley || sinvert) && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER)) {
111:     PetscObjectTypeCompare((PetscObject)eps,EPSLYAPII,&lyapii);
112:     if (!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) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Shift-and-invert requires a target 'which' (see EPSSetWhichEigenpairs), for instance -st_type sinvert -eps_target 0 -eps_target_magnitude");
113:   }
114:   return(0);
115: }

117: /*@
118:    EPSSetUp - Sets up all the internal data structures necessary for the
119:    execution of the eigensolver. Then calls STSetUp() for any set-up
120:    operations associated to the ST object.

122:    Collective on eps

124:    Input Parameter:
125: .  eps   - eigenproblem solver context

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

132:    Level: developer

134: .seealso: EPSCreate(), EPSSolve(), EPSDestroy(), STSetUp(), EPSSetInitialSpace()
135: @*/
136: PetscErrorCode EPSSetUp(EPS eps)
137: {
139:   Mat            A,B;
140:   SlepcSC        sc;
141:   PetscInt       k,nmat;
142:   PetscBool      flg,istrivial;

146:   if (eps->state) return(0);
147:   PetscLogEventBegin(EPS_SetUp,eps,0,0,0);

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

152:   /* Set default solver type (EPSSetFromOptions was not called) */
153:   if (!((PetscObject)eps)->type_name) {
154:     EPSSetType(eps,EPSKRYLOVSCHUR);
155:   }
156:   if (!eps->st) { EPSGetST(eps,&eps->st); }
157:   EPSSetDefaultST(eps);

159:   STSetTransform(eps->st,PETSC_TRUE);
160:   if (eps->useds && !eps->ds) { EPSGetDS(eps,&eps->ds); }
161:   if (eps->twosided) {
162:     if (!eps->hasts) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing left eigenvectors (no two-sided variant)");
163:     if (eps->ishermitian) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided methods are not intended for symmetric problems");
164:     if (!eps->dsts) { DSDuplicate(eps->ds,&eps->dsts); }
165:   }
166:   if (!eps->rg) { EPSGetRG(eps,&eps->rg); }
167:   if (!((PetscObject)eps->rg)->type_name) {
168:     RGSetType(eps->rg,RGINTERVAL);
169:   }

171:   /* Set problem dimensions */
172:   STGetNumMatrices(eps->st,&nmat);
173:   if (!nmat) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSSetOperators must be called first");
174:   STMatGetSize(eps->st,&eps->n,NULL);
175:   STMatGetLocalSize(eps->st,&eps->nloc,NULL);

177:   /* Set default problem type */
178:   if (!eps->problem_type) {
179:     if (nmat==1) {
180:       EPSSetProblemType(eps,EPS_NHEP);
181:     } else {
182:       EPSSetProblemType(eps,EPS_GNHEP);
183:     }
184:   } else if (nmat==1 && eps->isgeneralized) {
185:     PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n");
186:     eps->isgeneralized = PETSC_FALSE;
187:     eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;
188:   } else if (nmat>1 && !eps->isgeneralized) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Inconsistent EPS state");

190:   if (eps->nev > eps->n) eps->nev = eps->n;
191:   if (eps->ncv > eps->n) eps->ncv = eps->n;

193:   /* initialization of matrix norms */
194:   if (eps->conv==EPS_CONV_NORM) {
195:     if (!eps->nrma) {
196:       STGetMatrix(eps->st,0,&A);
197:       MatNorm(A,NORM_INFINITY,&eps->nrma);
198:     }
199:     if (nmat>1 && !eps->nrmb) {
200:       STGetMatrix(eps->st,1,&B);
201:       MatNorm(B,NORM_INFINITY,&eps->nrmb);
202:     }
203:   }

205:   /* call specific solver setup */
206:   (*eps->ops->setup)(eps);

208:   /* if purification is set, check that it really makes sense */
209:   if (eps->purify) {
210:     if (eps->categ==EPS_CATEGORY_PRECOND || eps->categ==EPS_CATEGORY_CONTOUR) eps->purify = PETSC_FALSE;
211:     else {
212:       if (!eps->isgeneralized) eps->purify = PETSC_FALSE;
213:       else if (!eps->ishermitian && !eps->ispositive) eps->purify = PETSC_FALSE;
214:       else {
215:         PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg);
216:         if (flg) eps->purify = PETSC_FALSE;
217:       }
218:     }
219:   }

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

224:   /* fill sorting criterion context */
225:   switch (eps->which) {
226:     case EPS_LARGEST_MAGNITUDE:
227:       eps->sc->comparison    = SlepcCompareLargestMagnitude;
228:       eps->sc->comparisonctx = NULL;
229:       break;
230:     case EPS_SMALLEST_MAGNITUDE:
231:       eps->sc->comparison    = SlepcCompareSmallestMagnitude;
232:       eps->sc->comparisonctx = NULL;
233:       break;
234:     case EPS_LARGEST_REAL:
235:       eps->sc->comparison    = SlepcCompareLargestReal;
236:       eps->sc->comparisonctx = NULL;
237:       break;
238:     case EPS_SMALLEST_REAL:
239:       eps->sc->comparison    = SlepcCompareSmallestReal;
240:       eps->sc->comparisonctx = NULL;
241:       break;
242:     case EPS_LARGEST_IMAGINARY:
243:       eps->sc->comparison    = SlepcCompareLargestImaginary;
244:       eps->sc->comparisonctx = NULL;
245:       break;
246:     case EPS_SMALLEST_IMAGINARY:
247:       eps->sc->comparison    = SlepcCompareSmallestImaginary;
248:       eps->sc->comparisonctx = NULL;
249:       break;
250:     case EPS_TARGET_MAGNITUDE:
251:       eps->sc->comparison    = SlepcCompareTargetMagnitude;
252:       eps->sc->comparisonctx = &eps->target;
253:       break;
254:     case EPS_TARGET_REAL:
255:       eps->sc->comparison    = SlepcCompareTargetReal;
256:       eps->sc->comparisonctx = &eps->target;
257:       break;
258:     case EPS_TARGET_IMAGINARY:
259: #if defined(PETSC_USE_COMPLEX)
260:       eps->sc->comparison    = SlepcCompareTargetImaginary;
261:       eps->sc->comparisonctx = &eps->target;
262: #endif
263:       break;
264:     case EPS_ALL:
265:       eps->sc->comparison    = SlepcCompareSmallestReal;
266:       eps->sc->comparisonctx = NULL;
267:       break;
268:     case EPS_WHICH_USER:
269:       break;
270:   }
271:   eps->sc->map    = NULL;
272:   eps->sc->mapobj = NULL;

274:   /* fill sorting criterion for DS */
275:   if (eps->useds && eps->which!=EPS_ALL) {
276:     DSGetSlepcSC(eps->ds,&sc);
277:     RGIsTrivial(eps->rg,&istrivial);
278:     sc->rg            = istrivial? NULL: eps->rg;
279:     sc->comparison    = eps->sc->comparison;
280:     sc->comparisonctx = eps->sc->comparisonctx;
281:     sc->map           = SlepcMap_ST;
282:     sc->mapobj        = (PetscObject)eps->st;
283:     if (eps->twosided) {
284:       DSSetSlepcSC(eps->dsts,sc);
285:     }
286:   }

288:   /* Build balancing matrix if required */
289:   if (eps->balance!=EPS_BALANCE_USER) {
290:     STSetBalanceMatrix(eps->st,NULL);
291:     if (!eps->ishermitian && (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE)) {
292:       if (!eps->D) {
293:         BVCreateVec(eps->V,&eps->D);
294:         PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->D);
295:       }
296:       EPSBuildBalance_Krylov(eps);
297:       STSetBalanceMatrix(eps->st,eps->D);
298:     }
299:   }

301:   /* Setup ST */
302:   STSetUp(eps->st);
303:   EPSCheckCompatibleST(eps);

305:   /* process deflation and initial vectors */
306:   if (eps->nds<0) {
307:     k = -eps->nds;
308:     BVInsertConstraints(eps->V,&k,eps->defl);
309:     SlepcBasisDestroy_Private(&eps->nds,&eps->defl);
310:     eps->nds = k;
311:     STCheckNullSpace(eps->st,eps->V);
312:   }
313:   if (eps->nini<0) {
314:     k = -eps->nini;
315:     if (k>eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The number of initial vectors is larger than ncv");
316:     BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE);
317:     SlepcBasisDestroy_Private(&eps->nini,&eps->IS);
318:     eps->nini = k;
319:   }
320:   if (eps->twosided && eps->ninil<0) {
321:     k = -eps->ninil;
322:     if (k>eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The number of left initial vectors is larger than ncv");
323:     BVInsertVecs(eps->W,0,&k,eps->ISL,PETSC_TRUE);
324:     SlepcBasisDestroy_Private(&eps->ninil,&eps->ISL);
325:     eps->ninil = k;
326:   }

328:   PetscLogEventEnd(EPS_SetUp,eps,0,0,0);
329:   eps->state = EPS_STATE_SETUP;
330:   return(0);
331: }

333: /*@
334:    EPSSetOperators - Sets the matrices associated with the eigenvalue problem.

336:    Collective on eps

338:    Input Parameters:
339: +  eps - the eigenproblem solver context
340: .  A  - the matrix associated with the eigensystem
341: -  B  - the second matrix in the case of generalized eigenproblems

343:    Notes:
344:    To specify a standard eigenproblem, use NULL for parameter B.

346:    It must be called before EPSSetUp(). If it is called again after EPSSetUp() then
347:    the EPS object is reset.

349:    Level: beginner

351: .seealso: EPSSolve(), EPSSetUp(), EPSReset(), EPSGetST(), STGetMatrix()
352: @*/
353: PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)
354: {
356:   PetscInt       m,n,m0,nmat;
357:   Mat            mat[2];


366:   /* Check for square matrices */
367:   MatGetSize(A,&m,&n);
368:   if (m!=n) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A is a non-square matrix");
369:   if (B) {
370:     MatGetSize(B,&m0,&n);
371:     if (m0!=n) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B is a non-square matrix");
372:     if (m!=m0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Dimensions of A and B do not match");
373:   }
374:   if (eps->state && n!=eps->n) { EPSReset(eps); }
375:   eps->nrma = 0.0;
376:   eps->nrmb = 0.0;
377:   if (!eps->st) { EPSGetST(eps,&eps->st); }
378:   mat[0] = A;
379:   if (B) {
380:     mat[1] = B;
381:     nmat = 2;
382:   } else nmat = 1;
383:   STSetMatrices(eps->st,nmat,mat);
384:   eps->state = EPS_STATE_INITIAL;
385:   return(0);
386: }

388: /*@
389:    EPSGetOperators - Gets the matrices associated with the eigensystem.

391:    Collective on eps

393:    Input Parameter:
394: .  eps - the EPS context

396:    Output Parameters:
397: +  A  - the matrix associated with the eigensystem
398: -  B  - the second matrix in the case of generalized eigenproblems

400:    Level: intermediate

402: .seealso: EPSSolve(), EPSGetST(), STGetMatrix(), STSetMatrices()
403: @*/
404: PetscErrorCode EPSGetOperators(EPS eps,Mat *A,Mat *B)
405: {
407:   ST             st;
408:   PetscInt       k;

412:   EPSGetST(eps,&st);
413:   if (A) { STGetMatrix(st,0,A); }
414:   if (B) {
415:     STGetNumMatrices(st,&k);
416:     if (k==1) B = NULL;
417:     else {
418:       STGetMatrix(st,1,B);
419:     }
420:   }
421:   return(0);
422: }

424: /*@C
425:    EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
426:    space.

428:    Collective on eps

430:    Input Parameter:
431: +  eps - the eigenproblem solver context
432: .  n   - number of vectors
433: -  v   - set of basis vectors of the deflation space

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

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

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

447:    Level: intermediate

449: .seealso: EPSSetInitialSpace()
450: @*/
451: PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec v[])
452: {

458:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n out of range");
459:   if (n>0) {
462:   }
463:   SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl);
464:   if (n>0) eps->state = EPS_STATE_INITIAL;
465:   return(0);
466: }

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

472:    Collective on eps

474:    Input Parameter:
475: +  eps - the eigenproblem solver context
476: .  n   - number of vectors
477: -  is  - set of basis vectors of the initial space

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

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

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

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

492:    Level: intermediate

494: .seealso: EPSSetLeftInitialSpace(), EPSSetDeflationSpace()
495: @*/
496: PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[])
497: {

503:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
504:   if (n>0) {
507:   }
508:   SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS);
509:   if (n>0) eps->state = EPS_STATE_INITIAL;
510:   return(0);
511: }

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

517:    Collective on eps

519:    Input Parameter:
520: +  eps - the eigenproblem solver context
521: .  n   - number of vectors
522: -  isl - set of basis vectors of the left initial space

524:    Notes:
525:    Left initial vectors are used to initiate the left search space in two-sided
526:    eigensolvers. Users should pass here an approximation of the left eigenspace,
527:    if available.

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

531:    Level: intermediate

533: .seealso: EPSSetInitialSpace(), EPSSetTwoSided()
534: @*/
535: PetscErrorCode EPSSetLeftInitialSpace(EPS eps,PetscInt n,Vec isl[])
536: {

542:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
543:   if (n>0) {
546:   }
547:   SlepcBasisReference_Private(n,isl,&eps->ninil,&eps->ISL);
548:   if (n>0) eps->state = EPS_STATE_INITIAL;
549:   return(0);
550: }

552: /*
553:   EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
554:   by the user. This is called at setup.
555:  */
556: PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
557: {
559:   PetscBool      krylov;

562:   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
563:     PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,"");
564:     if (krylov) {
565:       if (*ncv<nev+1 && !(*ncv==nev && *ncv==eps->n)) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev+1");
566:     } else {
567:       if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev");
568:     }
569:   } else if (*mpd!=PETSC_DEFAULT) { /* mpd set */
570:     *ncv = PetscMin(eps->n,nev+(*mpd));
571:   } else { /* neither set: defaults depend on nev being small or large */
572:     if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
573:     else {
574:       *mpd = 500;
575:       *ncv = PetscMin(eps->n,nev+(*mpd));
576:     }
577:   }
578:   if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
579:   return(0);
580: }

582: /*@
583:    EPSAllocateSolution - Allocate memory storage for common variables such
584:    as eigenvalues and eigenvectors.

586:    Collective on eps

588:    Input Parameters:
589: +  eps   - eigensolver context
590: -  extra - number of additional positions, used for methods that require a
591:            working basis slightly larger than ncv

593:    Developers Note:
594:    This is SLEPC_EXTERN because it may be required by user plugin EPS
595:    implementations.

597:    Level: developer
598: @*/
599: PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
600: {
602:   PetscInt       oldsize,newc,requested;
603:   PetscLogDouble cnt;
604:   Vec            t;

607:   requested = eps->ncv + extra;

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

613:   /* allocate space for eigenvalues and friends */
614:   if (requested != oldsize || !eps->eigr) {
615:     PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm);
616:     PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm);
617:     cnt = 2*newc*sizeof(PetscScalar) + 2*newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
618:     PetscLogObjectMemory((PetscObject)eps,cnt);
619:   }

621:   /* workspace for the case of arbitrary selection */
622:   if (eps->arbitrary) {
623:     if (eps->rr) {
624:       PetscFree2(eps->rr,eps->ri);
625:     }
626:     PetscMalloc2(requested,&eps->rr,requested,&eps->ri);
627:     PetscLogObjectMemory((PetscObject)eps,2*newc*sizeof(PetscScalar));
628:   }

630:   /* allocate V */
631:   if (!eps->V) { EPSGetBV(eps,&eps->V); }
632:   if (!oldsize) {
633:     if (!((PetscObject)(eps->V))->type_name) {
634:       BVSetType(eps->V,BVSVEC);
635:     }
636:     STMatCreateVecsEmpty(eps->st,&t,NULL);
637:     BVSetSizesFromVec(eps->V,t,requested);
638:     VecDestroy(&t);
639:   } else {
640:     BVResize(eps->V,requested,PETSC_FALSE);
641:   }

643:   /* allocate W */
644:   if (eps->twosided) {
645:     BVDestroy(&eps->W);
646:     BVDuplicate(eps->V,&eps->W);
647:   }
648:   return(0);
649: }