Actual source code: epssetup.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:    EPS routines related to problem setup
 12: */

 14: #include <slepc/private/epsimpl.h>       /*I "slepceps.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 EPSSetFromOptions (before STSetFromOptions)
 20:    and also at EPSSetUp (in case EPSSetFromOptions was not called).
 21: */
 22: PetscErrorCode EPSSetDefaultST(EPS eps)
 23: {

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

 34: /*
 35:    This is done by preconditioned eigensolvers that use the PC only.
 36:    It sets STPRECOND with KSPPREONLY.
 37: */
 38: PetscErrorCode EPSSetDefaultST_Precond(EPS eps)
 39: {
 41:   KSP            ksp;

 44:   if (!((PetscObject)eps->st)->type_name) {
 45:     STSetType(eps->st,STPRECOND);
 46:     STPrecondSetKSPHasMat(eps->st,PETSC_TRUE);
 47:   }
 48:   STGetKSP(eps->st,&ksp);
 49:   if (!((PetscObject)ksp)->type_name) {
 50:     KSPSetType(ksp,KSPPREONLY);
 51:   }
 52:   return(0);
 53: }

 55: /*
 56:    This is done by preconditioned eigensolvers that can also use the KSP.
 57:    It sets STPRECOND with the default KSP (GMRES) and maxit=5.
 58: */
 59: PetscErrorCode EPSSetDefaultST_GMRES(EPS eps)
 60: {
 62:   KSP            ksp;

 65:   if (!((PetscObject)eps->st)->type_name) {
 66:     STSetType(eps->st,STPRECOND);
 67:     STPrecondSetKSPHasMat(eps->st,PETSC_TRUE);
 68:   }
 69:   STGetKSP(eps->st,&ksp);
 70:   if (!((PetscObject)ksp)->type_name) {
 71:     KSPSetType(ksp,KSPGMRES);
 72:     KSPSetTolerances(ksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,5);
 73:   }
 74:   return(0);
 75: }

 77: /*@
 78:    EPSSetUp - Sets up all the internal data structures necessary for the
 79:    execution of the eigensolver. Then calls STSetUp() for any set-up
 80:    operations associated to the ST object.

 82:    Collective on EPS

 84:    Input Parameter:
 85: .  eps   - eigenproblem solver context

 87:    Notes:
 88:    This function need not be called explicitly in most cases, since EPSSolve()
 89:    calls it. It can be useful when one wants to measure the set-up time
 90:    separately from the solve time.

 92:    Level: developer

 94: .seealso: EPSCreate(), EPSSolve(), EPSDestroy(), STSetUp(), EPSSetInitialSpace()
 95: @*/
 96: PetscErrorCode EPSSetUp(EPS eps)
 97: {
 99:   Mat            A,B;
100:   SlepcSC        sc;
101:   PetscInt       k,nmat;
102:   PetscBool      flg,istrivial,precond;
103: #if defined(PETSC_USE_COMPLEX)
104:   PetscScalar    sigma;
105: #endif

109:   if (eps->state) return(0);
110:   PetscLogEventBegin(EPS_SetUp,eps,0,0,0);

112:   /* reset the convergence flag from the previous solves */
113:   eps->reason = EPS_CONVERGED_ITERATING;

115:   /* Set default solver type (EPSSetFromOptions was not called) */
116:   if (!((PetscObject)eps)->type_name) {
117:     EPSSetType(eps,EPSKRYLOVSCHUR);
118:   }
119:   if (!eps->st) { EPSGetST(eps,&eps->st); }
120:   EPSSetDefaultST(eps);

122:   PetscObjectTypeCompare((PetscObject)eps->st,STPRECOND,&precond);
123:   if (eps->categ==EPS_CATEGORY_PRECOND && !precond) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver requires ST=PRECOND");
124:   if (eps->categ!=EPS_CATEGORY_PRECOND && precond) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"STPRECOND is intended for preconditioned eigensolvers only");

126:   STSetTransform(eps->st,PETSC_TRUE);
127:   if (eps->useds && !eps->ds) { EPSGetDS(eps,&eps->ds); }
128:   if (eps->twosided) {
129:     if (!eps->hasts) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"This solver does not support computing left eigenvectors (no two-sided variant)");
130:     if (eps->ishermitian) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Two-sided methods are not intended for symmetric problems");
131:     if (!eps->dsts) { DSDuplicate(eps->ds,&eps->dsts); }
132:   }
133:   if (!eps->rg) { EPSGetRG(eps,&eps->rg); }
134:   if (!((PetscObject)eps->rg)->type_name) {
135:     RGSetType(eps->rg,RGINTERVAL);
136:   }

