Actual source code: epssetup.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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: #if defined(SLEPC_HAVE_SCALAPACK) || defined(SLEPC_HAVE_ELPA) || defined(SLEPC_HAVE_ELEMENTAL) || defined(SLEPC_HAVE_EVSL)
 77: /*
 78:    This is for direct eigensolvers that work with A and B directly, so
 79:    no need to factorize B.
 80: */
 81: PetscErrorCode EPSSetDefaultST_NoFactor(EPS eps)
 82: {
 84:   KSP            ksp;
 85:   PC             pc;

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

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

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

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

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

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

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

136:   /* make sure that the user does not specify smallest magnitude with shift-and-invert */
137:   if ((cayley || sinvert) && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER)) {
138:     PetscObjectTypeCompare((PetscObject)eps,EPSLYAPII,&lyapii);
139:     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),PETSC_ERR_USER_INPUT,"Shift-and-invert requires a target 'which' (see EPSSetWhichEigenpairs), for instance -st_type sinvert -eps_target 0 -eps_target_magnitude");
140:   }
141:   return(0);
142: }

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

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

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

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

229:    Collective on eps

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

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

239:    Level: developer

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

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

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

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

265:   STSetTransform(eps->st,PETSC_TRUE);
266:   if (eps->useds && !eps->ds) { EPSGetDS(eps,&eps->ds); }
267:   if (eps->twosided) {
268:     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);
269:   }
270:   if (!eps->rg) { EPSGetRG(eps,&eps->rg); }
271:   if (!((PetscObject)eps->rg)->type_name) {
272:     RGSetType(eps->rg,RGINTERVAL);
273:   }

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

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

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

297:   /* check some combinations of eps->which */
298:   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");

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

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

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

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

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

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

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

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

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

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

384:    Collective on eps

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

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

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

397:    Level: beginner

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


414:   /* Check matrix sizes */
415:   MatGetSize(A,&m,&n);
416:   MatGetLocalSize(A,&mloc,&nloc);
417:   if (m!=n) SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A is a non-square matrix (%D rows, %D cols)",m,n);
418:   if (mloc!=nloc) SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A does not have equal row and column sizes (%D, %D)",mloc,nloc);
419:   if (B) {
420:     MatGetSize(B,&m0,&n);
421:     MatGetLocalSize(B,&mloc0,&nloc);
422:     if (m0!=n) SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B is a non-square matrix (%D rows, %D cols)",m0,n);
423:     if (mloc0!=nloc) SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B does not have equal row and column local sizes (%D, %D)",mloc0,nloc);
424:     if (m!=m0) SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Dimensions of A and B do not match (%D, %D)",m,m0);
425:     if (mloc!=mloc0) SETERRQ2(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Local dimensions of A and B do not match (%D, %D)",mloc,mloc0);
426:   }
427:   if (eps->state && (n!=eps->n || nloc!=eps->nloc)) { EPSReset(eps); }
428:   eps->nrma = 0.0;
429:   eps->nrmb = 0.0;
430:   if (!eps->st) { EPSGetST(eps,&eps->st); }
431:   mat[0] = A;
432:   if (B) {
433:     mat[1] = B;
434:     nmat = 2;
435:   } else nmat = 1;
436:   STSetMatrices(eps->st,nmat,mat);
437:   eps->state = EPS_STATE_INITIAL;
438:   return(0);
439: }

441: /*@
442:    EPSGetOperators - Gets the matrices associated with the eigensystem.

444:    Collective on eps

446:    Input Parameter:
447: .  eps - the EPS context

449:    Output Parameters:
450: +  A  - the matrix associated with the eigensystem
451: -  B  - the second matrix in the case of generalized eigenproblems

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

456:    Level: intermediate

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

468:   EPSGetST(eps,&st);
469:   STGetNumMatrices(st,&k);
470:   if (A) {
471:     if (k<1) *A = NULL;
472:     else {
473:       STGetMatrix(st,0,A);
474:     }
475:   }
476:   if (B) {
477:     if (k<2) *B = NULL;
478:     else {
479:       STGetMatrix(st,1,B);
480:     }
481:   }
482:   return(0);
483: }

485: /*@C
486:    EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
487:    space.

489:    Collective on eps

491:    Input Parameters:
492: +  eps - the eigenproblem solver context
493: .  n   - number of vectors
494: -  v   - set of basis vectors of the deflation space

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

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

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

508:    Level: intermediate

510: .seealso: EPSSetInitialSpace()
511: @*/
512: PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec v[])
513: {

519:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n out of range");
520:   if (n>0) {
523:   }
524:   SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl);
525:   if (n>0) eps->state = EPS_STATE_INITIAL;
526:   return(0);
527: }

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

