Actual source code: pepsetup.c

slepc-3.18.2 2023-01-26
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:   PetscTryTypeMethod(pep,setdefaultst);
 25:   if (!((PetscObject)pep->st)->type_name) STSetType(pep->st,STSHIFT);
 26:   return 0;
 27: }

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

 39: /*@
 40:    PEPSetUp - Sets up all the internal data structures necessary for the
 41:    execution of the PEP solver.

 43:    Collective on pep

 45:    Input Parameter:
 46: .  pep   - solver context

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

 53:    Level: developer

 55: .seealso: PEPCreate(), PEPSolve(), PEPDestroy()
 56: @*/
 57: PetscErrorCode PEPSetUp(PEP pep)
 58: {
 59:   SlepcSC        sc;
 60:   PetscBool      istrivial,flg;
 61:   PetscInt       k;
 62:   KSP            ksp;
 63:   PC             pc;
 64:   PetscMPIInt    size;
 65:   MatSolverType  stype;

 68:   if (pep->state) return 0;
 69:   PetscLogEventBegin(PEP_SetUp,pep,0,0,0);

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

 74:   /* set default solver type (PEPSetFromOptions was not called) */
 75:   if (!((PetscObject)pep)->type_name) PEPSetType(pep,PEPTOAR);
 76:   if (!pep->st) PEPGetST(pep,&pep->st);
 77:   PEPSetDefaultST(pep);
 78:   if (!pep->ds) PEPGetDS(pep,&pep->ds);
 79:   if (!pep->rg) PEPGetRG(pep,&pep->rg);
 80:   if (!((PetscObject)pep->rg)->type_name) RGSetType(pep->rg,RGINTERVAL);

 82:   /* check matrices, transfer them to ST */
 84:   STSetMatrices(pep->st,pep->nmat,pep->A);

 86:   /* set problem dimensions */
 87:   MatGetSize(pep->A[0],&pep->n,NULL);
 88:   MatGetLocalSize(pep->A[0],&pep->nloc,NULL);

 90:   /* set default problem type */
 91:   if (!pep->problem_type) PEPSetProblemType(pep,PEP_GENERAL);
 92:   if (pep->nev > (pep->nmat-1)*pep->n) pep->nev = (pep->nmat-1)*pep->n;
 93:   if (pep->ncv > (pep->nmat-1)*pep->n) pep->ncv = (pep->nmat-1)*pep->n;

 95:   /* check consistency of refinement options */
 96:   if (pep->refine) {
 97:     if (!pep->scheme) {  /* set default scheme */
 98:       PEPRefineGetKSP(pep,&ksp);
 99:       KSPGetPC(ksp,&pc);
100:       PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
101:       if (flg) PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
102:       pep->scheme = flg? PEP_REFINE_SCHEME_MBE: PEP_REFINE_SCHEME_SCHUR;
103:     }
104:     if (pep->scheme==PEP_REFINE_SCHEME_MBE) {
105:       PEPRefineGetKSP(pep,&ksp);
106:       KSPGetPC(ksp,&pc);
107:       PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
108:       if (flg) PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
110:       MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
111:       if (size>1) {   /* currently selected PC is a factorization */
112:         PCFactorGetMatSolverType(pc,&stype);
113:         PetscStrcmp(stype,MATSOLVERPETSC,&flg);
115:       }
116:     }
117:     if (pep->scheme==PEP_REFINE_SCHEME_SCHUR) {
119:     }
120:   }
121:   /* call specific solver setup */
122:   PetscUseTypeMethod(pep,setup);

124:   /* set tolerance if not yet set */
125:   if (pep->tol==PETSC_DEFAULT) pep->tol = SLEPC_DEFAULT_TOL;
126:   if (pep->refine) {
127:     if (pep->rtol==PETSC_DEFAULT) pep->rtol = PetscMax(pep->tol/1000,PETSC_MACHINE_EPSILON);
128:     if (pep->rits==PETSC_DEFAULT) pep->rits = (pep->refine==PEP_REFINE_SIMPLE)? 10: 1;
129:   }

131:   /* set default extraction */
132:   if (!pep->extract) {
133:     pep->extract = (pep->basis==PEP_BASIS_MONOMIAL)? PEP_EXTRACT_NORM: PEP_EXTRACT_NONE;
134:   }

136:   /* fill sorting criterion context */
137:   switch (pep->which) {
138:     case PEP_LARGEST_MAGNITUDE:
139:       pep->sc->comparison    = SlepcCompareLargestMagnitude;
140:       pep->sc->comparisonctx = NULL;
141:       break;
142:     case PEP_SMALLEST_MAGNITUDE:
143:       pep->sc->comparison    = SlepcCompareSmallestMagnitude;
144:       pep->sc->comparisonctx = NULL;
145:       break;
146:     case PEP_LARGEST_REAL:
147:       pep->sc->comparison    = SlepcCompareLargestReal;
148:       pep->sc->comparisonctx = NULL;
149:       break;
150:     case PEP_SMALLEST_REAL:
151:       pep->sc->comparison    = SlepcCompareSmallestReal;
152:       pep->sc->comparisonctx = NULL;
153:       break;
154:     case PEP_LARGEST_IMAGINARY:
155:       pep->sc->comparison    = SlepcCompareLargestImaginary;
156:       pep->sc->comparisonctx = NULL;
157:       break;
158:     case PEP_SMALLEST_IMAGINARY:
159:       pep->sc->comparison    = SlepcCompareSmallestImaginary;
160:       pep->sc->comparisonctx = NULL;
161:       break;
162:     case PEP_TARGET_MAGNITUDE:
163:       pep->sc->comparison    = SlepcCompareTargetMagnitude;
164:       pep->sc->comparisonctx = &pep->target;
165:       break;
166:     case PEP_TARGET_REAL:
167:       pep->sc->comparison    = SlepcCompareTargetReal;
168:       pep->sc->comparisonctx = &pep->target;
169:       break;
170:     case PEP_TARGET_IMAGINARY:
171: #if defined(PETSC_USE_COMPLEX)
172:       pep->sc->comparison    = SlepcCompareTargetImaginary;
173:       pep->sc->comparisonctx = &pep->target;
174: #endif
175:       break;
176:     case PEP_ALL:
177:       pep->sc->comparison    = SlepcCompareSmallestReal;
178:       pep->sc->comparisonctx = NULL;
179:       break;
180:     case PEP_WHICH_USER:
181:       break;
182:   }
183:   pep->sc->map    = NULL;
184:   pep->sc->mapobj = NULL;

186:   /* fill sorting criterion for DS */
187:   if (pep->which!=PEP_ALL) {
188:     DSGetSlepcSC(pep->ds,&sc);
189:     RGIsTrivial(pep->rg,&istrivial);
190:     sc->rg            = istrivial? NULL: pep->rg;
191:     sc->comparison    = pep->sc->comparison;
192:     sc->comparisonctx = pep->sc->comparisonctx;
193:     sc->map           = SlepcMap_ST;
194:     sc->mapobj        = (PetscObject)pep->st;
195:   }
196:   /* setup ST */
197:   STSetUp(pep->st);

199:   /* compute matrix coefficients */
200:   STGetTransform(pep->st,&flg);
201:   if (!flg) {
202:     if (pep->which!=PEP_ALL && pep->solvematcoeffs) STMatSetUp(pep->st,1.0,pep->solvematcoeffs);
203:   } else {
205:   }

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

210:   /* build balancing matrix if required */
211:   if (pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) {
212:     if (!pep->Dl) BVCreateVec(pep->V,&pep->Dl);
213:     if (!pep->Dr) BVCreateVec(pep->V,&pep->Dr);
214:     PEPBuildDiagonalScaling(pep);
215:   }

217:   /* process initial vectors */
218:   if (pep->nini<0) {
219:     k = -pep->nini;
221:     BVInsertVecs(pep->V,0,&k,pep->IS,PETSC_TRUE);
222:     SlepcBasisDestroy_Private(&pep->nini,&pep->IS);
223:     pep->nini = k;
224:   }
225:   PetscLogEventEnd(PEP_SetUp,pep,0,0,0);
226:   pep->state = PEP_STATE_SETUP;
227:   return 0;
228: }

