Actual source code: nepsetup.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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 NEP

282:    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 NEP
289:    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: }