Actual source code: pepsetup.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:    PEP routines related to problem setup
 12: */

 14: #include <slepc/private/pepimpl.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 PEPSetFromOptions (before STSetFromOptions)
 20:    and also at PEPSetUp (in case PEPSetFromOptions was not called).
 21: */
 22: PetscErrorCode PEPSetDefaultST(PEP pep)
 23: {
 24:   PetscFunctionBegin;
 25:   PetscTryTypeMethod(pep,setdefaultst);
 26:   if (!((PetscObject)pep->st)->type_name) PetscCall(STSetType(pep->st,STSHIFT));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

 30: /*
 31:    This is used in Q-Arnoldi and STOAR to set the transform flag by
 32:    default, otherwise the user has to explicitly run with -st_transform
 33: */
 34: PetscErrorCode PEPSetDefaultST_Transform(PEP pep)
 35: {
 36:   PetscFunctionBegin;
 37:   PetscCall(STSetTransform(pep->st,PETSC_TRUE));
 38:   PetscFunctionReturn(PETSC_SUCCESS);
 39: }

 41: /*@
 42:    PEPSetDSType - Sets the type of the internal DS object based on the current
 43:    settings of the polynomial eigensolver.

 45:    Collective

 47:    Input Parameter:
 48: .  pep - polynomial eigensolver context

 50:    Note:
 51:    This function need not be called explicitly, since it will be called at
 52:    both PEPSetFromOptions() and PEPSetUp().

 54:    Level: developer

 56: .seealso: PEPSetFromOptions(), PEPSetUp()
 57: @*/
 58: PetscErrorCode PEPSetDSType(PEP pep)
 59: {
 60:   PetscFunctionBegin;
 62:   PetscTryTypeMethod(pep,setdstype);
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: /*@
 67:    PEPSetUp - Sets up all the internal data structures necessary for the
 68:    execution of the PEP solver.

 70:    Collective

 72:    Input Parameter:
 73: .  pep   - solver context

 75:    Notes:
 76:    This function need not be called explicitly in most cases, since PEPSolve()
 77:    calls it. It can be useful when one wants to measure the set-up time
 78:    separately from the solve time.

 80:    Level: developer

 82: .seealso: PEPCreate(), PEPSolve(), PEPDestroy()
 83: @*/
 84: PetscErrorCode PEPSetUp(PEP pep)
 85: {
 86:   SlepcSC        sc;
 87:   PetscBool      istrivial,flg;
 88:   PetscInt       k;
 89:   KSP            ksp;
 90:   PC             pc;
 91:   PetscMPIInt    size;
 92:   MatSolverType  stype;

 94:   PetscFunctionBegin;
 96:   if (pep->state) PetscFunctionReturn(PETSC_SUCCESS);
 97:   PetscCall(PetscLogEventBegin(PEP_SetUp,pep,0,0,0));

 99:   /* reset the convergence flag from the previous solves */
100:   pep->reason = PEP_CONVERGED_ITERATING;

102:   /* set default solver type (PEPSetFromOptions was not called) */
103:   if (!((PetscObject)pep)->type_name) PetscCall(PEPSetType(pep,PEPTOAR));
104:   if (!pep->st) PetscCall(PEPGetST(pep,&pep->st));
105:   PetscCall(PEPSetDefaultST(pep));
106:   if (!pep->ds) PetscCall(PEPGetDS(pep,&pep->ds));
107:   PetscCall(PEPSetDSType(pep));
108:   if (!pep->rg) PetscCall(PEPGetRG(pep,&pep->rg));
109:   if (!((PetscObject)pep->rg)->type_name) PetscCall(RGSetType(pep->rg,RGINTERVAL));

111:   /* check matrices, transfer them to ST */
112:   PetscCheck(pep->A,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"PEPSetOperators must be called first");
113:   PetscCall(STSetMatrices(pep->st,pep->nmat,pep->A));

115:   /* set problem dimensions */
116:   PetscCall(MatGetSize(pep->A[0],&pep->n,NULL));
117:   PetscCall(MatGetLocalSize(pep->A[0],&pep->nloc,NULL));

119:   /* set default problem type */
120:   if (!pep->problem_type) PetscCall(PEPSetProblemType(pep,PEP_GENERAL));
121:   if (pep->nev > (pep->nmat-1)*pep->n) pep->nev = (pep->nmat-1)*pep->n;
122:   if (pep->ncv > (pep->nmat-1)*pep->n) pep->ncv = (pep->nmat-1)*pep->n;

124:   /* check consistency of refinement options */
125:   if (pep->refine) {
126:     if (!pep->scheme) {  /* set default scheme */
127:       PetscCall(PEPRefineGetKSP(pep,&ksp));
128:       PetscCall(KSPGetPC(ksp,&pc));
129:       PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg));
130:       if (flg) PetscCall(PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,""));
131:       pep->scheme = flg? PEP_REFINE_SCHEME_MBE: PEP_REFINE_SCHEME_SCHUR;
132:     }
133:     if (pep->scheme==PEP_REFINE_SCHEME_MBE) {
134:       PetscCall(PEPRefineGetKSP(pep,&ksp));
135:       PetscCall(KSPGetPC(ksp,&pc));
136:       PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg));
137:       if (flg) PetscCall(PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,""));
138:       PetscCheck(flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The MBE scheme for refinement requires a direct solver in KSP");
139:       PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size));
140:       if (size>1) {   /* currently selected PC is a factorization */
141:         PetscCall(PCFactorGetMatSolverType(pc,&stype));
142:         PetscCall(PetscStrcmp(stype,MATSOLVERPETSC,&flg));
143:         PetscCheck(!flg,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"For Newton refinement, you chose to solve linear systems with a factorization, but in parallel runs you need to select an external package");
144:       }
145:     }
146:     if (pep->scheme==PEP_REFINE_SCHEME_SCHUR) {
147:       PetscCheck(pep->npart==1,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The Schur scheme for refinement does not support subcommunicators");
148:     }
149:   }
150:   /* call specific solver setup */
151:   PetscUseTypeMethod(pep,setup);

