1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcnep.h" I*/
16: /*@
17: NEPSetUp - Sets up all the internal data structures necessary for the
18: execution of the NEP solver.
20: Collective on NEP 22: Input Parameter:
23: . nep - solver context
25: Notes:
26: This function need not be called explicitly in most cases, since NEPSolve()
27: calls it. It can be useful when one wants to measure the set-up time
28: separately from the solve time.
30: Level: developer
32: .seealso: NEPCreate(), NEPSolve(), NEPDestroy()
33: @*/
34: PetscErrorCode NEPSetUp(NEP nep) 35: {
37: PetscInt k;
38: SlepcSC sc;
39: Mat T;
40: PetscBool flg;
41: KSP ksp;
42: PC pc;
43: PetscMPIInt size;
44: MatSolverType stype;
48: NEPCheckProblem(nep,1);
49: if (nep->state) return(0);
50: PetscLogEventBegin(NEP_SetUp,nep,0,0,0);
52: /* reset the convergence flag from the previous solves */
53: nep->reason = NEP_CONVERGED_ITERATING;
55: /* set default solver type (NEPSetFromOptions was not called) */
56: if (!((PetscObject)nep)->type_name) {
57: NEPSetType(nep,NEPRII);
58: }
59: if (nep->useds && !nep->ds) { NEPGetDS(nep,&nep->ds); }
60: if (!nep->rg) { NEPGetRG(nep,&nep->rg); }
61: if (!((PetscObject)nep->rg)->type_name) {
62: RGSetType(nep->rg,RGINTERVAL);
63: }
64: if (nep->twosided && !nep->hasts) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"This solver does not support computing left eigenvectors (no two-sided variant)");
66: /* set problem dimensions */
67: switch (nep->fui) {
68: case NEP_USER_INTERFACE_CALLBACK:
69: NEPGetFunction(nep,&T,NULL,NULL,NULL);
70: MatGetSize(T,&nep->n,NULL);
71: MatGetLocalSize(T,&nep->nloc,NULL);
72: break;
73: case NEP_USER_INTERFACE_SPLIT:
74: MatDuplicate(nep->A[0],MAT_DO_NOT_COPY_VALUES,&nep->function);
75: MatDuplicate(nep->A[0],MAT_DO_NOT_COPY_VALUES,&nep->jacobian);
76: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->function);
77: PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->jacobian);
78: MatGetSize(nep->A[0],&nep->n,NULL);
79: MatGetLocalSize(nep->A[0],&nep->nloc,NULL);
80: break;
81: case NEP_USER_INTERFACE_DERIVATIVES:
82: NEPGetDerivatives(nep,&T,NULL,NULL);
83: MatGetSize(T,&nep->n,NULL);
84: MatGetLocalSize(T,&nep->nloc,NULL);
85: break;
86: }
88: /* set default problem type */
89: if (!nep->problem_type) {
90: NEPSetProblemType(nep,NEP_GENERAL);
91: }
93: /* check consistency of refinement options */
94: if (nep->refine) {
95: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Iterative refinement only implemented in split form");
96: if (!nep->scheme) { /* set default scheme */
97: NEPRefineGetKSP(nep,&ksp);
98: KSPGetPC(ksp,&pc);
99: PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
100: if (flg) {
101: PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
102: }
103: nep->scheme = flg? NEP_REFINE_SCHEME_MBE: NEP_REFINE_SCHEME_SCHUR;
104: }
105: if (nep->scheme==NEP_REFINE_SCHEME_MBE) {
106: NEPRefineGetKSP(nep,&ksp);
107: KSPGetPC(ksp,&pc);
108: PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
109: if (flg) {
110: PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
111: }
112: if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The MBE scheme for refinement requires a direct solver in KSP");
113: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
114: if (size>1) { /* currently selected PC is a factorization */
115: PCFactorGetMatSolverType(pc,&stype);
116: PetscStrcmp(stype,MATSOLVERPETSC,&flg);
117: if (flg) SETERRQ(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");
118: }
119: }
120: if (nep->scheme==NEP_REFINE_SCHEME_SCHUR) {
121: if (nep->npart>1) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The Schur scheme for refinement does not support subcommunicators");
122: }
123: }
124: /* call specific solver setup */
125: (*nep->ops->setup)(nep);
127: /* by default, compute eigenvalues close to target */
128: /* nep->target should contain the initial guess for the eigenvalue */
129: if (!nep->which) nep->which = NEP_TARGET_MAGNITUDE;
131: /* set tolerance if not yet set */
132: if (nep->tol==PETSC_DEFAULT) nep->tol = SLEPC_DEFAULT_TOL;
133: if (nep->refine) {
134: if (nep->rtol==PETSC_DEFAULT) nep->rtol = PetscMax(nep->tol/1000,PETSC_MACHINE_EPSILON);
135: if (nep->rits==PETSC_DEFAULT) nep->rits = (nep->refine==NEP_REFINE_SIMPLE)? 10: 1;
136: }
138: /* fill sorting criterion context */
139: switch (nep->which) {
140: case NEP_LARGEST_MAGNITUDE:
141: nep->sc->comparison = SlepcCompareLargestMagnitude;
142: nep->sc->comparisonctx = NULL;
143: break;
144: case NEP_SMALLEST_MAGNITUDE:
145: nep->sc->comparison = SlepcCompareSmallestMagnitude;
146: nep->sc->comparisonctx = NULL;
147: break;
148: case NEP_LARGEST_REAL:
149: nep->sc->comparison = SlepcCompareLargestReal;
150: nep->sc->comparisonctx = NULL;
151: break;
152: case NEP_SMALLEST_REAL:
153: nep->sc->comparison = SlepcCompareSmallestReal;
154: nep->sc->comparisonctx = NULL;
155: break;
156: case NEP_LARGEST_IMAGINARY:
157: nep->sc->comparison = SlepcCompareLargestImaginary;
158: nep->sc->comparisonctx = NULL;
159: break;
160: case NEP_SMALLEST_IMAGINARY:
161: nep->sc->comparison = SlepcCompareSmallestImaginary;
162: nep->sc->comparisonctx = NULL;
163: break;
164: case NEP_TARGET_MAGNITUDE:
165: nep->sc->comparison = SlepcCompareTargetMagnitude;
166: nep->sc->comparisonctx = &nep->target;
167: break;
168: case NEP_TARGET_REAL:
169: nep->sc->comparison = SlepcCompareTargetReal;
170: nep->sc->comparisonctx = &nep->target;
171: break;
172: case NEP_TARGET_IMAGINARY:
173: #if defined(PETSC_USE_COMPLEX)
174: nep->sc->comparison = SlepcCompareTargetImaginary;
175: nep->sc->comparisonctx = &nep->target;
176: #endif
177: break;
178: case NEP_ALL:
179: nep->sc->comparison = SlepcCompareSmallestReal;
180: nep->sc->comparisonctx = NULL;
181: break;
182: case NEP_WHICH_USER:
183: break;
184: }
186: nep->sc->map = NULL;
187: nep->sc->mapobj = NULL;
189: /* fill sorting criterion for DS */
190: if (nep->useds) {
191: DSGetSlepcSC(nep->ds,&sc);
192: sc->comparison = nep->sc->comparison;
193: sc->comparisonctx = nep->sc->comparisonctx;
194: PetscObjectTypeCompare((PetscObject)nep,NEPNLEIGS,&flg);
195: if (!flg) {
196: sc->map = NULL;
197: sc->mapobj = NULL;
198: }
199: }
200: if (nep->nev > nep->ncv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"nev bigger than ncv");
202: /* process initial vectors */
203: if (nep->nini<0) {
204: k = -nep->nini;
205: if (k>nep->ncv) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The number of initial vectors is larger than ncv");
206: BVInsertVecs(nep->V,0,&k,nep->IS,PETSC_TRUE);
207: SlepcBasisDestroy_Private(&nep->nini,&nep->IS);
208: nep->nini = k;
209: }
210: PetscLogEventEnd(NEP_SetUp,nep,0,0,0);
211: nep->state = NEP_STATE_SETUP;
212: return(0);
213: }
215: /*@C
216: NEPSetInitialSpace - Specify a basis of vectors that constitute the initial
217: space, that is, the subspace from which the solver starts to iterate.
219: Collective on NEP and Vec
221: Input Parameter:
222: + nep - the nonlinear eigensolver context
223: . n - number of vectors
224: - is - set of basis vectors of the initial space
226: Notes:
227: Some solvers start to iterate on a single vector (initial vector). In that case,
228: the other vectors are ignored.
230: These vectors do not persist from one NEPSolve() call to the other, so the
231: initial space should be set every time.
233: The vectors do not need to be mutually orthonormal, since they are explicitly
234: orthonormalized internally.
236: Common usage of this function is when the user can provide a rough approximation
237: of the wanted eigenspace. Then, convergence may be faster.
239: Level: intermediate
240: @*/
241: PetscErrorCode NEPSetInitialSpace(NEP nep,PetscInt n,Vec *is)242: {
248: if (n<0) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
249: SlepcBasisReference_Private(n,is,&nep->nini,&nep->IS);
250: if (n>0) nep->state = NEP_STATE_INITIAL;
251: return(0);
252: }
254: /*
255: NEPSetDimensions_Default - Set reasonable values for ncv, mpd if not set
256: by the user. This is called at setup.
257: */
258: PetscErrorCode NEPSetDimensions_Default(NEP nep,PetscInt nev,PetscInt *ncv,PetscInt *mpd)259: {
261: if (*ncv) { /* ncv set */
262: if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)nep),1,"The value of ncv must be at least nev");
263: } else if (*mpd) { /* mpd set */
264: *ncv = PetscMin(nep->n,nev+(*mpd));
265: } else { /* neither set: defaults depend on nev being small or large */
266: if (nev<500) *ncv = PetscMin(nep->n,PetscMax(2*nev,nev+15));
267: else {
268: *mpd = 500;
269: *ncv = PetscMin(nep->n,nev+(*mpd));
270: }
271: }
272: if (!*mpd) *mpd = *ncv;
273: return(0);
274: }
276: /*@
277: NEPAllocateSolution - Allocate memory storage for common variables such
278: as eigenvalues and eigenvectors.
280: Collective on NEP282: Input Parameters:
283: + nep - eigensolver context
284: - extra - number of additional positions, used for methods that require a
285: working basis slightly larger than ncv
287: Developers Note:
288: This is SLEPC_EXTERN because it may be required by user plugin NEP289: implementations.
291: Level: developer
292: @*/
293: PetscErrorCode NEPAllocateSolution(NEP nep,PetscInt extra)294: {
296: PetscInt oldsize,newc,requested;
297: PetscLogDouble cnt;
298: Mat T;
299: Vec t;
302: requested = nep->ncv + extra;
304: /* oldsize is zero if this is the first time setup is called */
305: BVGetSizes(nep->V,NULL,NULL,&oldsize);
306: newc = PetscMax(0,requested-oldsize);
308: /* allocate space for eigenvalues and friends */
309: if (requested != oldsize || !nep->eigr) {
310: PetscFree4(nep->eigr,nep->eigi,nep->errest,nep->perm);
311: PetscMalloc4(requested,&nep->eigr,requested,&nep->eigi,requested,&nep->errest,requested,&nep->perm);
312: cnt = newc*sizeof(PetscScalar) + newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
313: PetscLogObjectMemory((PetscObject)nep,cnt);
314: }
316: /* allocate V */
317: if (!nep->V) { NEPGetBV(nep,&nep->V); }
318: if (!oldsize) {
319: if (!((PetscObject)(nep->V))->type_name) {
320: BVSetType(nep->V,BVSVEC);
321: }
322: if (nep->fui==NEP_USER_INTERFACE_SPLIT) T = nep->A[0];
323: else {
324: NEPGetFunction(nep,&T,NULL,NULL,NULL);
325: }
326: MatCreateVecsEmpty(T,&t,NULL);
327: BVSetSizesFromVec(nep->V,t,requested);
328: VecDestroy(&t);
329: } else {
330: BVResize(nep->V,requested,PETSC_FALSE);
331: }
333: /* allocate W */
334: if (nep->twosided) {
335: BVDestroy(&nep->W);
336: BVDuplicate(nep->V,&nep->W);
337: }
338: return(0);
339: }