Actual source code: pepsetup.c

slepc-3.9.2 2018-07-02
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2018, 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:    PEP routines related to problem setup
 12: */

 14: #include <slepc/private/pepimpl.h>       /*I "slepcpep.h" I*/

 16: /*
 17:    Let the solver choose the ST type that should be used by default,
 18:    otherwise set it to SHIFT.
 19:    This is called at PEPSetFromOptions (before STSetFromOptions)
 20:    and also at PEPSetUp (in case PEPSetFromOptions was not called).
 21: */
 22: PetscErrorCode PEPSetDefaultST(PEP pep)
 23: {

 27:   if (pep->ops->setdefaultst) { (*pep->ops->setdefaultst)(pep); }
 28:   if (!((PetscObject)pep->st)->type_name) {
 29:     STSetType(pep->st,STSHIFT);
 30:   }
 31:   return(0);
 32: }

 34: /*
 35:    This is used in Q-Arnoldi and STOAR to set the transform flag by
 36:    default, otherwise the user has to explicitly run with -st_transform
 37: */
 38: PetscErrorCode PEPSetDefaultST_Transform(PEP pep)
 39: {

 43:   STSetTransform(pep->st,PETSC_TRUE);
 44:   return(0);
 45: }

 47: /*@
 48:    PEPSetUp - Sets up all the internal data structures necessary for the
 49:    execution of the PEP solver.

 51:    Collective on PEP

 53:    Input Parameter:
 54: .  pep   - solver context

 56:    Notes:
 57:    This function need not be called explicitly in most cases, since PEPSolve()
 58:    calls it. It can be useful when one wants to measure the set-up time
 59:    separately from the solve time.

 61:    Level: developer

 63: .seealso: PEPCreate(), PEPSolve(), PEPDestroy()
 64: @*/
 65: PetscErrorCode PEPSetUp(PEP pep)
 66: {
 68:   SlepcSC        sc;
 69:   PetscBool      istrivial,flg;
 70:   PetscInt       k;
 71:   KSP            ksp;
 72:   PC             pc;
 73:   PetscMPIInt    size;
 74:   MatSolverType  stype;

 78:   if (pep->state) return(0);
 79:   PetscLogEventBegin(PEP_SetUp,pep,0,0,0);

 81:   /* reset the convergence flag from the previous solves */
 82:   pep->reason = PEP_CONVERGED_ITERATING;

 84:   /* set default solver type (PEPSetFromOptions was not called) */
 85:   if (!((PetscObject)pep)->type_name) {
 86:     PEPSetType(pep,PEPTOAR);
 87:   }
 88:   if (!pep->st) { PEPGetST(pep,&pep->st); }
 89:   PEPSetDefaultST(pep);
 90:   if (!pep->ds) { PEPGetDS(pep,&pep->ds); }
 91:   if (!pep->rg) { PEPGetRG(pep,&pep->rg); }
 92:   if (!((PetscObject)pep->rg)->type_name) {
 93:     RGSetType(pep->rg,RGINTERVAL);
 94:   }

 96:   /* check matrices, transfer them to ST */
 97:   if (!pep->A) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONGSTATE,"PEPSetOperators must be called first");
 98:   STSetMatrices(pep->st,pep->nmat,pep->A);

100:   /* set problem dimensions */
101:   MatGetSize(pep->A[0],&pep->n,NULL);
102:   MatGetLocalSize(pep->A[0],&pep->nloc,NULL);

104:   /* set default problem type */
105:   if (!pep->problem_type) {
106:     PEPSetProblemType(pep,PEP_GENERAL);
107:   }

109:   /* check consistency of refinement options */
110:   if (pep->refine) {
111:     if (!pep->scheme) {  /* set default scheme */
112:       PEPRefineGetKSP(pep,&ksp);
113:       KSPGetPC(ksp,&pc);
114:       PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
115:       if (flg) {
116:         PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
117:       }
118:       pep->scheme = flg? PEP_REFINE_SCHEME_MBE: PEP_REFINE_SCHEME_SCHUR;
119:     }
120:     if (pep->scheme==PEP_REFINE_SCHEME_MBE) {
121:       PEPRefineGetKSP(pep,&ksp);
122:       KSPGetPC(ksp,&pc);
123:       PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
124:       if (flg) {
125:         PetscObjectTypeCompareAny((PetscObject)pc,&flg,PCLU,PCCHOLESKY,"");
126:       }
127:       if (!flg) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The MBE scheme for refinement requires a direct solver in KSP");
128:       MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
129:       if (size>1) {   /* currently selected PC is a factorization */
130:         PCFactorGetMatSolverType(pc,&stype);
131:         PetscStrcmp(stype,MATSOLVERPETSC,&flg);
132:         if (flg) SETERRQ(PetscObjectComm((PetscObject)pep),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");
133:       }
134:     }
135:     if (pep->scheme==PEP_REFINE_SCHEME_SCHUR) {
136:       if (pep->npart>1) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"The Schur scheme for refinement does not support subcommunicators");
137:     }
138:   }
139:   /* call specific solver setup */
140:   (*pep->ops->setup)(pep);

