Actual source code: epssetup.c

slepc-3.22.1 2024-10-28
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:   PetscFunctionBegin;
 25:   PetscTryTypeMethod(eps,setdefaultst);
 26:   if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSHIFT));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

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

 38:   PetscFunctionBegin;
 39:   if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STPRECOND));
 40:   PetscCall(STGetKSP(eps->st,&ksp));
 41:   if (!((PetscObject)ksp)->type_name) PetscCall(KSPSetType(ksp,KSPPREONLY));
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

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

 53:   PetscFunctionBegin;
 54:   if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STPRECOND));
 55:   PetscCall(STPrecondSetKSPHasMat(eps->st,PETSC_TRUE));
 56:   PetscCall(STGetKSP(eps->st,&ksp));
 57:   if (!((PetscObject)ksp)->type_name) {
 58:     PetscCall(KSPSetType(ksp,KSPGMRES));
 59:     PetscCall(KSPSetTolerances(ksp,PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT,5));
 60:   }
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

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

 74:   PetscFunctionBegin;
 75:   if (!((PetscObject)eps->st)->type_name) PetscCall(STSetType(eps->st,STSHIFT));
 76:   PetscCall(STGetKSP(eps->st,&ksp));
 77:   if (!((PetscObject)ksp)->type_name) PetscCall(KSPSetType(ksp,KSPPREONLY));
 78:   PetscCall(KSPGetPC(ksp,&pc));
 79:   if (!((PetscObject)pc)->type_name) PetscCall(PCSetType(pc,PCNONE));
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }
 82: #endif

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

 94:   PetscFunctionBegin;
 95:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STPRECOND,&precond));
 96:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSHIFT,&shift));
 97:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STSINVERT,&sinvert));
 98:   PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&cayley));