138:   /* Set problem dimensions */
139:   STGetNumMatrices(eps->st,&nmat);
140:   if (!nmat) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"EPSSetOperators must be called first");
141:   STMatGetSize(eps->st,&eps->n,NULL);
142:   STMatGetLocalSize(eps->st,&eps->nloc,NULL);

144:   /* Set default problem type */
145:   if (!eps->problem_type) {
146:     if (nmat==1) {
147:       EPSSetProblemType(eps,EPS_NHEP);
148:     } else {
149:       EPSSetProblemType(eps,EPS_GNHEP);
150:     }
151:   } else if (nmat==1 && eps->isgeneralized) {
152:     PetscInfo(eps,"Eigenproblem set as generalized but no matrix B was provided; reverting to a standard eigenproblem\n");
153:     eps->isgeneralized = PETSC_FALSE;
154:     eps->problem_type = eps->ishermitian? EPS_HEP: EPS_NHEP;
155:   } else if (nmat>1 && !eps->isgeneralized) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Inconsistent EPS state");

157:   if (eps->nev > eps->n) eps->nev = eps->n;
158:   if (eps->ncv > eps->n) eps->ncv = eps->n;

160:   /* initialization of matrix norms */
161:   if (eps->conv==EPS_CONV_NORM) {
162:     if (!eps->nrma) {
163:       STGetMatrix(eps->st,0,&A);
164:       MatNorm(A,NORM_INFINITY,&eps->nrma);
165:     }
166:     if (nmat>1 && !eps->nrmb) {
167:       STGetMatrix(eps->st,1,&B);
168:       MatNorm(B,NORM_INFINITY,&eps->nrmb);
169:     }
170:   }

172:   /* call specific solver setup */
173:   (*eps->ops->setup)(eps);

175:   /* if purification is set, check that it really makes sense */
176:   if (eps->purify) {
177:     if (eps->categ==EPS_CATEGORY_PRECOND || eps->categ==EPS_CATEGORY_CONTOUR) eps->purify = PETSC_FALSE;
178:     else {
179:       if (!eps->isgeneralized) eps->purify = PETSC_FALSE;
180:       else if (!eps->ishermitian && !eps->ispositive) eps->purify = PETSC_FALSE;
181:       else {
182:         PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg);
183:         if (flg) eps->purify = PETSC_FALSE;
184:       }
185:     }
186:   }

188:   /* check extraction */
189:   PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STPRECOND,STSHIFT,"");
190:   if (!flg && eps->extraction && eps->extraction!=EPS_RITZ) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cannot use a spectral transformation combined with harmonic extraction");

192:   /* set tolerance if not yet set */
193:   if (eps->tol==PETSC_DEFAULT) eps->tol = SLEPC_DEFAULT_TOL;

195:   /* fill sorting criterion context */
196:   switch (eps->which) {
197:     case EPS_LARGEST_MAGNITUDE:
198:       eps->sc->comparison    = SlepcCompareLargestMagnitude;
199:       eps->sc->comparisonctx = NULL;
200:       break;
201:     case EPS_SMALLEST_MAGNITUDE:
202:       eps->sc->comparison    = SlepcCompareSmallestMagnitude;
203:       eps->sc->comparisonctx = NULL;
204:       break;
205:     case EPS_LARGEST_REAL:
206:       eps->sc->comparison    = SlepcCompareLargestReal;
207:       eps->sc->comparisonctx = NULL;
208:       break;
209:     case EPS_SMALLEST_REAL:
210:       eps->sc->comparison    = SlepcCompareSmallestReal;
211:       eps->sc->comparisonctx = NULL;
212:       break;
213:     case EPS_LARGEST_IMAGINARY:
214:       eps->sc->comparison    = SlepcCompareLargestImaginary;
215:       eps->sc->comparisonctx = NULL;
216:       break;
217:     case EPS_SMALLEST_IMAGINARY:
218:       eps->sc->comparison    = SlepcCompareSmallestImaginary;
219:       eps->sc->comparisonctx = NULL;
220:       break;
221:     case EPS_TARGET_MAGNITUDE:
222:       eps->sc->comparison    = SlepcCompareTargetMagnitude;
223:       eps->sc->comparisonctx = &eps->target;
224:       break;
225:     case EPS_TARGET_REAL:
226:       eps->sc->comparison    = SlepcCompareTargetReal;
227:       eps->sc->comparisonctx = &eps->target;
228:       break;
229:     case EPS_TARGET_IMAGINARY:
230: #if defined(PETSC_USE_COMPLEX)
231:       eps->sc->comparison    = SlepcCompareTargetImaginary;
232:       eps->sc->comparisonctx = &eps->target;
233: #endif
234:       break;
235:     case EPS_ALL:
236:       eps->sc->comparison    = SlepcCompareSmallestReal;
237:       eps->sc->comparisonctx = NULL;
238:       break;
239:     case EPS_WHICH_USER:
240:       break;
241:   }
242:   eps->sc->map    = NULL;
243:   eps->sc->mapobj = NULL;

