Actual source code: pepsetup.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:    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: {

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

 34: /*
 35:    This is used in Q-Arnoldi and STOAR to set the transform flag by
 36:    default, otherwise the user has to explicitly run with -st_transform
 37: */
 38: PetscErrorCode PEPSetDefaultST_Transform(PEP pep)
 39: {

 43:   STSetTransform(pep->st,PETSC_TRUE);
 44:   return(0);
 45: }

 47: /*@
 48:    PEPSetUp - Sets up all the internal data structures necessary for the
 49:    execution of the PEP solver.

 51:    Collective on pep

 53:    Input Parameter:
 54: .  pep   - solver context

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

 61:    Level: developer

 63: .seealso: PEPCreate(), PEPSolve(), PEPDestroy()
 64: @*/
 65: PetscErrorCode PEPSetUp(PEP pep)
 66: {
 68:   SlepcSC        sc;
 69:   PetscBool      istrivial,flg;
 70:   PetscInt       k;
 71:   KSP            ksp;
 72:   PC             pc;
 73:   PetscMPIInt    size;
 74:   MatSolverType  stype;

 78:   if (pep->state) return(0);
 79:   PetscLogEventBegin(PEP_SetUp,pep,0,0,0);

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

 84:   /* set default solver type (PEPSetFromOptions was not called) */
 85:   if (!((PetscObject)pep)->type_name) {
 86:     PEPSetType(pep,PEPTOAR);
 87:   }
 88:   if (!pep->st) { PEPGetST(pep,&pep->st); }
 89:   PEPSetDefaultST(pep);
 90:   if (!pep->ds) { PEPGetDS(pep,&pep->ds); }
 91:   if (!pep->rg) { PEPGetRG(pep,&pep->rg); }
 92:   if (!((PetscObject)pep->rg)->type_name) {
 93:     RGSetType(pep->rg,RGINTERVAL);
 94:   }

 96:   /* check matrices, transfer them to ST */
 97:   if (!pep->A) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"PEPSetOperators must be called first");
 98:   STSetMatrices(pep->st,pep->nmat,pep->A);

100:   /* set problem dimensions */
101:   MatGetSize(pep->A[0],&pep->n,NULL);
102:   MatGetLocalSize(pep->A[0],&pep->nloc,NULL);

104:   /* set default problem type */
105:   if (!pep->problem_type) {
106:     PEPSetProblemType(pep,PEP_GENERAL);
107:   }
108:   if (pep->nev > (pep->nmat-1)*pep->n) pep->nev = (pep->nmat-1)*pep->n;
109:   if (pep->ncv > (pep->nmat-1)*pep->n) pep->ncv = (pep->nmat-1)*pep->n;

111:   /* check consistency of refinement options */
112:   if (pep->refine) {
113:     if (!pep->scheme) {  /* set default scheme */
114:       PEPRefineGetKSP(pep,&ksp);
115:       KSPGetPC(ksp,&pc);
116:       PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
117:       if (flg) {
118:         PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
119:       }
120:       pep->scheme = flg? PEP_REFINE_SCHEME_MBE: PEP_REFINE_SCHEME_SCHUR;
121:     }
122:     if (pep->scheme==PEP_REFINE_SCHEME_MBE) {
123:       PEPRefineGetKSP(pep,&ksp);
124:       KSPGetPC(ksp,&pc);
125:       PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
126:       if (flg) {
127:         PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
128:       }
129:       if (!flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The MBE scheme for refinement requires a direct solver in KSP");
130:       MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
131:       if (size>1) {   /* currently selected PC is a factorization */
132:         PCFactorGetMatSolverType(pc,&stype);
133:         PetscStrcmp(stype,MATSOLVERPETSC,&flg);
134:         if (flg) SETERRQ(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");
135:       }
136:     }
137:     if (pep->scheme==PEP_REFINE_SCHEME_SCHUR) {
138:       if (pep->npart>1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The Schur scheme for refinement does not support subcommunicators");
139:     }
140:   }
141:   /* call specific solver setup */
142:   (*pep->ops->setup)(pep);

144:   /* set tolerance if not yet set */
145:   if (pep->tol==PETSC_DEFAULT) pep->tol = SLEPC_DEFAULT_TOL;
146:   if (pep->refine) {
147:     if (pep->rtol==PETSC_DEFAULT) pep->rtol = PetscMax(pep->tol/1000,PETSC_MACHINE_EPSILON);
148:     if (pep->rits==PETSC_DEFAULT) pep->rits = (pep->refine==PEP_REFINE_SIMPLE)? 10: 1;
149:   }

151:   /* set default extraction */
152:   if (!pep->extract) {
153:     pep->extract = (pep->basis==PEP_BASIS_MONOMIAL)? PEP_EXTRACT_NORM: PEP_EXTRACT_NONE;
154:   }

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

206:   /* fill sorting criterion for DS */
207:   if (pep->which!=PEP_ALL) {
208:     DSGetSlepcSC(pep->ds,&sc);
209:     RGIsTrivial(pep->rg,&istrivial);
210:     sc->rg            = istrivial? NULL: pep->rg;
211:     sc->comparison    = pep->sc->comparison;
212:     sc->comparisonctx = pep->sc->comparisonctx;
213:     sc->map           = SlepcMap_ST;
214:     sc->mapobj        = (PetscObject)pep->st;
215:   }
216:   /* setup ST */
217:   STSetUp(pep->st);

219:   /* compute matrix coefficients */
220:   STGetTransform(pep->st,&flg);
221:   if (!flg) {
222:     if (pep->which!=PEP_ALL && pep->solvematcoeffs) { STMatSetUp(pep->st,1.0,pep->solvematcoeffs); }
223:   } else {
224:     if (pep->basis!=PEP_BASIS_MONOMIAL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Cannot use ST-transform with non-monomial basis in PEP");
225:   }

227:   /* compute scale factor if no set by user */
228:   PEPComputeScaleFactor(pep);

230:   /* build balancing matrix if required */
231:   if (pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) {
232:     if (!pep->Dl) {
233:       BVCreateVec(pep->V,&pep->Dl);
234:       PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->Dl);
235:     }
236:     if (!pep->Dr) {
237:       BVCreateVec(pep->V,&pep->Dr);
238:       PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->Dr);
239:     }
240:     PEPBuildDiagonalScaling(pep);
241:   }

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

256: /*@
257:    PEPSetOperators - Sets the coefficient matrices associated with the polynomial
258:    eigenvalue problem.

260:    Collective on pep

262:    Input Parameters:
263: +  pep  - the eigenproblem solver context
264: .  nmat - number of matrices in array A
265: -  A    - the array of matrices associated with the eigenproblem

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

273:    Level: beginner

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

285:   if (nmat <= 0) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Non-positive value of nmat: %D",nmat);
286:   if (nmat <= 2) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Cannot solve linear eigenproblems with PEP; use EPS instead");

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

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

309:   PetscMalloc1(nmat,&pep->A);
310:   PetscCalloc2(3*nmat,&pep->pbc,nmat,&pep->nrma);
311:   PetscLogObjectMemory((PetscObject)pep,nmat*sizeof(Mat)+4*nmat*sizeof(PetscReal));
312:   for (i=0;i<nmat;i++) {
313:     pep->A[i]   = A[i];
314:     pep->pbc[i] = 1.0;  /* default to monomial basis */
315:   }
316:   pep->nmat = nmat;
317:   pep->state = PEP_STATE_INITIAL;
318:   return(0);
319: }

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

324:    Not collective, though parallel Mats are returned if the PEP is parallel

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

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

333:    Level: intermediate

335: .seealso: PEPSolve(), PEPSetOperators(), PEPGetNumMatrices()
336: @*/
337: PetscErrorCode PEPGetOperators(PEP pep,PetscInt k,Mat *A)
338: {
342:   if (k<0 || k>=pep->nmat) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",pep->nmat-1);
343:   *A = pep->A[k];
344:   return(0);
345: }

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

350:    Not collective

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

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

358:    Level: intermediate

360: .seealso: PEPSetOperators()
361: @*/
362: PetscErrorCode PEPGetNumMatrices(PEP pep,PetscInt *nmat)
363: {
367:   *nmat = pep->nmat;
368:   return(0);
369: }

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

375:    Collective on pep

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

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

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

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

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

395:    Level: intermediate
396: @*/
397: PetscErrorCode PEPSetInitialSpace(PEP pep,PetscInt n,Vec is[])
398: {

404:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
405:   if (n>0) {
408:   }
409:   SlepcBasisReference_Private(n,is,&pep->nini,&pep->IS);
410:   if (n>0) pep->state = PEP_STATE_INITIAL;
411:   return(0);
412: }

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

425:   PetscObjectTypeCompareAny((PetscObject)pep,&krylov,PEPTOAR,PEPSTOAR,PEPQARNOLDI,"");
426:   dim = (pep->nmat-1)*pep->n;
427:   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
428:     if (krylov) {
429:       if (*ncv<nev+1 && !(*ncv==nev && *ncv==dim)) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev+1");
430:     } else {
431:       if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
432:     }
433:   } else if (*mpd!=PETSC_DEFAULT) { /* 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_DEFAULT) *mpd = *ncv;
443:   return(0);
444: }

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

450:    Collective on pep

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:    Developers Note:
458:    This is SLEPC_EXTERN because it may be required by user plugin PEP
459:    implementations.

461:    Level: developer
462: @*/
463: PetscErrorCode PEPAllocateSolution(PEP pep,PetscInt extra)
464: {
466:   PetscInt       oldsize,newc,requested,requestedbv;
467:   PetscLogDouble cnt;
468:   Vec            t;

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:   BVGetSizes(pep->V,NULL,NULL,&oldsize);

477:   /* allocate space for eigenvalues and friends */
478:   if (requested != oldsize || !pep->eigr) {
479:     PetscFree4(pep->eigr,pep->eigi,pep->errest,pep->perm);
480:     PetscMalloc4(requested,&pep->eigr,requested,&pep->eigi,requested,&pep->errest,requested,&pep->perm);
481:     newc = PetscMax(0,requested-oldsize);
482:     cnt = 2*newc*sizeof(PetscScalar) + newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
483:     PetscLogObjectMemory((PetscObject)pep,cnt);
484:   }

486:   /* allocate V */
487:   if (!pep->V) { PEPGetBV(pep,&pep->V); }
488:   if (!oldsize) {
489:     if (!((PetscObject)(pep->V))->type_name) {
490:       BVSetType(pep->V,BVSVEC);
491:     }
492:     STMatCreateVecsEmpty(pep->st,&t,NULL);
493:     BVSetSizesFromVec(pep->V,t,requestedbv);
494:     VecDestroy(&t);
495:   } else {
496:     BVResize(pep->V,requestedbv,PETSC_FALSE);
497:   }
498:   return(0);
499: }