230: /*@
231:    PEPSetOperators - Sets the coefficient matrices associated with the polynomial
232:    eigenvalue problem.

234:    Collective on pep

236:    Input Parameters:
237: +  pep  - the eigenproblem solver context
238: .  nmat - number of matrices in array A
239: -  A    - the array of matrices associated with the eigenproblem

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

247:    Level: beginner

249: .seealso: PEPSolve(), PEPGetOperators(), PEPGetNumMatrices(), PEPSetBasis()
250: @*/
251: PetscErrorCode PEPSetOperators(PEP pep,PetscInt nmat,Mat A[])
252: {
253:   PetscInt       i,n=0,m,m0=0,mloc,nloc,mloc0=0;


261:   for (i=0;i<nmat;i++) {
264:     MatGetSize(A[i],&m,&n);
265:     MatGetLocalSize(A[i],&mloc,&nloc);
268:     if (!i) { m0 = m; mloc0 = mloc; }
271:     PetscObjectReference((PetscObject)A[i]);
272:   }

274:   if (pep->state && (n!=pep->n || nloc!=pep->nloc)) PEPReset(pep);
275:   else if (pep->nmat) {
276:     MatDestroyMatrices(pep->nmat,&pep->A);
277:     PetscFree2(pep->pbc,pep->nrma);
278:     PetscFree(pep->solvematcoeffs);
279:   }

281:   PetscMalloc1(nmat,&pep->A);
282:   PetscCalloc2(3*nmat,&pep->pbc,nmat,&pep->nrma);
283:   for (i=0;i<nmat;i++) {
284:     pep->A[i]   = A[i];
285:     pep->pbc[i] = 1.0;  /* default to monomial basis */
286:   }
287:   pep->nmat = nmat;
288:   pep->state = PEP_STATE_INITIAL;
289:   return 0;
290: }

292: /*@
293:    PEPGetOperators - Gets the matrices associated with the polynomial eigensystem.

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

297:    Input Parameters:
298: +  pep - the PEP context
299: -  k   - the index of the requested matrix (starting in 0)

301:    Output Parameter:
302: .  A - the requested matrix

304:    Level: intermediate

306: .seealso: PEPSolve(), PEPSetOperators(), PEPGetNumMatrices()
307: @*/
308: PetscErrorCode PEPGetOperators(PEP pep,PetscInt k,Mat *A)
309: {
313:   *A = pep->A[k];
314:   return 0;
315: }

317: /*@
318:    PEPGetNumMatrices - Returns the number of matrices stored in the PEP.

320:    Not collective

322:    Input Parameter:
323: .  pep - the PEP context

325:    Output Parameters:
326: .  nmat - the number of matrices passed in PEPSetOperators()

328:    Level: intermediate

330: .seealso: PEPSetOperators()
331: @*/
332: PetscErrorCode PEPGetNumMatrices(PEP pep,PetscInt *nmat)
333: {
336:   *nmat = pep->nmat;
337:   return 0;
338: }

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

344:    Collective on pep

346:    Input Parameters:
347: +  pep   - the polynomial eigensolver context
348: .  n     - number of vectors
349: -  is    - set of basis vectors of the initial space

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

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

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

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

364:    Level: intermediate

366: .seealso: PEPSetUp()
367: @*/
368: PetscErrorCode PEPSetInitialSpace(PEP pep,PetscInt n,Vec is[])
369: {
373:   if (n>0) {
376:   }
377:   SlepcBasisReference_Private(n,is,&pep->nini,&pep->IS);
378:   if (n>0) pep->state = PEP_STATE_INITIAL;
379:   return 0;
380: }

382: /*
383:   PEPSetDimensions_Default - Set reasonable values for ncv, mpd if not set
384:   by the user. This is called at setup.
385:  */
386: PetscErrorCode PEPSetDimensions_Default(PEP pep,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
387: {
388:   PetscBool      krylov;
389:   PetscInt       dim;

391:   PetscObjectTypeCompareAny((PetscObject)pep,&krylov,PEPTOAR,PEPSTOAR,PEPQARNOLDI,"");
392:   dim = (pep->nmat-1)*pep->n;
393:   if (*ncv!=PETSC_DEFAULT) { /* ncv set */
394:     if (krylov) {
396:     } else {
398:     }
399:   } else if (*mpd!=PETSC_DEFAULT) { /* mpd set */
400:     *ncv = PetscMin(dim,nev+(*mpd));
401:   } else { /* neither set: defaults depend on nev being small or large */
402:     if (nev<500) *ncv = PetscMin(dim,PetscMax(2*nev,nev+15));
403:     else {
404:       *mpd = 500;
405:       *ncv = PetscMin(dim,nev+(*mpd));
406:     }
407:   }
408:   if (*mpd==PETSC_DEFAULT) *mpd = *ncv;
409:   return 0;
410: }

412: /*@
413:    PEPAllocateSolution - Allocate memory storage for common variables such
414:    as eigenvalues and eigenvectors.

416:    Collective on pep

418:    Input Parameters:
419: +  pep   - eigensolver context
420: -  extra - number of additional positions, used for methods that require a
421:            working basis slightly larger than ncv

423:    Developer Notes:
424:    This is SLEPC_EXTERN because it may be required by user plugin PEP
425:    implementations.

427:    Level: developer

429: .seealso: PEPSetUp()
430: @*/
431: PetscErrorCode PEPAllocateSolution(PEP pep,PetscInt extra)
432: {
433:   PetscInt       oldsize,requested,requestedbv;
434:   Vec            t;

436:   requested = (pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1)) + extra;
437:   requestedbv = pep->ncv + extra;

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

442:   /* allocate space for eigenvalues and friends */
443:   if (requested != oldsize || !pep->eigr) {
444:     PetscFree4(pep->eigr,pep->eigi,pep->errest,pep->perm);
445:     PetscMalloc4(requested,&pep->eigr,requested,&pep->eigi,requested,&pep->errest,requested,&pep->perm);
446:   }

448:   /* allocate V */
449:   if (!pep->V) PEPGetBV(pep,&pep->V);
450:   if (!oldsize) {
451:     if (!((PetscObject)(pep->V))->type_name) BVSetType(pep->V,BVSVEC);
452:     STMatCreateVecsEmpty(pep->st,&t,NULL);
453:     BVSetSizesFromVec(pep->V,t,requestedbv);
454:     VecDestroy(&t);
455:   } else BVResize(pep->V,requestedbv,PETSC_FALSE);
456:   return 0;
457: }