Actual source code: nepsetup.c
slepc-3.22.1 2024-10-28
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: NEP routines related to problem setup
12: */
14: #include <slepc/private/nepimpl.h>
16: /*@
17: NEPSetDSType - Sets the type of the internal DS object based on the current
18: settings of the nonlinear eigensolver.
20: Collective
22: Input Parameter:
23: . nep - nonlinear eigensolver context
25: Note:
26: This function need not be called explicitly, since it will be called at
27: both NEPSetFromOptions() and NEPSetUp().
29: Level: developer
31: .seealso: NEPSetFromOptions(), NEPSetUp()
32: @*/
33: PetscErrorCode NEPSetDSType(NEP nep)
34: {
35: PetscFunctionBegin;
37: PetscTryTypeMethod(nep,setdstype);
38: PetscFunctionReturn(PETSC_SUCCESS);
39: }
41: /*@
42: NEPSetUp - Sets up all the internal data structures necessary for the
43: execution of the NEP solver.
45: Collective
47: Input Parameter:
48: . nep - solver context
50: Notes:
51: This function need not be called explicitly in most cases, since NEPSolve()
52: calls it. It can be useful when one wants to measure the set-up time
53: separately from the solve time.
55: Level: developer
57: .seealso: NEPCreate(), NEPSolve(), NEPDestroy()
58: @*/
59: PetscErrorCode NEPSetUp(NEP nep)
60: {
61: PetscInt k;
62: SlepcSC sc;
63: Mat T;
64: PetscBool flg;
65: KSP ksp;
66: PC pc;
67: PetscMPIInt size;
68: MatSolverType stype;
70: PetscFunctionBegin;
72: NEPCheckProblem(nep,1);
73: if (nep->state) PetscFunctionReturn(PETSC_SUCCESS);
74: PetscCall(PetscLogEventBegin(NEP_SetUp,nep,0,0,0));
76: /* reset the convergence flag from the previous solves */
77: nep->reason = NEP_CONVERGED_ITERATING;
79: /* set default solver type (NEPSetFromOptions was not called) */
80: if (!((PetscObject)nep)->type_name) PetscCall(NEPSetType(nep,NEPRII));
81: if (nep->useds && !nep->ds) PetscCall(NEPGetDS(nep,&nep->ds));
82: if (nep->useds) PetscCall(NEPSetDSType(nep));
83: if (!nep->rg) PetscCall(NEPGetRG(nep,&nep->rg));
84: if (!((PetscObject)nep->rg)->type_name) PetscCall(RGSetType(nep->rg,RGINTERVAL));
86: /* set problem dimensions */
87: switch (nep->fui) {
88: case NEP_USER_INTERFACE_CALLBACK:
89: PetscCall(NEPGetFunction(nep,&T,NULL,NULL,NULL));
90: PetscCall(MatGetSize(T,&nep->n,NULL));
91: PetscCall(MatGetLocalSize(T,&nep->nloc,NULL));
92: break;
93: case NEP_USER_INTERFACE_SPLIT:
94: PetscCall(MatDuplicate(nep->A[0],MAT_DO_NOT_COPY_VALUES,&nep->function));
95: if (nep->P) PetscCall(MatDuplicate(nep->P[0],MAT_DO_NOT_COPY_VALUES,&nep->function_pre));
96: PetscCall(MatDuplicate(nep->A[0],MAT_DO_NOT_COPY_VALUES,&nep->jacobian));
97: PetscCall(MatGetSize(nep->A[0],&nep->n,NULL));
98: PetscCall(MatGetLocalSize(nep->A[0],&nep->nloc,NULL));
99: break;
100: }
102: /* set default problem type */
103: if (!nep->problem_type) PetscCall(NEPSetProblemType(nep,NEP_GENERAL));
105: /* check consistency of refinement options */
106: if (nep->refine) {
107: PetscCheck(nep->fui==NEP_USER_INTERFACE_SPLIT,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Iterative refinement only implemented in split form");
108: if (!nep->scheme) { /* set default scheme */
109: PetscCall(NEPRefineGetKSP(nep,&ksp));
110: PetscCall(KSPGetPC(ksp,&pc));
111: PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg));
112: if (flg) PetscCall(PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,""));
113: nep->scheme = flg? NEP_REFINE_SCHEME_MBE: NEP_REFINE_SCHEME_SCHUR;
114: }
115: if (nep->scheme==NEP_REFINE_SCHEME_MBE) {
116: PetscCall(NEPRefineGetKSP(nep,&ksp));
117: PetscCall(KSPGetPC(ksp,&pc));
118: PetscCall(PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg));
119: if (flg) PetscCall(PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,""));
120: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The MBE scheme for refinement requires a direct solver in KSP");
121: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size));
122: if (size>1) { /* currently selected PC is a factorization */
123: PetscCall(PCFactorGetMatSolverType(pc,&stype));
124: PetscCall(PetscStrcmp(stype,MATSOLVERPETSC,&flg));
125: PetscCheck(!flg,PetscObjectComm((PetscObject)nep),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");
126: }
127: }
128: if (nep->scheme==NEP_REFINE_SCHEME_SCHUR) {
129: PetscCheck(nep->npart==1,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The Schur scheme for refinement does not support subcommunicators");
130: }
131: }
132: /* call specific solver setup */
133: PetscUseTypeMethod(nep,setup);
135: /* set tolerance if not yet set */
136: if (nep->tol==(PetscReal)PETSC_DETERMINE) nep->tol = SLEPC_DEFAULT_TOL;
137: if (nep->refine) {
138: if (nep->rtol==(PetscReal)PETSC_DETERMINE) nep->rtol = PetscMax(nep->tol/1000,PETSC_MACHINE_EPSILON);
139: if (nep->rits==(PetscReal)PETSC_DETERMINE) nep->rits = (nep->refine==NEP_REFINE_SIMPLE)? 10: 1;
140: }
142: /* fill sorting criterion context */
143: switch (nep->which) {
144: case NEP_LARGEST_MAGNITUDE:
145: nep->sc->comparison = SlepcCompareLargestMagnitude;
146: nep->sc->comparisonctx = NULL;
147: break;
148: case NEP_SMALLEST_MAGNITUDE:
149: nep->sc->comparison = SlepcCompareSmallestMagnitude;
150: nep->sc->comparisonctx = NULL;
151: break;
152: case NEP_LARGEST_REAL:
153: nep->sc->comparison = SlepcCompareLargestReal;
154: nep->sc->comparisonctx = NULL;
155: break;
156: case NEP_SMALLEST_REAL:
157: nep->sc->comparison = SlepcCompareSmallestReal;
158: nep->sc->comparisonctx = NULL;
159: break;
160: case NEP_LARGEST_IMAGINARY:
161: nep->sc->comparison = SlepcCompareLargestImaginary;
162: nep->sc->comparisonctx = NULL;
163: break;
164: case NEP_SMALLEST_IMAGINARY:
165: nep->sc->comparison = SlepcCompareSmallestImaginary;
166: nep->sc->comparisonctx = NULL;
167: break;
168: case NEP_TARGET_MAGNITUDE:
169: nep->sc->comparison = SlepcCompareTargetMagnitude;
170: nep->sc->comparisonctx = &nep->target;
171: break;
172: case NEP_TARGET_REAL:
173: nep->sc->comparison = SlepcCompareTargetReal;
174: nep->sc->comparisonctx = &nep->target;
175: break;
176: case NEP_TARGET_IMAGINARY:
177: #if defined(PETSC_USE_COMPLEX)
178: nep->sc->comparison = SlepcCompareTargetImaginary;
179: nep->sc->comparisonctx = &nep->target;
180: #endif
181: break;
182: case NEP_ALL:
183: nep->sc->comparison = SlepcCompareSmallestReal;
184: nep->sc->comparisonctx = NULL;
185: break;
186: case NEP_WHICH_USER:
187: break;
188: }
190: nep->sc->map = NULL;
191: nep->sc->mapobj = NULL;
193: /* fill sorting criterion for DS */
194: if (nep->useds) {
195: PetscCall(DSGetSlepcSC(nep->ds,&sc));
196: sc->comparison = nep->sc->comparison;
197: sc->comparisonctx = nep->sc->comparisonctx;
198: PetscCall(PetscObjectTypeCompare((PetscObject)nep,NEPNLEIGS,&flg));
199: if (!flg) {
200: sc->map = NULL;
201: sc->mapobj = NULL;
202: }
203: }
204: PetscCheck(nep->nev<=nep->ncv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"nev bigger than ncv");
206: /* process initial vectors */
207: if (nep->nini<0) {
208: k = -nep->nini;
209: PetscCheck(k<=nep->ncv,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The number of initial vectors is larger than ncv");
210: PetscCall(BVInsertVecs(nep->V,0,&k,nep->IS,PETSC_TRUE));
211: PetscCall(SlepcBasisDestroy_Private(&nep->nini,&nep->IS));
212: nep->nini = k;
213: }
214: PetscCall(PetscLogEventEnd(NEP_SetUp,nep,0,0,0));
215: nep->state = NEP_STATE_SETUP;
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: /*@
220: NEPSetInitialSpace - Specify a basis of vectors that constitute the initial
221: space, that is, the subspace from which the solver starts to iterate.
223: Collective
225: Input Parameters:
226: + nep - the nonlinear eigensolver context
227: . n - number of vectors
228: - is - set of basis vectors of the initial space
230: Notes:
231: Some solvers start to iterate on a single vector (initial vector). In that case,
232: the other vectors are ignored.
234: These vectors do not persist from one NEPSolve() call to the other, so the
235: initial space should be set every time.
237: The vectors do not need to be mutually orthonormal, since they are explicitly
238: orthonormalized internally.
240: Common usage of this function is when the user can provide a rough approximation
241: of the wanted eigenspace. Then, convergence may be faster.
243: Level: intermediate
245: .seealso: NEPSetUp()
246: @*/
247: PetscErrorCode NEPSetInitialSpace(NEP nep,PetscInt n,Vec is[])
248: {
249: PetscFunctionBegin;
252: PetscCheck(n>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
253: if (n>0) {
254: PetscAssertPointer(is,3);
256: }
257: PetscCall(SlepcBasisReference_Private(n,is,&nep->nini,&nep->IS));
258: if (n>0) nep->state = NEP_STATE_INITIAL;
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: /*
263: NEPSetDimensions_Default - Set reasonable values for ncv, mpd if not set
264: by the user. This is called at setup.
265: */
266: PetscErrorCode NEPSetDimensions_Default(NEP nep,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
267: {
268: PetscFunctionBegin;
269: if (*ncv!=PETSC_DETERMINE) { /* ncv set */
270: PetscCheck(*ncv>=nev,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER_INPUT,"The value of ncv must be at least nev");
271: } else if (*mpd!=PETSC_DETERMINE) { /* mpd set */
272: *ncv = PetscMin(nep->n,nev+(*mpd));
273: } else { /* neither set: defaults depend on nev being small or large */
274: if (nev<500) *ncv = PetscMin(nep->n,PetscMax(2*nev,nev+15));
275: else {
276: *mpd = 500;
277: *ncv = PetscMin(nep->n,nev+(*mpd));
278: }
279: }
280: if (*mpd==PETSC_DETERMINE) *mpd = *ncv;
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*@
285: NEPAllocateSolution - Allocate memory storage for common variables such
286: as eigenvalues and eigenvectors.
288: Collective
290: Input Parameters:
291: + nep - eigensolver context
292: - extra - number of additional positions, used for methods that require a
293: working basis slightly larger than ncv
295: Developer Notes:
296: This is SLEPC_EXTERN because it may be required by user plugin NEP
297: implementations.
299: Level: developer
301: .seealso: PEPSetUp()
302: @*/
303: PetscErrorCode NEPAllocateSolution(NEP nep,PetscInt extra)
304: {
305: PetscInt oldsize,requested;
306: PetscRandom rand;
307: Mat T;
308: Vec t;
310: PetscFunctionBegin;
311: requested = nep->ncv + extra;
313: /* oldsize is zero if this is the first time setup is called */
314: PetscCall(BVGetSizes(nep->V,NULL,NULL,&oldsize));
316: /* allocate space for eigenvalues and friends */
317: if (requested != oldsize || !nep->eigr) {
318: PetscCall(PetscFree4(nep->eigr,nep->eigi,nep->errest,nep->perm));
319: PetscCall(PetscMalloc4(requested,&nep->eigr,requested,&nep->eigi,requested,&nep->errest,requested,&nep->perm));
320: }
322: /* allocate V */
323: if (!nep->V) PetscCall(NEPGetBV(nep,&nep->V));
324: if (!oldsize) {
325: if (!((PetscObject)nep->V)->type_name) PetscCall(BVSetType(nep->V,BVMAT));
326: if (nep->fui==NEP_USER_INTERFACE_SPLIT) T = nep->A[0];
327: else PetscCall(NEPGetFunction(nep,&T,NULL,NULL,NULL));
328: PetscCall(MatCreateVecsEmpty(T,&t,NULL));
329: PetscCall(BVSetSizesFromVec(nep->V,t,requested));
330: PetscCall(VecDestroy(&t));
331: } else PetscCall(BVResize(nep->V,requested,PETSC_FALSE));
333: /* allocate W */
334: if (nep->twosided) {
335: PetscCall(BVGetRandomContext(nep->V,&rand)); /* make sure the random context is available when duplicating */
336: PetscCall(BVDestroy(&nep->W));
337: PetscCall(BVDuplicate(nep->V,&nep->W));
338: }
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }