Actual source code: nepsetup.c

slepc-main 2024-11-15
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:    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: }