142:   /* set tolerance if not yet set */
143:   if (pep->tol==PETSC_DEFAULT) pep->tol = SLEPC_DEFAULT_TOL;
144:   if (pep->refine) {
145:     if (pep->rtol==PETSC_DEFAULT) pep->rtol = PetscMax(pep->tol/1000,PETSC_MACHINE_EPSILON);
146:     if (pep->rits==PETSC_DEFAULT) pep->rits = (pep->refine==PEP_REFINE_SIMPLE)? 10: 1;
147:   }

149:   /* set default extraction */
150:   if (!pep->extract) {
151:     pep->extract = (pep->basis==PEP_BASIS_MONOMIAL)? PEP_EXTRACT_NORM: PEP_EXTRACT_NONE;
152:   }

154:   /* fill sorting criterion context */
155:   switch (pep->which) {
156:     case PEP_LARGEST_MAGNITUDE:
157:       pep->sc->comparison    = SlepcCompareLargestMagnitude;
158:       pep->sc->comparisonctx = NULL;
159:       break;
160:     case PEP_SMALLEST_MAGNITUDE:
161:       pep->sc->comparison    = SlepcCompareSmallestMagnitude;
162:       pep->sc->comparisonctx = NULL;
163:       break;
164:     case PEP_LARGEST_REAL:
165:       pep->sc->comparison    = SlepcCompareLargestReal;
166:       pep->sc->comparisonctx = NULL;
167:       break;
168:     case PEP_SMALLEST_REAL:
169:       pep->sc->comparison    = SlepcCompareSmallestReal;
170:       pep->sc->comparisonctx = NULL;
171:       break;
172:     case PEP_LARGEST_IMAGINARY:
173:       pep->sc->comparison    = SlepcCompareLargestImaginary;
174:       pep->sc->comparisonctx = NULL;
175:       break;
176:     case PEP_SMALLEST_IMAGINARY:
177:       pep->sc->comparison    = SlepcCompareSmallestImaginary;
178:       pep->sc->comparisonctx = NULL;
179:       break;
180:     case PEP_TARGET_MAGNITUDE:
181:       pep->sc->comparison    = SlepcCompareTargetMagnitude;
182:       pep->sc->comparisonctx = &pep->target;
183:       break;
184:     case PEP_TARGET_REAL:
185:       pep->sc->comparison    = SlepcCompareTargetReal;
186:       pep->sc->comparisonctx = &pep->target;
187:       break;
188:     case PEP_TARGET_IMAGINARY:
189: #if defined(PETSC_USE_COMPLEX)
190:       pep->sc->comparison    = SlepcCompareTargetImaginary;
191:       pep->sc->comparisonctx = &pep->target;
192: #endif
193:       break;
194:     case PEP_ALL:
195:       pep->sc->comparison    = SlepcCompareSmallestReal;
196:       pep->sc->comparisonctx = NULL;
197:       break;
198:     case PEP_WHICH_USER:
199:       break;
200:   }
201:   pep->sc->map    = NULL;
202:   pep->sc->mapobj = NULL;

204:   /* fill sorting criterion for DS */
205:   if (pep->which!=PEP_ALL) {
206:     DSGetSlepcSC(pep->ds,&sc);
207:     RGIsTrivial(pep->rg,&istrivial);
208:     sc->rg            = istrivial? NULL: pep->rg;
209:     sc->comparison    = pep->sc->comparison;
210:     sc->comparisonctx = pep->sc->comparisonctx;
211:     sc->map           = SlepcMap_ST;
212:     sc->mapobj        = (PetscObject)pep->st;
213:   }
214:   /* setup ST */
215:   STSetUp(pep->st);

217:   /* compute matrix coefficients */
218:   STGetTransform(pep->st,&flg);
219:   if (!flg) {
220:     if (pep->solvematcoeffs) { STMatSetUp(pep->st,1.0,pep->solvematcoeffs); }
221:   } else {
222:     if (pep->basis!=PEP_BASIS_MONOMIAL) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_SUP,"Cannot use ST-transform with non-monomial basis in PEP");
223:   }

225:   /* compute scale factor if no set by user */
226:   PEPComputeScaleFactor(pep);