100:   /* preconditioned eigensolvers */
101:   PetscCheck(eps->categ!=EPS_CATEGORY_PRECOND || precond,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires ST=PRECOND");
102:   PetscCheck(eps->categ==EPS_CATEGORY_PRECOND || !precond,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"STPRECOND is intended for preconditioned eigensolvers only");

104:   /* harmonic extraction */
105:   PetscCheck(precond || shift || !eps->extraction || eps->extraction==EPS_RITZ,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot use a spectral transformation combined with harmonic extraction");

107:   /* real shifts in Hermitian problems */
108: #if defined(PETSC_USE_COMPLEX)
109:   PetscCall(STGetShift(eps->st,&sigma));
110:   PetscCheck(!eps->ishermitian || PetscImaginaryPart(sigma)==0.0,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Hermitian problems are not compatible with complex shifts");
111: #endif

113:   /* Cayley with PGNHEP */
114:   PetscCheck(!cayley || eps->problem_type!=EPS_PGNHEP,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cayley spectral transformation is not compatible with PGNHEP");

116:   /* make sure that the user does not specify smallest magnitude with shift-and-invert */
117:   if ((cayley || sinvert) && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER)) {
118:     PetscCall(PetscObjectTypeCompare((PetscObject)eps,EPSLYAPII,&lyapii));
119:     PetscCheck(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,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");
120:   }
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }

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

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

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

220: /*
221:    EPSSetUpSort_Default: configure both EPS and DS sorting criterion
222: */
223: PetscErrorCode EPSSetUpSort_Default(EPS eps)
224: {
225:   SlepcSC        sc;
226:   PetscBool      istrivial;

228:   PetscFunctionBegin;
229:   /* fill sorting criterion context */
230:   PetscCall(EPSSetUpSort_Basic(eps));
231:   /* fill sorting criterion for DS */
232:   PetscCall(DSGetSlepcSC(eps->ds,&sc));
233:   PetscCall(RGIsTrivial(eps->rg,&istrivial));
234:   sc->rg            = istrivial? NULL: eps->rg;
235:   sc->comparison    = eps->sc->comparison;
236:   sc->comparisonctx = eps->sc->comparisonctx;
237:   sc->map           = SlepcMap_ST;
238:   sc->mapobj        = (PetscObject)eps->st;
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

242: /*@
243:    EPSSetDSType - Sets the type of the internal DS object based on the current
244:    settings of the eigenvalue solver.

246:    Collective

248:    Input Parameter:
249: .  eps - eigenproblem solver context

251:    Note:
252:    This function need not be called explicitly, since it will be called at
253:    both EPSSetFromOptions() and EPSSetUp().

255:    Level: developer

257: .seealso: EPSSetFromOptions(), EPSSetUp()
258: @*/
259: PetscErrorCode EPSSetDSType(EPS eps)
260: {
261:   PetscFunctionBegin;
263:   PetscTryTypeMethod(eps,setdstype);
264:   PetscFunctionReturn(PETSC_SUCCESS);
265: }

267: /*@
268:    EPSSetUp - Sets up all the internal data structures necessary for the
269:    execution of the eigensolver. Then calls STSetUp() for any set-up
270:    operations associated to the ST object.

272:    Collective

274:    Input Parameter:
275: .  eps   - eigenproblem solver context

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

282:    Level: developer

284: .seealso: EPSCreate(), EPSSolve(), EPSDestroy(), STSetUp(), EPSSetInitialSpace()
285: @*/
286: PetscErrorCode EPSSetUp(EPS eps)
287: {
288:   Mat            A,B;
289:   PetscInt       k,nmat;
290:   PetscBool      flg;

292:   PetscFunctionBegin;
294:   if (eps->state) PetscFunctionReturn(PETSC_SUCCESS);
295:   PetscCall(PetscLogEventBegin(EPS_SetUp,eps,0,0,0));

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

300:   /* Set default solver type (EPSSetFromOptions was not called) */
301:   if (!((PetscObject)eps)->type_name) PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
302:   if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
303:   PetscCall(EPSSetDefaultST(eps));

305:   PetscCall(STSetTransform(eps->st,PETSC_TRUE));
306:   PetscCall(STSetStructured(eps->st,PETSC_FALSE));
307:   if (eps->useds && !eps->ds) PetscCall(EPSGetDS(eps,&eps->ds));
308:   if (eps->useds) PetscCall(EPSSetDSType(eps));
309:   if (eps->twosided) {
310:     PetscCheck(!eps->ishermitian || (eps->isgeneralized && !eps->ispositive),PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided methods are not intended for %s problems",SLEPC_STRING_HERMITIAN);
311:   }
312:   if (!eps->rg) PetscCall(EPSGetRG(eps,&eps->rg));
313:   if (!((PetscObject)eps->rg)->type_name) PetscCall(RGSetType(eps->rg,RGINTERVAL));

315:   /* Set problem dimensions */
316:   PetscCall(STGetNumMatrices(eps->st,&nmat));
317:   PetscCheck(nmat,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSSetOperators must be called first");
318:   PetscCall(STMatGetSize(eps->st,&eps->n,NULL));
319:   PetscCall(STMatGetLocalSize(eps->st,&eps->nloc,NULL));

321:   /* Set default problem type */
322:   if (!eps->problem_type) {
323:     if (nmat==1) PetscCall(EPSSetProblemType(eps,EPS_NHEP));
324:     else PetscCall(EPSSetProblemType(eps,EPS_GNHEP));
325:   } else if (nmat==1 && eps->isgeneralized) {
326:     PetscCall(PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n"));
327:     eps->isgeneralized = PETSC_FALSE;
328:     eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;
329:   } else PetscCheck(nmat==1 || eps->isgeneralized,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Inconsistent EPS state: the problem type does not match the number of matrices");

331:   if (eps->isstructured) {
332:     /* make sure the user has set the appropriate matrix */
333:     PetscCall(STGetMatrix(eps->st,0,&A));
334:     if (eps->problem_type==EPS_BSE) PetscCall(SlepcCheckMatStruct(A,SLEPC_MAT_STRUCT_BSE,NULL));
335:   }

337:   if (eps->nev > eps->n) eps->nev = eps->n;
338:   if (eps->ncv > eps->n) eps->ncv = eps->n;

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

343:   /* initialization of matrix norms */
344:   if (eps->conv==EPS_CONV_NORM) {
345:     if (!eps->nrma) {
346:       PetscCall(STGetMatrix(eps->st,0,&A));
347:       PetscCall(MatNorm(A,NORM_INFINITY,&eps->nrma));
348:     }
349:     if (nmat>1 && !eps->nrmb) {
350:       PetscCall(STGetMatrix(eps->st,1,&B));
351:       PetscCall(MatNorm(B,NORM_INFINITY,&eps->nrmb));
352:     }
353:   }

355:   /* call specific solver setup */
356:   PetscUseTypeMethod(eps,setup);

358:   /* if purification is set, check that it really makes sense */
359:   if (eps->purify) {
360:     if (eps->categ==EPS_CATEGORY_PRECOND || eps->categ==EPS_CATEGORY_CONTOUR) eps->purify = PETSC_FALSE;
361:     else {
362:       if (!eps->isgeneralized) eps->purify = PETSC_FALSE;
363:       else if (!eps->ishermitian && !eps->ispositive) eps->purify = PETSC_FALSE;
364:       else {
365:         PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg));
366:         if (flg) eps->purify = PETSC_FALSE;
367:       }
368:     }
369:   }

371:   /* set tolerance if not yet set */
372:   if (eps->tol==(PetscReal)PETSC_DETERMINE) eps->tol = SLEPC_DEFAULT_TOL;

374:   /* set up sorting criterion */
375:   PetscTryTypeMethod(eps,setupsort);

377:   /* Build balancing matrix if required */
378:   if (eps->balance!=EPS_BALANCE_USER) {
379:     PetscCall(STSetBalanceMatrix(eps->st,NULL));
380:     if (!eps->ishermitian && (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE)) {
381:       if (!eps->D) PetscCall(BVCreateVec(eps->V,&eps->D));
382:       PetscCall(EPSBuildBalance_Krylov(eps));
383:       PetscCall(STSetBalanceMatrix(eps->st,eps->D));
384:     }
385:   }

387:   /* Setup ST */
388:   PetscCall(STSetUp(eps->st));
389:   PetscCall(EPSCheckCompatibleST(eps));

391:   /* process deflation and initial vectors */
392:   if (eps->nds<0) {
393:     PetscCheck(!eps->isstructured,PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Deflation space is not supported in structured eigenproblems");
394:     k = -eps->nds;
395:     PetscCall(BVInsertConstraints(eps->V,&k,eps->defl));
396:     PetscCall(SlepcBasisDestroy_Private(&eps->nds,&eps->defl));
397:     eps->nds = k;
398:     PetscCall(STCheckNullSpace(eps->st,eps->V));
399:   }
400:   if (eps->nini<0) {
401:     k = -eps->nini;
402:     PetscCheck(k<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
403:     PetscCall(BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE));
404:     PetscCall(SlepcBasisDestroy_Private(&eps->nini,&eps->IS));
405:     eps->nini = k;
406:   }
407:   if (eps->twosided && eps->ninil<0) {
408:     k = -eps->ninil;
409:     PetscCheck(k<=eps->ncv,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The number of left initial vectors is larger than ncv");
410:     PetscCall(BVInsertVecs(eps->W,0,&k,eps->ISL,PETSC_TRUE));
411:     PetscCall(SlepcBasisDestroy_Private(&eps->ninil,&eps->ISL));
412:     eps->ninil = k;
413:   }

415:   PetscCall(PetscLogEventEnd(EPS_SetUp,eps,0,0,0));
416:   eps->state = EPS_STATE_SETUP;
417:   PetscFunctionReturn(PETSC_SUCCESS);
418: }

420: /*@
421:    EPSSetOperators - Sets the matrices associated with the eigenvalue problem.

423:    Collective

425:    Input Parameters:
426: +  eps - the eigenproblem solver context
427: .  A  - the matrix associated with the eigensystem
428: -  B  - the second matrix in the case of generalized eigenproblems

430:    Notes:
431:    To specify a standard eigenproblem, use NULL for parameter B.

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

436:    For structured eigenproblem types such as EPS_BSE (see EPSSetProblemType()), the
437:    provided matrices must have been created with the corresponding helper function,
438:    i.e., MatCreateBSE().

440:    Level: beginner

442: .seealso: EPSSolve(), EPSSetUp(), EPSReset(), EPSGetST(), STGetMatrix(), EPSSetProblemType()
443: @*/
444: PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)
445: {
446:   PetscInt       m,n,m0,mloc,nloc,mloc0,nmat;
447:   Mat            mat[2];

449:   PetscFunctionBegin;
453:   PetscCheckSameComm(eps,1,A,2);
454:   if (B) PetscCheckSameComm(eps,1,B,3);

456:   /* Check matrix sizes */
457:   PetscCall(MatGetSize(A,&m,&n));
458:   PetscCall(MatGetLocalSize(A,&mloc,&nloc));
459:   PetscCheck(m==n,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m,n);
460:   PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A does not have equal row and column sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",mloc,nloc);
461:   if (B) {
462:     PetscCall(MatGetSize(B,&m0,&n));
463:     PetscCall(MatGetLocalSize(B,&mloc0,&nloc));
464:     PetscCheck(m0==n,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",m0,n);
465:     PetscCheck(mloc0==nloc,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",mloc0,nloc);
466:     PetscCheck(m==m0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Dimensions of A and B do not match (%" PetscInt_FMT ", %" PetscInt_FMT ")",m,m0);
467:     PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Local dimensions of A and B do not match (%" PetscInt_FMT ", %" PetscInt_FMT ")",mloc,mloc0);
468:   }
469:   if (eps->state && (n!=eps->n || nloc!=eps->nloc)) PetscCall(EPSReset(eps));
470:   eps->nrma = 0.0;
471:   eps->nrmb = 0.0;
472:   if (!eps->st) PetscCall(EPSGetST(eps,&eps->st));
473:   mat[0] = A;
474:   if (B) {
475:     mat[1] = B;
476:     nmat = 2;
477:   } else nmat = 1;
478:   PetscCall(STSetMatrices(eps->st,nmat,mat));
479:   eps->state = EPS_STATE_INITIAL;
480:   PetscFunctionReturn(PETSC_SUCCESS);
481: }

483: /*@
484:    EPSGetOperators - Gets the matrices associated with the eigensystem.

486:    Collective

488:    Input Parameter:
489: .  eps - the EPS context

491:    Output Parameters:
492: +  A  - the matrix associated with the eigensystem
493: -  B  - the second matrix in the case of generalized eigenproblems

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

498:    Level: intermediate

500: .seealso: EPSSolve(), EPSGetST(), STGetMatrix(), STSetMatrices()
501: @*/
502: PetscErrorCode EPSGetOperators(EPS eps,Mat *A,Mat *B)
503: {
504:   ST             st;
505:   PetscInt       k;

507:   PetscFunctionBegin;
509:   PetscCall(EPSGetST(eps,&st));
510:   PetscCall(STGetNumMatrices(st,&k));
511:   if (A) {
512:     if (k<1) *A = NULL;
513:     else PetscCall(STGetMatrix(st,0,A));
514:   }
515:   if (B) {
516:     if (k<2) *B = NULL;
517:     else PetscCall(STGetMatrix(st,1,B));
518:   }
519:   PetscFunctionReturn(PETSC_SUCCESS);
520: }

522: /*@
523:    EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
524:    space.

526:    Collective

528:    Input Parameters:
529: +  eps - the eigenproblem solver context
530: .  n   - number of vectors
531: -  v   - set of basis vectors of the deflation space

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

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

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

545:    Level: intermediate

547: .seealso: EPSSetInitialSpace()
548: @*/
549: PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec v[])
550: {
551:   PetscFunctionBegin;
554:   PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
555:   if (n>0) {
556:     PetscAssertPointer(v,3);
558:   }
559:   PetscCall(SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl));
560:   if (n>0) eps->state = EPS_STATE_INITIAL;
561:   PetscFunctionReturn(PETSC_SUCCESS);
562: }

564: /*@
565:    EPSSetInitialSpace - Specify a basis of vectors that constitute the initial
566:    space, that is, the subspace from which the solver starts to iterate.

568:    Collective

570:    Input Parameters:
571: +  eps - the eigenproblem solver context
572: .  n   - number of vectors
573: -  is  - set of basis vectors of the initial space

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

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

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

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

588:    Level: intermediate

590: .seealso: EPSSetLeftInitialSpace(), EPSSetDeflationSpace()
591: @*/
592: PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec is[])
593: {
594:   PetscFunctionBegin;
597:   PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
598:   if (n>0) {
599:     PetscAssertPointer(is,3);
601:   }
602:   PetscCall(SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS));
603:   if (n>0) eps->state = EPS_STATE_INITIAL;
604:   PetscFunctionReturn(PETSC_SUCCESS);
605: }

607: /*@
608:    EPSSetLeftInitialSpace - Specify a basis of vectors that constitute the left
609:    initial space, used by two-sided solvers to start the left subspace.

611:    Collective

613:    Input Parameters:
614: +  eps - the eigenproblem solver context
615: .  n   - number of vectors
616: -  isl - set of basis vectors of the left initial space

618:    Notes:
619:    Left initial vectors are used to initiate the left search space in two-sided
620:    eigensolvers. Users should pass here an approximation of the left eigenspace,
621:    if available.

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

625:    Level: intermediate

627: .seealso: EPSSetInitialSpace(), EPSSetTwoSided()
628: @*/
629: PetscErrorCode EPSSetLeftInitialSpace(EPS eps,PetscInt n,Vec isl[])
630: {
631:   PetscFunctionBegin;
634:   PetscCheck(n>=0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
635:   if (n>0) {
636:     PetscAssertPointer(isl,3);
638:   }
639:   PetscCall(SlepcBasisReference_Private(n,isl,&eps->ninil,&eps->ISL));
640:   if (n>0) eps->state = EPS_STATE_INITIAL;
641:   PetscFunctionReturn(PETSC_SUCCESS);
642: }

644: /*
645:   EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
646:   by the user. This is called at setup.
647:  */
648: PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
649: {
650:   PetscBool      krylov;

652:   PetscFunctionBegin;
653:   if (*ncv!=PETSC_DETERMINE) { /* ncv set */
654:     PetscCall(PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,""));
655:     if (krylov) {
656:       PetscCheck(*ncv>=nev+1 || (*ncv==nev && *ncv==eps->n),PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
657:     } else {
658:       PetscCheck(*ncv>=nev,PetscObjectComm((PetscObject)eps),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
659:     }
660:   } else if (*mpd!=PETSC_DETERMINE) { /* mpd set */
661:     *ncv = PetscMin(eps->n,nev+(*mpd));
662:   } else { /* neither set: defaults depend on nev being small or large */
663:     if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
664:     else {
665:       *mpd = 500;
666:       *ncv = PetscMin(eps->n,nev+(*mpd));
667:     }
668:   }
669:   if (*mpd==PETSC_DETERMINE) *mpd = *ncv;
670:   PetscFunctionReturn(PETSC_SUCCESS);
671: }

673: /*@
674:    EPSAllocateSolution - Allocate memory storage for common variables such
675:    as eigenvalues and eigenvectors.

677:    Collective

679:    Input Parameters:
680: +  eps   - eigensolver context
681: -  extra - number of additional positions, used for methods that require a
682:            working basis slightly larger than ncv

684:    Developer Notes:
685:    This is SLEPC_EXTERN because it may be required by user plugin EPS
686:    implementations.

688:    Level: developer

690: .seealso: EPSSetUp()
691: @*/
692: PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
693: {
694:   PetscInt       oldsize,requested;
695:   PetscRandom    rand;
696:   Vec            t;

698:   PetscFunctionBegin;
699:   requested = eps->ncv + extra;

701:   /* oldsize is zero if this is the first time setup is called */
702:   PetscCall(BVGetSizes(eps->V,NULL,NULL,&oldsize));

704:   /* allocate space for eigenvalues and friends */
705:   if (requested != oldsize || !eps->eigr) {
706:     PetscCall(PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm));
707:     PetscCall(PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm));
708:   }

710:   /* workspace for the case of arbitrary selection */
711:   if (eps->arbitrary) {
712:     if (eps->rr) PetscCall(PetscFree2(eps->rr,eps->ri));
713:     PetscCall(PetscMalloc2(requested,&eps->rr,requested,&eps->ri));
714:   }

716:   /* allocate V */
717:   if (!eps->V) PetscCall(EPSGetBV(eps,&eps->V));
718:   if (!oldsize) {
719:     if (!((PetscObject)eps->V)->type_name) PetscCall(BVSetType(eps->V,BVMAT));
720:     PetscCall(STMatCreateVecsEmpty(eps->st,&t,NULL));
721:     PetscCall(BVSetSizesFromVec(eps->V,t,requested));
722:     PetscCall(VecDestroy(&t));
723:   } else PetscCall(BVResize(eps->V,requested,PETSC_FALSE));

725:   /* allocate W */
726:   if (eps->twosided) {
727:     PetscCall(BVGetRandomContext(eps->V,&rand));  /* make sure the random context is available when duplicating */
728:     PetscCall(BVDestroy(&eps->W));
729:     PetscCall(BVDuplicate(eps->V,&eps->W));
730:   }
731:   PetscFunctionReturn(PETSC_SUCCESS);
732: }