Actual source code: epssetup.c

slepc-3.14.0 2020-09-30
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:    This is for direct eigensolvers that work with A and B directly, so
 78:    no need to factorize B.
 79: */
 80: PetscErrorCode EPSSetDefaultST_NoFactor(EPS eps)
 81: {
 83:   KSP            ksp;
 84:   PC             pc;

 87:   if (!((PetscObject)eps->st)->type_name) {
 88:     STSetType(eps->st,STSHIFT);
 89:   }
 90:   STGetKSP(eps->st,&ksp);
 91:   if (!((PetscObject)ksp)->type_name) {
 92:     KSPSetType(ksp,KSPPREONLY);
 93:   }
 94:   KSPGetPC(ksp,&pc);
 95:   if (!((PetscObject)pc)->type_name) {
 96:     PCSetType(pc,PCNONE);
 97:   }
 98:   return(0);
 99: }

101: /*
102:    Check that the ST selected by the user is compatible with the EPS solver and options
103: */
104: PetscErrorCode EPSCheckCompatibleST(EPS eps)
105: {
107:   PetscBool      precond,shift,sinvert,cayley,lyapii;
108: #if defined(PETSC_USE_COMPLEX)
109:   PetscScalar    sigma;
110: #endif

113:   PetscObjectTypeCompare((PetscObject)eps->st,STPRECOND,&precond);
114:   PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&shift);
115:   PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert);
116:   PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&cayley);

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

122:   /* harmonic extraction */
123:   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");

125:   /* real shifts in Hermitian problems */
126: #if defined(PETSC_USE_COMPLEX)
127:   STGetShift(eps->st,&sigma);
128:   if (eps->ishermitian && PetscImaginaryPart(sigma) != 0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Hermitian problems are not compatible with complex shifts");
129: #endif

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

134:   /* make sure that the user does not specify smallest magnitude with shift-and-invert */
135:   if ((cayley || sinvert) && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER)) {
136:     PetscObjectTypeCompare((PetscObject)eps,EPSLYAPII,&lyapii);
137:     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");
138:   }
139:   return(0);
140: }

142: /*
143:    EPSSetUpSort_Basic: configure the EPS sorting criterion according to 'which'
144: */
145: PetscErrorCode EPSSetUpSort_Basic(EPS eps)
146: {
148:   switch (eps->which) {
149:     case EPS_LARGEST_MAGNITUDE:
150:       eps->sc->comparison    = SlepcCompareLargestMagnitude;
151:       eps->sc->comparisonctx = NULL;
152:       break;
153:     case EPS_SMALLEST_MAGNITUDE:
154:       eps->sc->comparison    = SlepcCompareSmallestMagnitude;
155:       eps->sc->comparisonctx = NULL;
156:       break;
157:     case EPS_LARGEST_REAL:
158:       eps->sc->comparison    = SlepcCompareLargestReal;
159:       eps->sc->comparisonctx = NULL;
160:       break;
161:     case EPS_SMALLEST_REAL:
162:       eps->sc->comparison    = SlepcCompareSmallestReal;
163:       eps->sc->comparisonctx = NULL;
164:       break;
165:     case EPS_LARGEST_IMAGINARY:
166:       eps->sc->comparison    = SlepcCompareLargestImaginary;
167:       eps->sc->comparisonctx = NULL;
168:       break;
169:     case EPS_SMALLEST_IMAGINARY:
170:       eps->sc->comparison    = SlepcCompareSmallestImaginary;
171:       eps->sc->comparisonctx = NULL;
172:       break;
173:     case EPS_TARGET_MAGNITUDE:
174:       eps->sc->comparison    = SlepcCompareTargetMagnitude;
175:       eps->sc->comparisonctx = &eps->target;
176:       break;
177:     case EPS_TARGET_REAL:
178:       eps->sc->comparison    = SlepcCompareTargetReal;
179:       eps->sc->comparisonctx = &eps->target;
180:       break;
181:     case EPS_TARGET_IMAGINARY:
182: #if defined(PETSC_USE_COMPLEX)
183:       eps->sc->comparison    = SlepcCompareTargetImaginary;
184:       eps->sc->comparisonctx = &eps->target;
185: #endif
186:       break;
187:     case EPS_ALL:
188:       eps->sc->comparison    = SlepcCompareSmallestReal;
189:       eps->sc->comparisonctx = NULL;
190:       break;
191:     case EPS_WHICH_USER:
192:       break;
193:   }
194:   eps->sc->map    = NULL;
195:   eps->sc->mapobj = NULL;
196:   return(0);
197: }