153:   /* set tolerance if not yet set */
154:   if (pep->tol==(PetscReal)PETSC_DETERMINE) pep->tol = SLEPC_DEFAULT_TOL;
155:   if (pep->refine) {
156:     if (pep->rtol==(PetscReal)PETSC_DETERMINE) pep->rtol = PetscMax(pep->tol/1000,PETSC_MACHINE_EPSILON);
157:     if (pep->rits==PETSC_DETERMINE) pep->rits = (pep->refine==PEP_REFINE_SIMPLE)? 10: 1;
158:   }

160:   /* set default extraction */
161:   if (!pep->extract) {
162:     pep->extract = (pep->basis==PEP_BASIS_MONOMIAL)? PEP_EXTRACT_NORM: PEP_EXTRACT_NONE;
163:   }

165:   /* fill sorting criterion context */
166:   switch (pep->which) {
167:     case PEP_LARGEST_MAGNITUDE:
168:       pep->sc->comparison    = SlepcCompareLargestMagnitude;
169:       pep->sc->comparisonctx = NULL;
170:       break;
171:     case PEP_SMALLEST_MAGNITUDE:
172:       pep->sc->comparison    = SlepcCompareSmallestMagnitude;
173:       pep->sc->comparisonctx = NULL;
174:       break;
175:     case PEP_LARGEST_REAL:
176:       pep->sc->comparison    = SlepcCompareLargestReal;
177:       pep->sc->comparisonctx = NULL;
178:       break;
179:     case PEP_SMALLEST_REAL:
180:       pep->sc->comparison    = SlepcCompareSmallestReal;
181:       pep->sc->comparisonctx = NULL;
182:       break;
183:     case PEP_LARGEST_IMAGINARY:
184:       pep->sc->comparison    = SlepcCompareLargestImaginary;
185:       pep->sc->comparisonctx = NULL;
186:       break;
187:     case PEP_SMALLEST_IMAGINARY:
188:       pep->sc->comparison    = SlepcCompareSmallestImaginary;
189:       pep->sc->comparisonctx = NULL;
190:       break;
191:     case PEP_TARGET_MAGNITUDE:
192:       pep->sc->comparison    = SlepcCompareTargetMagnitude;
193:       pep->sc->comparisonctx = &pep->target;
194:       break;
195:     case PEP_TARGET_REAL:
196:       pep->sc->comparison    = SlepcCompareTargetReal;
197:       pep->sc->comparisonctx = &pep->target;
198:       break;
199:     case PEP_TARGET_IMAGINARY:
200: #if defined(PETSC_USE_COMPLEX)
201:       pep->sc->comparison    = SlepcCompareTargetImaginary;
202:       pep->sc->comparisonctx = &pep->target;
203: #endif
204:       break;
205:     case PEP_ALL:
206:       pep->sc->comparison    = SlepcCompareSmallestReal;
207:       pep->sc->comparisonctx = NULL;
208:       break;
209:     case PEP_WHICH_USER:
210:       break;
211:   }
212:   pep->sc->map    = NULL;
213:   pep->sc->mapobj = NULL;