245:   /* fill sorting criterion for DS */
246:   if (eps->useds && eps->which!=EPS_ALL) {
247:     DSGetSlepcSC(eps->ds,&sc);
248:     RGIsTrivial(eps->rg,&istrivial);
249:     sc->rg            = istrivial? NULL: eps->rg;
250:     sc->comparison    = eps->sc->comparison;
251:     sc->comparisonctx = eps->sc->comparisonctx;
252:     sc->map           = SlepcMap_ST;
253:     sc->mapobj        = (PetscObject)eps->st;
254:     if (eps->twosided) {
255:       DSSetSlepcSC(eps->dsts,sc);
256:     }
257:   }

259:   /* Build balancing matrix if required */
260:   if (eps->balance!=EPS_BALANCE_USER) {
261:     STSetBalanceMatrix(eps->st,NULL);
262:     if (!eps->ishermitian && (eps->balance==EPS_BALANCE_ONESIDE || eps->balance==EPS_BALANCE_TWOSIDE)) {
263:       if (!eps->D) {
264:         BVCreateVec(eps->V,&eps->D);
265:         PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->D);
266:       }
267:       EPSBuildBalance_Krylov(eps);
268:       STSetBalanceMatrix(eps->st,eps->D);
269:     }
270:   }

272:   /* Setup ST */
273:   STSetUp(eps->st);

275: #if defined(PETSC_USE_COMPLEX)
276:   STGetShift(eps->st,&sigma);
277:   if (eps->ishermitian && PetscImaginaryPart(sigma) != 0.0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Hermitian problems are not compatible with complex shifts");
278: #endif
279:   PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&flg);
280:   if (flg && eps->problem_type == EPS_PGNHEP) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_SUP,"Cayley spectral transformation is not compatible with PGNHEP");
281:   PetscObjectTypeCompareAny((PetscObject)eps->st,&flg,STSINVERT,STCAYLEY,"");
282:   if (flg && (eps->categ==EPS_CATEGORY_KRYLOV || eps->categ==EPS_CATEGORY_OTHER) && (eps->which!=EPS_TARGET_MAGNITUDE && eps->which!=EPS_TARGET_REAL && eps->which!=EPS_TARGET_IMAGINARY && eps->which!=EPS_ALL && eps->which!=EPS_WHICH_USER)) SETERRQ(PetscObjectComm((PetscObject)eps),1,"Shift-and-invert requires a target 'which' (see EPSSetWhichEigenpairs), for instance -st_type sinvert -eps_target 0 -eps_target_magnitude");

284:   /* process deflation and initial vectors */
285:   if (eps->nds<0) {
286:     k = -eps->nds;
287:     BVInsertConstraints(eps->V,&k,eps->defl);
288:     SlepcBasisDestroy_Private(&eps->nds,&eps->defl);
289:     eps->nds = k;
290:     STCheckNullSpace(eps->st,eps->V);
291:   }
292:   if (eps->nini<0) {
293:     k = -eps->nini;
294:     if (k>eps->ncv) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The number of initial vectors is larger than ncv");
295:     BVInsertVecs(eps->V,0,&k,eps->IS,PETSC_TRUE);
296:     SlepcBasisDestroy_Private(&eps->nini,&eps->IS);
297:     eps->nini = k;
298:   }

300:   PetscLogEventEnd(EPS_SetUp,eps,0,0,0);
301:   eps->state = EPS_STATE_SETUP;
302:   return(0);
303: }