228:   /* build balancing matrix if required */
229:   if (pep->scale==PEP_SCALE_DIAGONAL || pep->scale==PEP_SCALE_BOTH) {
230:     if (!pep->Dl) {
231:       BVCreateVec(pep->V,&pep->Dl);
232:       PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->Dl);
233:     }
234:     if (!pep->Dr) {
235:       BVCreateVec(pep->V,&pep->Dr);
236:       PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->Dr);
237:     }
238:     PEPBuildDiagonalScaling(pep);
239:   }

241:   /* process initial vectors */
242:   if (pep->nini<0) {
243:     k = -pep->nini;
244:     if (k>pep->ncv) SETERRQ(PetscObjectComm((PetscObject)pep),1,"The number of initial vectors is larger than ncv");
245:     BVInsertVecs(pep->V,0,&k,pep->IS,PETSC_TRUE);
246:     SlepcBasisDestroy_Private(&pep->nini,&pep->IS);
247:     pep->nini = k;
248:   }
249:   PetscLogEventEnd(PEP_SetUp,pep,0,0,0);
250:   pep->state = PEP_STATE_SETUP;
251:   return(0);
252: }

254: /*@
255:    PEPSetOperators - Sets the coefficient matrices associated with the polynomial
256:    eigenvalue problem.

258:    Collective on PEP and Mat

260:    Input Parameters:
261: +  pep  - the eigenproblem solver context
262: .  nmat - number of matrices in array A
263: -  A    - the array of matrices associated with the eigenproblem

265:    Notes:
266:    The polynomial eigenproblem is defined as P(l)*x=0, where l is
267:    the eigenvalue, x is the eigenvector, and P(l) is defined as
268:    P(l) = A_0 + l*A_1 + ... + l^d*A_d, with d=nmat-1 (the degree of P).
269:    For non-monomial bases, this expression is different.

271:    Level: beginner

273: .seealso: PEPSolve(), PEPGetOperators(), PEPGetNumMatrices(), PEPSetBasis()
274: @*/
275: PetscErrorCode PEPSetOperators(PEP pep,PetscInt nmat,Mat A[])
276: {
278:   PetscInt       i,n,m,m0=0;

283:   if (nmat <= 0) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Non-positive value of nmat: %D",nmat);
284:   if (nmat <= 2) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Cannot solve linear eigenproblems with PEP; use EPS instead");

288:   MatGetSize(A[0],&m,&n);
289:   if (pep->state && n!=pep->n) { PEPReset(pep); }
290:   else if (pep->nmat) {
291:     MatDestroyMatrices(pep->nmat,&pep->A);
292:     PetscFree2(pep->pbc,pep->nrma);
293:     PetscFree(pep->solvematcoeffs);
294:   }

296:   PetscMalloc1(nmat,&pep->A);
297:   PetscCalloc2(3*nmat,&pep->pbc,nmat,&pep->nrma);
298:   for (i=0;i<nmat;i++) pep->pbc[i] = 1.0;  /* default to monomial basis */
299:   PetscLogObjectMemory((PetscObject)pep,nmat*sizeof(Mat)+4*nmat*sizeof(PetscReal));
300:   for (i=0;i<nmat;i++) {
303:     MatGetSize(A[i],&m,&n);
304:     if (m!=n) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"A[%D] is a non-square matrix",i);
305:     if (!i) m0 = m;
306:     if (m!=m0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_INCOMP,"Dimensions of matrices do not match with each other");
307:     PetscObjectReference((PetscObject)A[i]);
308:     pep->A[i] = A[i];
309:   }
310:   pep->nmat = nmat;
311:   pep->state = PEP_STATE_INITIAL;
312:   return(0);
313: }

315: /*@
316:    PEPGetOperators - Gets the matrices associated with the polynomial eigensystem.

318:    Not collective, though parallel Mats are returned if the PEP is parallel

320:    Input Parameters:
321: +  pep - the PEP context
322: -  k   - the index of the requested matrix (starting in 0)

324:    Output Parameter:
325: .  A - the requested matrix

327:    Level: intermediate

329: .seealso: PEPSolve(), PEPSetOperators(), PEPGetNumMatrices()
330: @*/
331: PetscErrorCode PEPGetOperators(PEP pep,PetscInt k,Mat *A)
332: {
336:   if (k<0 || k>=pep->nmat) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",pep->nmat-1);
337:   *A = pep->A[k];
338:   return(0);
339: }

341: /*@
342:    PEPGetNumMatrices - Returns the number of matrices stored in the PEP.

344:    Not collective

346:    Input Parameter:
347: .  pep - the PEP context

349:    Output Parameters:
350: .  nmat - the number of matrices passed in PEPSetOperators()

352:    Level: intermediate

354: .seealso: PEPSetOperators()
355: @*/
356: PetscErrorCode PEPGetNumMatrices(PEP pep,PetscInt *nmat)
357: {
361:   *nmat = pep->nmat;
362:   return(0);
363: }