215:   /* fill sorting criterion for DS */
216:   if (pep->which!=PEP_ALL) {
217:     PetscCall(DSGetSlepcSC(pep->ds,&sc));
218:     PetscCall(RGIsTrivial(pep->rg,&istrivial));
219:     sc->rg            = istrivial? NULL: pep->rg;
220:     sc->comparison    = pep->sc->comparison;
221:     sc->comparisonctx = pep->sc->comparisonctx;
222:     sc->map           = SlepcMap_ST;
223:     sc->mapobj        = (PetscObject)pep->st;
224:   }
225:   /* setup ST */
226:   PetscCall(STSetUp(pep->st));

228:   /* compute matrix coefficients */
229:   PetscCall(STGetTransform(pep->st,&flg));
230:   if (!flg) {
231:     if (pep->which!=PEP_ALL && pep->solvematcoeffs) PetscCall(STMatSetUp(pep->st,1.0,pep->solvematcoeffs));
232:   } else {
233:     PetscCheck(pep->basis==PEP_BASIS_MONOMIAL,PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Cannot use ST-transform with non-monomial basis in PEP");
234:   }

236:   /* compute scale factor if no set by user */
237:   PetscCall(PEPComputeScaleFactor(pep));

239:   /* build balancing matrix if required */
240:   if (pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) {
241:     if (!pep->Dl) PetscCall(BVCreateVec(pep->V,&pep->Dl));
242:     if (!pep->Dr) PetscCall(BVCreateVec(pep->V,&pep->Dr));
243:     PetscCall(PEPBuildDiagonalScaling(pep));
244:   }

246:   /* process initial vectors */
247:   if (pep->nini<0) {
248:     k = -pep->nini;
249:     PetscCheck(k<=pep->ncv,PetscObjectComm((PetscObject)pep),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
250:     PetscCall(BVInsertVecs(pep->V,0,&k,pep->IS,PETSC_TRUE));
251:     PetscCall(SlepcBasisDestroy_Private(&pep->nini,&pep->IS));
252:     pep->nini = k;
253:   }
254:   PetscCall(PetscLogEventEnd(PEP_SetUp,pep,0,0,0));
255:   pep->state = PEP_STATE_SETUP;
256:   PetscFunctionReturn(PETSC_SUCCESS);
257: }

259: /*@
260:    PEPSetOperators - Sets the coefficient matrices associated with the polynomial
261:    eigenvalue problem.

263:    Collective

265:    Input Parameters:
266: +  pep  - the eigenproblem solver context
267: .  nmat - number of matrices in array A
268: -  A    - the array of matrices associated with the eigenproblem

270:    Notes:
271:    The polynomial eigenproblem is defined as P(l)*x=0, where l is
272:    the eigenvalue, x is the eigenvector, and P(l) is defined as
273:    P(l) = A_0 + l*A_1 + ... + l^d*A_d, with d=nmat-1 (the degree of P).
274:    For non-monomial bases, this expression is different.

276:    Level: beginner

278: .seealso: PEPSolve(), PEPGetOperators(), PEPGetNumMatrices(), PEPSetBasis()
279: @*/
280: PetscErrorCode PEPSetOperators(PEP pep,PetscInt nmat,Mat A[])
281: {
282:   PetscInt       i,n=0,m,m0=0,mloc,nloc,mloc0=0;

284:   PetscFunctionBegin;
287:   PetscCheck(nmat>0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Non-positive value of nmat: %" PetscInt_FMT,nmat);
288:   PetscCheck(nmat>2,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Cannot solve linear eigenproblems with PEP; use EPS instead");
289:   PetscAssertPointer(A,3);

291:   for (i=0;i<nmat;i++) {
293:     PetscCheckSameComm(pep,1,A[i],3);
294:     PetscCall(MatGetSize(A[i],&m,&n));
295:     PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
296:     PetscCheck(m==n,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,m,n);
297:     PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
298:     if (!i) { m0 = m; mloc0 = mloc; }
299:     PetscCheck(m==m0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_INCOMP,"Dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,m,m0);
300:     PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_INCOMP,"Local dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
301:     PetscCall(PetscObjectReference((PetscObject)A[i]));
302:   }

304:   if (pep->state && (n!=pep->n || nloc!=pep->nloc)) PetscCall(PEPReset(pep));
305:   else if (pep->nmat) {
306:     PetscCall(MatDestroyMatrices(pep->nmat,&pep->A));
307:     PetscCall(PetscFree2(pep->pbc,pep->nrma));
308:     PetscCall(PetscFree(pep->solvematcoeffs));
309:   }

311:   PetscCall(PetscMalloc1(nmat,&pep->A));
312:   PetscCall(PetscCalloc2(3*nmat,&pep->pbc,nmat,&pep->nrma));
313:   for (i=0;i<nmat;i++) {
314:     pep->A[i]   = A[i];
315:     pep->pbc[i] = 1.0;  /* default to monomial basis */
316:   }
317:   pep->nmat = nmat;
318:   pep->state = PEP_STATE_INITIAL;
319:   PetscFunctionReturn(PETSC_SUCCESS);
320: }

322: /*@
323:    PEPGetOperators - Gets the matrices associated with the polynomial eigensystem.

325:    Not Collective

327:    Input Parameters:
328: +  pep - the PEP context
329: -  k   - the index of the requested matrix (starting in 0)

331:    Output Parameter:
332: .  A - the requested matrix

334:    Level: intermediate

336: .seealso: PEPSolve(), PEPSetOperators(), PEPGetNumMatrices()
337: @*/
338: PetscErrorCode PEPGetOperators(PEP pep,PetscInt k,Mat *A)
339: {
340:   PetscFunctionBegin;
342:   PetscAssertPointer(A,3);
343:   PetscCheck(k>=0 && k<pep->nmat,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,pep->nmat-1);
344:   *A = pep->A[k];
345:   PetscFunctionReturn(PETSC_SUCCESS);
346: }

348: /*@
349:    PEPGetNumMatrices - Returns the number of matrices stored in the PEP.

351:    Not Collective

353:    Input Parameter:
354: .  pep - the PEP context

356:    Output Parameters:
357: .  nmat - the number of matrices passed in PEPSetOperators()

359:    Level: intermediate

361: .seealso: PEPSetOperators()
362: @*/
363: PetscErrorCode PEPGetNumMatrices(PEP pep,PetscInt *nmat)
364: {
365:   PetscFunctionBegin;
367:   PetscAssertPointer(nmat,2);
368:   *nmat = pep->nmat;
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }

372: /*@
373:    PEPSetInitialSpace - Specify a basis of vectors that constitute the initial
374:    space, that is, the subspace from which the solver starts to iterate.

376:    Collective

378:    Input Parameters:
379: +  pep   - the polynomial eigensolver context
380: .  n     - number of vectors
381: -  is    - set of basis vectors of the initial space

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

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

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

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

396:    Level: intermediate

398: .seealso: PEPSetUp()
399: @*/
400: PetscErrorCode PEPSetInitialSpace(PEP pep,PetscInt n,Vec is[])
401: {
402:   PetscFunctionBegin;
405:   PetscCheck(n>=0,PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
406:   if (n>0) {
407:     PetscAssertPointer(is,3);
409:   }
410:   PetscCall(SlepcBasisReference_Private(n,is,&pep->nini,&pep->IS));
411:   if (n>0) pep->state = PEP_STATE_INITIAL;
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }

415: /*
416:   PEPSetDimensions_Default - Set reasonable values for ncv, mpd if not set
417:   by the user. This is called at setup.
418:  */
419: PetscErrorCode PEPSetDimensions_Default(PEP pep,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
420: {
421:   PetscBool      krylov;
422:   PetscInt       dim;

424:   PetscFunctionBegin;
425:   PetscCall(PetscObjectTypeCompareAny((PetscObject)pep,&krylov,PEPTOAR,PEPSTOAR,PEPQARNOLDI,""));
426:   dim = (pep->nmat-1)*pep->n;
427:   if (*ncv!=PETSC_DETERMINE) { /* ncv set */
428:     if (krylov) {
429:       PetscCheck(*ncv>nev || (*ncv==nev && *ncv==dim),PetscObjectComm((PetscObject)pep),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
430:     } else {
431:       PetscCheck(*ncv>=nev,PetscObjectComm((PetscObject)pep),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
432:     }
433:   } else if (*mpd!=PETSC_DETERMINE) { /* mpd set */
434:     *ncv = PetscMin(dim,nev+(*mpd));
435:   } else { /* neither set: defaults depend on nev being small or large */
436:     if (nev<500) *ncv = PetscMin(dim,PetscMax(2*nev,nev+15));
437:     else {
438:       *mpd = 500;
439:       *ncv = PetscMin(dim,nev+(*mpd));
440:     }
441:   }
442:   if (*mpd==PETSC_DETERMINE) *mpd = *ncv;
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: /*@
447:    PEPAllocateSolution - Allocate memory storage for common variables such
448:    as eigenvalues and eigenvectors.

450:    Collective

452:    Input Parameters:
453: +  pep   - eigensolver context
454: -  extra - number of additional positions, used for methods that require a
455:            working basis slightly larger than ncv

457:    Developer Notes:
458:    This is SLEPC_EXTERN because it may be required by user plugin PEP
459:    implementations.

461:    Level: developer

463: .seealso: PEPSetUp()
464: @*/
465: PetscErrorCode PEPAllocateSolution(PEP pep,PetscInt extra)
466: {
467:   PetscInt       oldsize,requested,requestedbv;
468:   Vec            t;

470:   PetscFunctionBegin;
471:   requested = (pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1)) + extra;
472:   requestedbv = pep->ncv + extra;

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

477:   /* allocate space for eigenvalues and friends */
478:   if (requested != oldsize || !pep->eigr) {
479:     PetscCall(PetscFree4(pep->eigr,pep->eigi,pep->errest,pep->perm));
480:     PetscCall(PetscMalloc4(requested,&pep->eigr,requested,&pep->eigi,requested,&pep->errest,requested,&pep->perm));
481:   }

483:   /* allocate V */
484:   if (!pep->V) PetscCall(PEPGetBV(pep,&pep->V));
485:   if (!oldsize) {
486:     if (!((PetscObject)pep->V)->type_name) PetscCall(BVSetType(pep->V,BVMAT));
487:     PetscCall(STMatCreateVecsEmpty(pep->st,&t,NULL));
488:     PetscCall(BVSetSizesFromVec(pep->V,t,requestedbv));
489:     PetscCall(VecDestroy(&t));
490:   } else PetscCall(BVResize(pep->V,requestedbv,PETSC_FALSE));
491:   PetscFunctionReturn(PETSC_SUCCESS);
492: }