199: /*
200:    EPSSetUpSort_Default: configure both EPS and DS sorting criterion
201: */
202: PetscErrorCode EPSSetUpSort_Default(EPS eps)
203: {
205:   SlepcSC        sc;
206:   PetscBool      istrivial;

209:   /* fill sorting criterion context */
210:   EPSSetUpSort_Basic(eps);
211:   /* fill sorting criterion for DS */
212:   DSGetSlepcSC(eps->ds,&sc);
213:   RGIsTrivial(eps->rg,&istrivial);
214:   sc->rg            = istrivial? NULL: eps->rg;
215:   sc->comparison    = eps->sc->comparison;
216:   sc->comparisonctx = eps->sc->comparisonctx;
217:   sc->map           = SlepcMap_ST;
218:   sc->mapobj        = (PetscObject)eps->st;
219:   if (eps->twosided) {
220:     DSSetSlepcSC(eps->dsts,sc);
221:   }
222:   return(0);
223: }

225: /*@
226:    EPSSetUp - Sets up all the internal data structures necessary for the
227:    execution of the eigensolver. Then calls STSetUp() for any set-up
228:    operations associated to the ST object.

230:    Collective on eps

232:    Input Parameter:
233: .  eps   - eigenproblem solver context

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

240:    Level: developer

242: .seealso: EPSCreate(), EPSSolve(), EPSDestroy(), STSetUp(), EPSSetInitialSpace()
243: @*/
244: PetscErrorCode EPSSetUp(EPS eps)
245: {
247:   Mat            A,B;
248:   PetscInt       k,nmat;
249:   PetscBool      flg;

253:   if (eps->state) return(0);
254:   PetscLogEventBegin(EPS_SetUp,eps,0,0,0);

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

259:   /* Set default solver type (EPSSetFromOptions was not called) */
260:   if (!((PetscObject)eps)->type_name) {
261:     EPSSetType(eps,EPSKRYLOVSCHUR);
262:   }
263:   if (!eps->st) { EPSGetST(eps,&eps->st); }
264:   EPSSetDefaultST(eps);

266:   STSetTransform(eps->st,PETSC_TRUE);
267:   if (eps->useds && !eps->ds) { EPSGetDS(eps,&eps->ds); }
268:   if (eps->twosided) {
269:     if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive)) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided methods are not intended for %s problems",SLEPC_STRING_HERMITIAN);
270:     if (!eps->dsts) { DSDuplicate(eps->ds,&eps->dsts); }
271:   }
272:   if (!eps->rg) { EPSGetRG(eps,&eps->rg); }
273:   if (!((PetscObject)eps->rg)->type_name) {
274:     RGSetType(eps->rg,RGINTERVAL);
275:   }

277:   /* Set problem dimensions */
278:   STGetNumMatrices(eps->st,&nmat);
279:   if (!nmat) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSSetOperators must be called first");
280:   STMatGetSize(eps->st,&eps->n,NULL);
281:   STMatGetLocalSize(eps->st,&eps->nloc,NULL);