365: /*@C
366:    PEPSetInitialSpace - Specify a basis of vectors that constitute the initial
367:    space, that is, the subspace from which the solver starts to iterate.

369:    Collective on PEP and Vec

371:    Input Parameter:
372: +  pep   - the polynomial eigensolver context
373: .  n     - number of vectors
374: -  is    - set of basis vectors of the initial space

376:    Notes:
377:    Some solvers start to iterate on a single vector (initial vector). In that case,
378:    the other vectors are ignored.

380:    These vectors do not persist from one PEPSolve() call to the other, so the
381:    initial space should be set every time.

383:    The vectors do not need to be mutually orthonormal, since they are explicitly
384:    orthonormalized internally.

386:    Common usage of this function is when the user can provide a rough approximation
387:    of the wanted eigenspace. Then, convergence may be faster.

389:    Level: intermediate
390: @*/
391: PetscErrorCode PEPSetInitialSpace(PEP pep,PetscInt n,Vec *is)
392: {

398:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
399:   SlepcBasisReference_Private(n,is,&pep->nini,&pep->IS);
400:   if (n>0) pep->state = PEP_STATE_INITIAL;
401:   return(0);
402: }

404: /*
405:   PEPSetDimensions_Default - Set reasonable values for ncv, mpd if not set
406:   by the user. This is called at setup.
407:  */
408: PetscErrorCode PEPSetDimensions_Default(PEP pep,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
409: {
411:   PetscBool      krylov;
412:   PetscInt       dim;

415:   PetscObjectTypeCompareAny((PetscObject)pep,&krylov,PEPTOAR,PEPQARNOLDI,"");
416:   dim = krylov?(pep->nmat-1)*pep->n:pep->n;
417:   if (*ncv) { /* ncv set */
418:     if (krylov) {
419:       if (*ncv<nev+1 && !(*ncv==nev && *ncv==dim)) SETERRQ(PetscObjectComm((PetscObject)pep),1,"The value of ncv must be at least nev+1");
420:     } else {
421:       if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)pep),1,"The value of ncv must be at least nev");
422:     }
423:   } else if (*mpd) { /* mpd set */
424:     *ncv = PetscMin(dim,nev+(*mpd));
425:   } else { /* neither set: defaults depend on nev being small or large */
426:     if (nev<500) *ncv = PetscMin(dim,PetscMax(2*nev,nev+15));
427:     else {
428:       *mpd = 500;
429:       *ncv = PetscMin(dim,nev+(*mpd));
430:     }
431:   }
432:   if (!*mpd) *mpd = *ncv;
433:   return(0);
434: }

436: /*@
437:    PEPAllocateSolution - Allocate memory storage for common variables such
438:    as eigenvalues and eigenvectors.

440:    Collective on PEP

442:    Input Parameters:
443: +  pep   - eigensolver context
444: -  extra - number of additional positions, used for methods that require a
445:            working basis slightly larger than ncv

447:    Developers Note:
448:    This is PETSC_EXTERN because it may be required by user plugin PEP
449:    implementations.

451:    Level: developer
452: @*/
453: PetscErrorCode PEPAllocateSolution(PEP pep,PetscInt extra)
454: {
456:   PetscInt       oldsize,newc,requested,requestedbv;
457:   PetscLogDouble cnt;
458:   Vec            t;

461:   requested = (pep->lineariz? pep->ncv: pep->ncv*(pep->nmat-1)) + extra;
462:   requestedbv = pep->ncv + extra;

464:   /* oldsize is zero if this is the first time setup is called */
465:   BVGetSizes(pep->V,NULL,NULL,&oldsize);

467:   /* allocate space for eigenvalues and friends */
468:   if (requested != oldsize || !pep->eigr) {
469:     PetscFree4(pep->eigr,pep->eigi,pep->errest,pep->perm);
470:     PetscMalloc4(requested,&pep->eigr,requested,&pep->eigi,requested,&pep->errest,requested,&pep->perm);
471:     newc = PetscMax(0,requested-oldsize);
472:     cnt = 2*newc*sizeof(PetscScalar) + newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
473:     PetscLogObjectMemory((PetscObject)pep,cnt);
474:   }

476:   /* allocate V */
477:   if (!pep->V) { PEPGetBV(pep,&pep->V); }
478:   if (!oldsize) {
479:     if (!((PetscObject)(pep->V))->type_name) {
480:       BVSetType(pep->V,BVSVEC);
481:     }
482:     STMatCreateVecsEmpty(pep->st,&t,NULL);
483:     BVSetSizesFromVec(pep->V,t,requestedbv);
484:     VecDestroy(&t);
485:   } else {
486:     BVResize(pep->V,requestedbv,PETSC_FALSE);
487:   }
488:   return(0);
489: }