305: /*@
306:    EPSSetOperators - Sets the matrices associated with the eigenvalue problem.

308:    Collective on EPS and Mat

310:    Input Parameters:
311: +  eps - the eigenproblem solver context
312: .  A  - the matrix associated with the eigensystem
313: -  B  - the second matrix in the case of generalized eigenproblems

315:    Notes:
316:    To specify a standard eigenproblem, use NULL for parameter B.

318:    It must be called before EPSSetUp(). If it is called again after EPSSetUp() then
319:    the EPS object is reset.

321:    Level: beginner

323: .seealso: EPSSolve(), EPSSetUp(), EPSReset(), EPSGetST(), STGetMatrix()
324: @*/
325: PetscErrorCode EPSSetOperators(EPS eps,Mat A,Mat B)
326: {
328:   PetscInt       m,n,m0,nmat;
329:   Mat            mat[2];


338:   /* Check for square matrices */
339:   MatGetSize(A,&m,&n);
340:   if (m!=n) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"A is a non-square matrix");
341:   if (B) {
342:     MatGetSize(B,&m0,&n);
343:     if (m0!=n) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"B is a non-square matrix");
344:     if (m!=m0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_INCOMP,"Dimensions of A and B do not match");
345:   }
346:   if (eps->state && n!=eps->n) { EPSReset(eps); }
347:   eps->nrma = 0.0;
348:   eps->nrmb = 0.0;
349:   if (!eps->st) { EPSGetST(eps,&eps->st); }
350:   mat[0] = A;
351:   if (B) {
352:     mat[1] = B;
353:     nmat = 2;
354:   } else nmat = 1;
355:   STSetMatrices(eps->st,nmat,mat);
356:   eps->state = EPS_STATE_INITIAL;
357:   return(0);
358: }

360: /*@
361:    EPSGetOperators - Gets the matrices associated with the eigensystem.

363:    Collective on EPS and Mat

365:    Input Parameter:
366: .  eps - the EPS context

368:    Output Parameters:
369: +  A  - the matrix associated with the eigensystem
370: -  B  - the second matrix in the case of generalized eigenproblems

372:    Level: intermediate

374: .seealso: EPSSolve(), EPSGetST(), STGetMatrix(), STSetMatrices()
375: @*/
376: PetscErrorCode EPSGetOperators(EPS eps,Mat *A,Mat *B)
377: {
379:   ST             st;
380:   PetscInt       k;

384:   EPSGetST(eps,&st);
385:   if (A) { STGetMatrix(st,0,A); }
386:   if (B) {
387:     STGetNumMatrices(st,&k);
388:     if (k==1) B = NULL;
389:     else {
390:       STGetMatrix(st,1,B);
391:     }
392:   }
393:   return(0);
394: }

396: /*@
397:    EPSSetDeflationSpace - Specify a basis of vectors that constitute the deflation
398:    space.

400:    Collective on EPS and Vec

402:    Input Parameter:
403: +  eps - the eigenproblem solver context
404: .  n   - number of vectors
405: -  v   - set of basis vectors of the deflation space

407:    Notes:
408:    When a deflation space is given, the eigensolver seeks the eigensolution
409:    in the restriction of the problem to the orthogonal complement of this
410:    space. This can be used for instance in the case that an invariant
411:    subspace is known beforehand (such as the nullspace of the matrix).

413:    These vectors do not persist from one EPSSolve() call to the other, so the
414:    deflation space should be set every time.

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

419:    Level: intermediate

421: .seealso: EPSSetInitialSpace()
422: @*/
423: PetscErrorCode EPSSetDeflationSpace(EPS eps,PetscInt n,Vec *v)
424: {

430:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n out of range");
431:   SlepcBasisReference_Private(n,v,&eps->nds,&eps->defl);
432:   if (n>0) eps->state = EPS_STATE_INITIAL;
433:   return(0);
434: }

436: /*@C
437:    EPSSetInitialSpace - Specify a basis of vectors that constitute the initial
438:    space, that is, the subspace from which the solver starts to iterate.

440:    Collective on EPS and Vec

442:    Input Parameter:
443: +  eps - the eigenproblem solver context
444: .  n   - number of vectors
445: -  is  - set of basis vectors of the initial space

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

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

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

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

460:    Level: intermediate

462: .seealso: EPSSetDeflationSpace()
463: @*/
464: PetscErrorCode EPSSetInitialSpace(EPS eps,PetscInt n,Vec *is)
465: {

471:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"Argument n cannot be negative");
472:   SlepcBasisReference_Private(n,is,&eps->nini,&eps->IS);
473:   if (n>0) eps->state = EPS_STATE_INITIAL;
474:   return(0);
475: }