283:   /* Set default problem type */
284:   if (!eps->problem_type) {
285:     if (nmat==1) {
286:       EPSSetProblemType(eps,EPS_NHEP);
287:     } else {
288:       EPSSetProblemType(eps,EPS_GNHEP);
289:     }
290:   } else if (nmat==1 && eps->isgeneralized) {
291:     PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n");
292:     eps->isgeneralized = PETSC_FALSE;
293:     eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;
294:   } else if (nmat>1 && !eps->isgeneralized) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Inconsistent EPS state: the problem type does not match the number of matrices");

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

299:   /* check some combinations of eps->which */
300:   if (eps->ishermitian && (!eps->isgeneralized || eps->ispositive) && (eps->which==EPS_LARGEST_IMAGINARY || eps->which==EPS_SMALLEST_IMAGINARY || eps->which==EPS_TARGET_IMAGINARY)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Sorting the eigenvalues along the imaginary axis is not allowed when all eigenvalues are real");

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

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

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

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

333:   /* set up sorting criterion */
334:   if (eps->ops->setupsort) {
335:     (*eps->ops->setupsort)(eps);
336:   }

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

351:   /* Setup ST */
352:   STSetUp(eps->st);
353:   EPSCheckCompatibleST(eps);

355:   /* process deflation and initial vectors */
356:   if (eps->nds<0) {
357:     k = -eps->nds;
358:     BVInsertConstraints(eps->V,&k,eps->defl);
359:     SlepcBasisDestroy_Private(&eps->nds,&eps->defl);
360:     eps->nds = k;
361:     STCheckNullSpace(eps->st,eps->V);
362:   }
363:   if (eps->nini<0) {
364:     k = -eps->nini;
365:     if (k>eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The number of initial vectors is larger than ncv");
366:     BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE);
367:     SlepcBasisDestroy_Private(&eps->nini,&eps->IS);
368:     eps->nini = k;
369:   }
370:   if (eps->twosided && eps->ninil<0) {
371:     k = -eps->ninil;
372:     if (k>eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The number of left initial vectors is larger than ncv");
373:     BVInsertVecs(eps->W,0,&k,eps->ISL,PETSC_TRUE);
374:     SlepcBasisDestroy_Private(&eps->ninil,&eps->ISL);
375:     eps->ninil = k;
376:   }

378:   PetscLogEventEnd(EPS_SetUp,eps,0,0,0);
379:   eps->state = EPS_STATE_SETUP;
380:   return(0);
381: }

383: /*@
384:    EPSSetOperators - Sets the matrices associated with the eigenvalue problem.

386:    Collective on eps

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

393:    Notes:
394:    To specify a standard eigenproblem, use NULL for parameter B.

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

399:    Level: beginner

401: .seealso: EPSSolve(), EPSSetUp(), EPSReset(), EPSGetST(), STGetMatrix()
402: @*/
403: PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)
404: {
406:   PetscInt       m,n,m0,nmat;
407:   Mat            mat[2];


416:   /* Check for square matrices */
417:   MatGetSize(A,&m,&n);
418:   if (m!=n) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A is a non-square matrix");
419:   if (B) {
420:     MatGetSize(B,&m0,&n);
421:     if (m0!=n) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B is a non-square matrix");
422:     if (m!=m0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Dimensions of A and B do not match");
423:   }
424:   if (eps->state && n!=eps->n) { 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:   return(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:    Level: intermediate

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

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

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

478:    Collective on eps

480:    Input Parameter:
481: +  eps - the eigenproblem solver context
482: .  n   - number of vectors
483: -  v   - set of basis vectors of the deflation space

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

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

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

497:    Level: intermediate

499: .seealso: EPSSetInitialSpace()
500: @*/
501: PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec v[])
502: {

508:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n out of range");
509:   if (n>0) {
512:   }
513:   SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl);
514:   if (n>0) eps->state = EPS_STATE_INITIAL;
515:   return(0);
516: }

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

522:    Collective on eps

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

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

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

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

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

542:    Level: intermediate

544: .seealso: EPSSetLeftInitialSpace(), EPSSetDeflationSpace()
545: @*/
546: PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[])
547: {

553:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
554:   if (n>0) {
557:   }
558:   SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS);
559:   if (n>0) eps->state = EPS_STATE_INITIAL;
560:   return(0);
561: }

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