533:    Collective on eps

535:    Input Parameters:
536: +  eps - the eigenproblem solver context
537: .  n   - number of vectors
538: -  is  - set of basis vectors of the initial space

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

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

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

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

553:    Level: intermediate

555: .seealso: EPSSetLeftInitialSpace(), EPSSetDeflationSpace()
556: @*/
557: PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[])
558: {

564:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
565:   if (n>0) {
568:   }
569:   SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS);
570:   if (n>0) eps->state = EPS_STATE_INITIAL;
571:   return(0);
572: }

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

578:    Collective on eps

580:    Input Parameters:
581: +  eps - the eigenproblem solver context
582: .  n   - number of vectors
583: -  isl - set of basis vectors of the left initial space

585:    Notes:
586:    Left initial vectors are used to initiate the left search space in two-sided
587:    eigensolvers. Users should pass here an approximation of the left eigenspace,
588:    if available.

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

592:    Level: intermediate

594: .seealso: EPSSetInitialSpace(), EPSSetTwoSided()
595: @*/
596: PetscErrorCode EPSSetLeftInitialSpace(EPS eps,PetscInt n,Vec isl[])
597: {

603:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
604:   if (n>0) {
607:   }
608:   SlepcBasisReference_Private(n,isl,&eps->ninil,&eps->ISL);
609:   if (n>0) eps->state = EPS_STATE_INITIAL;
610:   return(0);
611: }

613: /*
614:   EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
615:   by the user. This is called at setup.
616:  */
617: PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
618: {
620:   PetscBool      krylov;

623:   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
624:     PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,"");
625:     if (krylov) {
626:       if (*ncv<nev+1 && !(*ncv==nev && *ncv==eps->n)) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
627:     } else {
628:       if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
629:     }
630:   } else if (*mpd!=PETSC_DEFAULT) { /* mpd set */
631:     *ncv = PetscMin(eps->n,nev+(*mpd));
632:   } else { /* neither set: defaults depend on nev being small or large */
633:     if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
634:     else {
635:       *mpd = 500;
636:       *ncv = PetscMin(eps->n,nev+(*mpd));
637:     }
638:   }
639:   if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
640:   return(0);
641: }

643: /*@
644:    EPSAllocateSolution - Allocate memory storage for common variables such
645:    as eigenvalues and eigenvectors.

647:    Collective on eps

649:    Input Parameters:
650: +  eps   - eigensolver context
651: -  extra - number of additional positions, used for methods that require a
652:            working basis slightly larger than ncv

654:    Developers Note:
655:    This is SLEPC_EXTERN because it may be required by user plugin EPS
656:    implementations.

658:    Level: developer
659: @*/
660: PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
661: {
663:   PetscInt       oldsize,newc,requested;
664:   PetscLogDouble cnt;
665:   PetscRandom    rand;
666:   Vec            t;

669:   requested = eps->ncv + extra;

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

675:   /* allocate space for eigenvalues and friends */
676:   if (requested != oldsize || !eps->eigr) {
677:     PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm);
678:     PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm);
679:     cnt = 2*newc*sizeof(PetscScalar) + 2*newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
680:     PetscLogObjectMemory((PetscObject)eps,cnt);
681:   }

683:   /* workspace for the case of arbitrary selection */
684:   if (eps->arbitrary) {
685:     if (eps->rr) {
686:       PetscFree2(eps->rr,eps->ri);
687:     }
688:     PetscMalloc2(requested,&eps->rr,requested,&eps->ri);
689:     PetscLogObjectMemory((PetscObject)eps,2*newc*sizeof(PetscScalar));
690:   }

692:   /* allocate V */
693:   if (!eps->V) { EPSGetBV(eps,&eps->V); }
694:   if (!oldsize) {
695:     if (!((PetscObject)(eps->V))->type_name) {
696:       BVSetType(eps->V,BVSVEC);
697:     }
698:     STMatCreateVecsEmpty(eps->st,&t,NULL);
699:     BVSetSizesFromVec(eps->V,t,requested);
700:     VecDestroy(&t);
701:   } else {
702:     BVResize(eps->V,requested,PETSC_FALSE);
703:   }

705:   /* allocate W */
706:   if (eps->twosided) {
707:     BVGetRandomContext(eps->V,&rand);  /* make sure the random context is available when duplicating */
708:     BVDestroy(&eps->W);
709:     BVDuplicate(eps->V,&eps->W);
710:   }
711:   return(0);
712: }