477: /*
478:   EPSSetDimensions_Default - Set reasonable values for ncv, mpd if not set
479:   by the user. This is called at setup.
480:  */
481: PetscErrorCode EPSSetDimensions_Default(EPS eps,PetscInt nev,PetscInt *ncv,PetscInt *mpd)
482: {
484:   PetscBool      krylov;

487:   if (*ncv) { /* ncv set */
488:     PetscObjectTypeCompareAny((PetscObject)eps,&krylov,EPSKRYLOVSCHUR,EPSARNOLDI,EPSLANCZOS,"");
489:     if (krylov) {
490:       if (*ncv<nev+1 && !(*ncv==nev && *ncv==eps->n)) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev+1");
491:     } else {
492:       if (*ncv<nev) SETERRQ(PetscObjectComm((PetscObject)eps),1,"The value of ncv must be at least nev");
493:     }
494:   } else if (*mpd) { /* mpd set */
495:     *ncv = PetscMin(eps->n,nev+(*mpd));
496:   } else { /* neither set: defaults depend on nev being small or large */
497:     if (nev<500) *ncv = PetscMin(eps->n,PetscMax(2*nev,nev+15));
498:     else {
499:       *mpd = 500;
500:       *ncv = PetscMin(eps->n,nev+(*mpd));
501:     }
502:   }
503:   if (!*mpd) *mpd = *ncv;
504:   return(0);
505: }

507: /*@
508:    EPSAllocateSolution - Allocate memory storage for common variables such
509:    as eigenvalues and eigenvectors.

511:    Collective on EPS

513:    Input Parameters:
514: +  eps   - eigensolver context
515: -  extra - number of additional positions, used for methods that require a
516:            working basis slightly larger than ncv

518:    Developers Note:
519:    This is SLEPC_EXTERN because it may be required by user plugin EPS
520:    implementations.

522:    Level: developer
523: @*/
524: PetscErrorCode EPSAllocateSolution(EPS eps,PetscInt extra)
525: {
527:   PetscInt       oldsize,newc,requested;
528:   PetscLogDouble cnt;
529:   Vec            t;

532:   requested = eps->ncv + extra;

534:   /* oldsize is zero if this is the first time setup is called */
535:   BVGetSizes(eps->V,NULL,NULL,&oldsize);
536:   newc = PetscMax(0,requested-oldsize);

538:   /* allocate space for eigenvalues and friends */
539:   if (requested != oldsize || !eps->eigr) {
540:     PetscFree4(eps->eigr,eps->eigi,eps->errest,eps->perm);
541:     PetscMalloc4(requested,&eps->eigr,requested,&eps->eigi,requested,&eps->errest,requested,&eps->perm);
542:     cnt = 2*newc*sizeof(PetscScalar) + 2*newc*sizeof(PetscReal) + newc*sizeof(PetscInt);
543:     PetscLogObjectMemory((PetscObject)eps,cnt);
544:   }

546:   /* workspace for the case of arbitrary selection */
547:   if (eps->arbitrary) {
548:     if (eps->rr) {
549:       PetscFree2(eps->rr,eps->ri);
550:     }
551:     PetscMalloc2(requested,&eps->rr,requested,&eps->ri);
552:     PetscLogObjectMemory((PetscObject)eps,2*newc*sizeof(PetscScalar));
553:   }

555:   /* allocate V */
556:   if (!eps->V) { EPSGetBV(eps,&eps->V); }
557:   if (!oldsize) {
558:     if (!((PetscObject)(eps->V))->type_name) {
559:       BVSetType(eps->V,BVSVEC);
560:     }
561:     STMatCreateVecsEmpty(eps->st,&t,NULL);
562:     BVSetSizesFromVec(eps->V,t,requested);
563:     VecDestroy(&t);
564:   } else {
565:     BVResize(eps->V,requested,PETSC_FALSE);
566:   }

568:   /* allocate W */
569:   if (eps->twosided) {
570:     BVDestroy(&eps->W);
571:     BVDuplicate(eps->V,&eps->W);
572:   }
573:   return(0);
574: }