567:    Collective on eps

569:    Input Parameter:
570: +  eps - the eigenproblem solver context
571: .  n   - number of vectors
572: -  isl - set of basis vectors of the left initial space

574:    Notes:
575:    Left initial vectors are used to initiate the left search space in two-sided
576:    eigensolvers. Users should pass here an approximation of the left eigenspace,
577:    if available.

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

581:    Level: intermediate

583: .seealso: EPSSetInitialSpace(), EPSSetTwoSided()
584: @*/
585: PetscErrorCode EPSSetLeftInitialSpace(EPS eps,PetscInt n,Vec isl[])
586: {

592:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
593:   if (n>0) {
596:   }
597:   SlepcBasisReference_Private(n,isl,&eps->ninil,&eps->ISL);
598:   if (n>0) eps->state = EPS_STATE_INITIAL;
599:   return(0);
600: }

602: /*
603:   EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
604:   by the user. This is called at setup.
605:  */
606: PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
607: {
609:   PetscBool      krylov;

612:   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
613:     PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,"");
614:     if (krylov) {
615:       if (*ncv<nev+1 && !(*ncv==nev && *ncv==eps->n)) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev+1");
616:     } else {
617:       if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev");
618:     }
619:   } else if (*mpd!=PETSC_DEFAULT) { /* mpd set */
620:     *ncv = PetscMin(eps->n,nev+(*mpd));
621:   } else { /* neither set: defaults depend on nev being small or large */
622:     if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
623:     else {
624:       *mpd = 500;
625:       *ncv = PetscMin(eps->n,nev+(*mpd));
626:     }
627:   }
628:   if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
629:   return(0);
630: }

632: /*@
633:    EPSAllocateSolution - Allocate memory storage for common variables such
634:    as eigenvalues and eigenvectors.

636:    Collective on eps

638:    Input Parameters:
639: +  eps   - eigensolver context
640: -  extra - number of additional positions, used for methods that require a
641:            working basis slightly larger than ncv

643:    Developers Note:
644:    This is SLEPC_EXTERN because it may be required by user plugin EPS
645:    implementations.

647:    Level: developer
648: @*/
649: PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
650: {
652:   PetscInt       oldsize,newc,requested;
653:   PetscLogDouble cnt;
654:   Vec            t;

657:   requested = eps->ncv + extra;

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

663:   /* allocate space for eigenvalues and friends */
664:   if (requested != oldsize || !eps->eigr) {
665:     PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm);
666:     PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm);
667:     cnt = 2*newc*sizeof(PetscScalar) + 2*newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
668:     PetscLogObjectMemory((PetscObject)eps,cnt);
669:   }

671:   /* workspace for the case of arbitrary selection */
672:   if (eps->arbitrary) {
673:     if (eps->rr) {
674:       PetscFree2(eps->rr,eps->ri);
675:     }
676:     PetscMalloc2(requested,&eps->rr,requested,&eps->ri);
677:     PetscLogObjectMemory((PetscObject)eps,2*newc*sizeof(PetscScalar));
678:   }

680:   /* allocate V */
681:   if (!eps->V) { EPSGetBV(eps,&eps->V); }
682:   if (!oldsize) {
683:     if (!((PetscObject)(eps->V))->type_name) {
684:       BVSetType(eps->V,BVSVEC);
685:     }
686:     STMatCreateVecsEmpty(eps->st,&t,NULL);
687:     BVSetSizesFromVec(eps->V,t,requested);
688:     VecDestroy(&t);
689:   } else {
690:     BVResize(eps->V,requested,PETSC_FALSE);
691:   }

693:   /* allocate W */
694:   if (eps->twosided) {
695:     BVDestroy(&eps->W);
696:     BVDuplicate(eps->V,&eps->W);
697:   }
698:   return(0);
699: }