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