Actual source code: pepsetup.c
slepc-3.21.0 2024-03-30
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_DEFAULT) pep->tol = SLEPC_DEFAULT_TOL;
155: if (pep->refine) {
156: if (pep->rtol==(PetscReal)PETSC_DEFAULT) pep->rtol = PetscMax(pep->tol/1000,PETSC_MACHINE_EPSILON);
157: if (pep->rits==PETSC_DEFAULT) 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: /*@C
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_DEFAULT) { /* 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_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: 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: }