Actual source code: epsbasic.c

slepc-3.13.3 2020-06-14
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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:    Basic EPS routines
 12: */

 14:  #include <slepc/private/epsimpl.h>

 16: PetscFunctionList EPSList = 0;
 17: PetscBool         EPSRegisterAllCalled = PETSC_FALSE;
 18: PetscClassId      EPS_CLASSID = 0;
 19: PetscLogEvent     EPS_SetUp = 0,EPS_Solve = 0;

 21: /*@
 22:    EPSCreate - Creates the default EPS context.

 24:    Collective

 26:    Input Parameter:
 27: .  comm - MPI communicator

 29:    Output Parameter:
 30: .  eps - location to put the EPS context

 32:    Note:
 33:    The default EPS type is EPSKRYLOVSCHUR

 35:    Level: beginner

 37: .seealso: EPSSetUp(), EPSSolve(), EPSDestroy(), EPS
 38: @*/
 39: PetscErrorCode EPSCreate(MPI_Comm comm,EPS *outeps)
 40: {
 42:   EPS            eps;

 46:   *outeps = 0;
 47:   EPSInitializePackage();
 48:   SlepcHeaderCreate(eps,EPS_CLASSID,"EPS","Eigenvalue Problem Solver","EPS",comm,EPSDestroy,EPSView);

 50:   eps->max_it          = PETSC_DEFAULT;
 51:   eps->nev             = 1;
 52:   eps->ncv             = PETSC_DEFAULT;
 53:   eps->mpd             = PETSC_DEFAULT;
 54:   eps->nini            = 0;
 55:   eps->nds             = 0;
 56:   eps->target          = 0.0;
 57:   eps->tol             = PETSC_DEFAULT;
 58:   eps->conv            = EPS_CONV_REL;
 59:   eps->stop            = EPS_STOP_BASIC;
 60:   eps->which           = (EPSWhich)0;
 61:   eps->inta            = 0.0;
 62:   eps->intb            = 0.0;
 63:   eps->problem_type    = (EPSProblemType)0;
 64:   eps->extraction      = EPS_RITZ;
 65:   eps->balance         = EPS_BALANCE_NONE;
 66:   eps->balance_its     = 5;
 67:   eps->balance_cutoff  = 1e-8;
 68:   eps->trueres         = PETSC_FALSE;
 69:   eps->trackall        = PETSC_FALSE;
 70:   eps->purify          = PETSC_TRUE;
 71:   eps->twosided        = PETSC_FALSE;

 73:   eps->converged       = EPSConvergedRelative;
 74:   eps->convergeduser   = NULL;
 75:   eps->convergeddestroy= NULL;
 76:   eps->stopping        = EPSStoppingBasic;
 77:   eps->stoppinguser    = NULL;
 78:   eps->stoppingdestroy = NULL;
 79:   eps->arbitrary       = NULL;
 80:   eps->convergedctx    = NULL;
 81:   eps->stoppingctx     = NULL;
 82:   eps->arbitraryctx    = NULL;
 83:   eps->numbermonitors  = 0;

 85:   eps->st              = NULL;
 86:   eps->ds              = NULL;
 87:   eps->dsts            = NULL;
 88:   eps->V               = NULL;
 89:   eps->W               = NULL;
 90:   eps->rg              = NULL;
 91:   eps->D               = NULL;
 92:   eps->IS              = NULL;
 93:   eps->ISL             = NULL;
 94:   eps->defl            = NULL;
 95:   eps->eigr            = NULL;
 96:   eps->eigi            = NULL;
 97:   eps->errest          = NULL;
 98:   eps->rr              = NULL;
 99:   eps->ri              = NULL;
100:   eps->perm            = NULL;
101:   eps->nwork           = 0;
102:   eps->work            = NULL;
103:   eps->data            = NULL;

105:   eps->state           = EPS_STATE_INITIAL;
106:   eps->categ           = EPS_CATEGORY_KRYLOV;
107:   eps->nconv           = 0;
108:   eps->its             = 0;
109:   eps->nloc            = 0;
110:   eps->nrma            = 0.0;
111:   eps->nrmb            = 0.0;
112:   eps->useds           = PETSC_FALSE;
113:   eps->hasts           = PETSC_FALSE;
114:   eps->isgeneralized   = PETSC_FALSE;
115:   eps->ispositive      = PETSC_FALSE;
116:   eps->ishermitian     = PETSC_FALSE;
117:   eps->reason          = EPS_CONVERGED_ITERATING;

119:   PetscNewLog(eps,&eps->sc);
120:   *outeps = eps;
121:   return(0);
122: }

124: /*@C
125:    EPSSetType - Selects the particular solver to be used in the EPS object.

127:    Logically Collective on eps

129:    Input Parameters:
130: +  eps  - the eigensolver context
131: -  type - a known method

133:    Options Database Key:
134: .  -eps_type <method> - Sets the method; use -help for a list
135:     of available methods

137:    Notes:
138:    See "slepc/include/slepceps.h" for available methods. The default
139:    is EPSKRYLOVSCHUR.

141:    Normally, it is best to use the EPSSetFromOptions() command and
142:    then set the EPS type from the options database rather than by using
143:    this routine.  Using the options database provides the user with
144:    maximum flexibility in evaluating the different available methods.
145:    The EPSSetType() routine is provided for those situations where it
146:    is necessary to set the iterative solver independently of the command
147:    line or options database.

149:    Level: intermediate

151: .seealso: STSetType(), EPSType
152: @*/
153: PetscErrorCode EPSSetType(EPS eps,EPSType type)
154: {
155:   PetscErrorCode ierr,(*r)(EPS);
156:   PetscBool      match;


162:   PetscObjectTypeCompare((PetscObject)eps,type,&match);
163:   if (match) return(0);

165:   PetscFunctionListFind(EPSList,type,&r);
166:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown EPS type given: %s",type);

168:   if (eps->ops->destroy) { (*eps->ops->destroy)(eps); }
169:   PetscMemzero(eps->ops,sizeof(struct _EPSOps));

171:   eps->state = EPS_STATE_INITIAL;
172:   PetscObjectChangeTypeName((PetscObject)eps,type);
173:   (*r)(eps);
174:   return(0);
175: }

177: /*@C
178:    EPSGetType - Gets the EPS type as a string from the EPS object.

180:    Not Collective

182:    Input Parameter:
183: .  eps - the eigensolver context

185:    Output Parameter:
186: .  name - name of EPS method

188:    Level: intermediate

190: .seealso: EPSSetType()
191: @*/
192: PetscErrorCode EPSGetType(EPS eps,EPSType *type)
193: {
197:   *type = ((PetscObject)eps)->type_name;
198:   return(0);
199: }

201: /*@C
202:    EPSRegister - Adds a method to the eigenproblem solver package.

204:    Not Collective

206:    Input Parameters:
207: +  name - name of a new user-defined solver
208: -  function - routine to create the solver context

210:    Notes:
211:    EPSRegister() may be called multiple times to add several user-defined solvers.

213:    Sample usage:
214: .vb
215:     EPSRegister("my_solver",MySolverCreate);
216: .ve

218:    Then, your solver can be chosen with the procedural interface via
219: $     EPSSetType(eps,"my_solver")
220:    or at runtime via the option
221: $     -eps_type my_solver

223:    Level: advanced

225: .seealso: EPSRegisterAll()
226: @*/
227: PetscErrorCode EPSRegister(const char *name,PetscErrorCode (*function)(EPS))
228: {

232:   EPSInitializePackage();
233:   PetscFunctionListAdd(&EPSList,name,function);
234:   return(0);
235: }

237: /*@
238:    EPSReset - Resets the EPS context to the initial state (prior to setup)
239:    and destroys any allocated Vecs and Mats.

241:    Collective on eps

243:    Input Parameter:
244: .  eps - eigensolver context obtained from EPSCreate()

246:    Note:
247:    This can be used when a problem of different matrix size wants to be solved.
248:    All options that have previously been set are preserved, so in a next use
249:    the solver configuration is the same, but new sizes for matrices and vectors
250:    are allowed.

252:    Level: advanced

254: .seealso: EPSDestroy()
255: @*/
256: PetscErrorCode EPSReset(EPS eps)
257: {

262:   if (!eps) return(0);
263:   if (eps->ops->reset) { (eps->ops->reset)(eps); }
264:   if (eps->st) { STReset(eps->st); }
265:   VecDestroy(&eps->D);
266:   BVDestroy(&eps->V);
267:   BVDestroy(&eps->W);
268:   VecDestroyVecs(eps->nwork,&eps->work);
269:   eps->nwork = 0;
270:   eps->state = EPS_STATE_INITIAL;
271:   return(0);
272: }

274: /*@
275:    EPSDestroy - Destroys the EPS context.

277:    Collective on eps

279:    Input Parameter:
280: .  eps - eigensolver context obtained from EPSCreate()

282:    Level: beginner

284: .seealso: EPSCreate(), EPSSetUp(), EPSSolve()
285: @*/
286: PetscErrorCode EPSDestroy(EPS *eps)
287: {

291:   if (!*eps) return(0);
293:   if (--((PetscObject)(*eps))->refct > 0) { *eps = 0; return(0); }
294:   EPSReset(*eps);
295:   if ((*eps)->ops->destroy) { (*(*eps)->ops->destroy)(*eps); }
296:   if ((*eps)->eigr) {
297:     PetscFree4((*eps)->eigr,(*eps)->eigi,(*eps)->errest,(*eps)->perm);
298:   }
299:   if ((*eps)->rr) {
300:     PetscFree2((*eps)->rr,(*eps)->ri);
301:   }
302:   STDestroy(&(*eps)->st);
303:   RGDestroy(&(*eps)->rg);
304:   DSDestroy(&(*eps)->ds);
305:   DSDestroy(&(*eps)->dsts);
306:   PetscFree((*eps)->sc);
307:   /* just in case the initial vectors have not been used */
308:   SlepcBasisDestroy_Private(&(*eps)->nds,&(*eps)->defl);
309:   SlepcBasisDestroy_Private(&(*eps)->nini,&(*eps)->IS);
310:   SlepcBasisDestroy_Private(&(*eps)->ninil,&(*eps)->ISL);
311:   if ((*eps)->convergeddestroy) {
312:     (*(*eps)->convergeddestroy)((*eps)->convergedctx);
313:   }
314:   EPSMonitorCancel(*eps);
315:   PetscHeaderDestroy(eps);
316:   return(0);
317: }

319: /*@
320:    EPSSetTarget - Sets the value of the target.

322:    Logically Collective on eps

324:    Input Parameters:
325: +  eps    - eigensolver context
326: -  target - the value of the target

328:    Options Database Key:
329: .  -eps_target <scalar> - the value of the target

331:    Notes:
332:    The target is a scalar value used to determine the portion of the spectrum
333:    of interest. It is used in combination with EPSSetWhichEigenpairs().

335:    In the case of complex scalars, a complex value can be provided in the
336:    command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
337:    -eps_target 1.0+2.0i

339:    Level: intermediate

341: .seealso: EPSGetTarget(), EPSSetWhichEigenpairs()
342: @*/
343: PetscErrorCode EPSSetTarget(EPS eps,PetscScalar target)
344: {

350:   eps->target = target;
351:   if (!eps->st) { EPSGetST(eps,&eps->st); }
352:   STSetDefaultShift(eps->st,target);
353:   return(0);
354: }

356: /*@
357:    EPSGetTarget - Gets the value of the target.

359:    Not Collective

361:    Input Parameter:
362: .  eps - eigensolver context

364:    Output Parameter:
365: .  target - the value of the target

367:    Note:
368:    If the target was not set by the user, then zero is returned.

370:    Level: intermediate

372: .seealso: EPSSetTarget()
373: @*/
374: PetscErrorCode EPSGetTarget(EPS eps,PetscScalar* target)
375: {
379:   *target = eps->target;
380:   return(0);
381: }

383: /*@
384:    EPSSetInterval - Defines the computational interval for spectrum slicing.

386:    Logically Collective on eps

388:    Input Parameters:
389: +  eps  - eigensolver context
390: .  inta - left end of the interval
391: -  intb - right end of the interval

393:    Options Database Key:
394: .  -eps_interval <a,b> - set [a,b] as the interval of interest

396:    Notes:
397:    Spectrum slicing is a technique employed for computing all eigenvalues of
398:    symmetric eigenproblems in a given interval. This function provides the
399:    interval to be considered. It must be used in combination with EPS_ALL, see
400:    EPSSetWhichEigenpairs().

402:    In the command-line option, two values must be provided. For an open interval,
403:    one can give an infinite, e.g., -eps_interval 1.0,inf or -eps_interval -inf,1.0.
404:    An open interval in the programmatic interface can be specified with
405:    PETSC_MAX_REAL and -PETSC_MAX_REAL.

407:    Level: intermediate

409: .seealso: EPSGetInterval(), EPSSetWhichEigenpairs()
410: @*/
411: PetscErrorCode EPSSetInterval(EPS eps,PetscReal inta,PetscReal intb)
412: {
417:   if (inta>intb) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
418:   if (eps->inta != inta || eps->intb != intb) {
419:     eps->inta = inta;
420:     eps->intb = intb;
421:     eps->state = EPS_STATE_INITIAL;
422:   }
423:   return(0);
424: }

426: /*@
427:    EPSGetInterval - Gets the computational interval for spectrum slicing.

429:    Not Collective

431:    Input Parameter:
432: .  eps - eigensolver context

434:    Output Parameters:
435: +  inta - left end of the interval
436: -  intb - right end of the interval

438:    Level: intermediate

440:    Note:
441:    If the interval was not set by the user, then zeros are returned.

443: .seealso: EPSSetInterval()
444: @*/
445: PetscErrorCode EPSGetInterval(EPS eps,PetscReal* inta,PetscReal* intb)
446: {
449:   if (inta) *inta = eps->inta;
450:   if (intb) *intb = eps->intb;
451:   return(0);
452: }

454: /*@
455:    EPSSetST - Associates a spectral transformation object to the eigensolver.

457:    Collective on eps

459:    Input Parameters:
460: +  eps - eigensolver context obtained from EPSCreate()
461: -  st   - the spectral transformation object

463:    Note:
464:    Use EPSGetST() to retrieve the spectral transformation context (for example,
465:    to free it at the end of the computations).

467:    Level: advanced

469: .seealso: EPSGetST()
470: @*/
471: PetscErrorCode EPSSetST(EPS eps,ST st)
472: {

479:   PetscObjectReference((PetscObject)st);
480:   STDestroy(&eps->st);
481:   eps->st = st;
482:   PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->st);
483:   return(0);
484: }

486: /*@
487:    EPSGetST - Obtain the spectral transformation (ST) object associated
488:    to the eigensolver object.

490:    Not Collective

492:    Input Parameters:
493: .  eps - eigensolver context obtained from EPSCreate()

495:    Output Parameter:
496: .  st - spectral transformation context

498:    Level: intermediate

500: .seealso: EPSSetST()
501: @*/
502: PetscErrorCode EPSGetST(EPS eps,ST *st)
503: {

509:   if (!eps->st) {
510:     STCreate(PetscObjectComm((PetscObject)eps),&eps->st);
511:     PetscObjectIncrementTabLevel((PetscObject)eps->st,(PetscObject)eps,0);
512:     PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->st);
513:     PetscObjectSetOptions((PetscObject)eps->st,((PetscObject)eps)->options);
514:   }
515:   *st = eps->st;
516:   return(0);
517: }

519: /*@
520:    EPSSetBV - Associates a basis vectors object to the eigensolver.

522:    Collective on eps

524:    Input Parameters:
525: +  eps - eigensolver context obtained from EPSCreate()
526: -  V   - the basis vectors object

528:    Level: advanced

530: .seealso: EPSGetBV()
531: @*/
532: PetscErrorCode EPSSetBV(EPS eps,BV V)
533: {

540:   PetscObjectReference((PetscObject)V);
541:   BVDestroy(&eps->V);
542:   eps->V = V;
543:   PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->V);
544:   return(0);
545: }

547: /*@
548:    EPSGetBV - Obtain the basis vectors object associated to the eigensolver object.

550:    Not Collective

552:    Input Parameters:
553: .  eps - eigensolver context obtained from EPSCreate()

555:    Output Parameter:
556: .  V - basis vectors context

558:    Level: advanced

560: .seealso: EPSSetBV()
561: @*/
562: PetscErrorCode EPSGetBV(EPS eps,BV *V)
563: {

569:   if (!eps->V) {
570:     BVCreate(PetscObjectComm((PetscObject)eps),&eps->V);
571:     PetscObjectIncrementTabLevel((PetscObject)eps->V,(PetscObject)eps,0);
572:     PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->V);
573:     PetscObjectSetOptions((PetscObject)eps->V,((PetscObject)eps)->options);
574:   }
575:   *V = eps->V;
576:   return(0);
577: }

579: /*@
580:    EPSSetRG - Associates a region object to the eigensolver.

582:    Collective on eps

584:    Input Parameters:
585: +  eps - eigensolver context obtained from EPSCreate()
586: -  rg  - the region object

588:    Note:
589:    Use EPSGetRG() to retrieve the region context (for example,
590:    to free it at the end of the computations).

592:    Level: advanced

594: .seealso: EPSGetRG()
595: @*/
596: PetscErrorCode EPSSetRG(EPS eps,RG rg)
597: {

604:   PetscObjectReference((PetscObject)rg);
605:   RGDestroy(&eps->rg);
606:   eps->rg = rg;
607:   PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->rg);
608:   return(0);
609: }

611: /*@
612:    EPSGetRG - Obtain the region object associated to the eigensolver.

614:    Not Collective

616:    Input Parameters:
617: .  eps - eigensolver context obtained from EPSCreate()

619:    Output Parameter:
620: .  rg - region context

622:    Level: advanced

624: .seealso: EPSSetRG()
625: @*/
626: PetscErrorCode EPSGetRG(EPS eps,RG *rg)
627: {

633:   if (!eps->rg) {
634:     RGCreate(PetscObjectComm((PetscObject)eps),&eps->rg);
635:     PetscObjectIncrementTabLevel((PetscObject)eps->rg,(PetscObject)eps,0);
636:     PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->rg);
637:     PetscObjectSetOptions((PetscObject)eps->rg,((PetscObject)eps)->options);
638:   }
639:   *rg = eps->rg;
640:   return(0);
641: }

643: /*@
644:    EPSSetDS - Associates a direct solver object to the eigensolver.

646:    Collective on eps

648:    Input Parameters:
649: +  eps - eigensolver context obtained from EPSCreate()
650: -  ds  - the direct solver object

652:    Note:
653:    Use EPSGetDS() to retrieve the direct solver context (for example,
654:    to free it at the end of the computations).

656:    Level: advanced

658: .seealso: EPSGetDS()
659: @*/
660: PetscErrorCode EPSSetDS(EPS eps,DS ds)
661: {

668:   PetscObjectReference((PetscObject)ds);
669:   DSDestroy(&eps->ds);
670:   eps->ds = ds;
671:   PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->ds);
672:   return(0);
673: }

675: /*@
676:    EPSGetDS - Obtain the direct solver object associated to the eigensolver object.

678:    Not Collective

680:    Input Parameters:
681: .  eps - eigensolver context obtained from EPSCreate()

683:    Output Parameter:
684: .  ds - direct solver context

686:    Level: advanced

688: .seealso: EPSSetDS()
689: @*/
690: PetscErrorCode EPSGetDS(EPS eps,DS *ds)
691: {

697:   if (!eps->ds) {
698:     DSCreate(PetscObjectComm((PetscObject)eps),&eps->ds);
699:     PetscObjectIncrementTabLevel((PetscObject)eps->ds,(PetscObject)eps,0);
700:     PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->ds);
701:     PetscObjectSetOptions((PetscObject)eps->ds,((PetscObject)eps)->options);
702:   }
703:   *ds = eps->ds;
704:   return(0);
705: }

707: /*@
708:    EPSIsGeneralized - Ask if the EPS object corresponds to a generalized
709:    eigenvalue problem.

711:    Not collective

713:    Input Parameter:
714: .  eps - the eigenproblem solver context

716:    Output Parameter:
717: .  is - the answer

719:    Level: intermediate

721: .seealso: EPSIsHermitian(), EPSIsPositive()
722: @*/
723: PetscErrorCode EPSIsGeneralized(EPS eps,PetscBool* is)
724: {
728:   *is = eps->isgeneralized;
729:   return(0);
730: }

732: /*@
733:    EPSIsHermitian - Ask if the EPS object corresponds to a Hermitian
734:    eigenvalue problem.

736:    Not collective

738:    Input Parameter:
739: .  eps - the eigenproblem solver context

741:    Output Parameter:
742: .  is - the answer

744:    Level: intermediate

746: .seealso: EPSIsGeneralized(), EPSIsPositive()
747: @*/
748: PetscErrorCode EPSIsHermitian(EPS eps,PetscBool* is)
749: {
753:   *is = eps->ishermitian;
754:   return(0);
755: }

757: /*@
758:    EPSIsPositive - Ask if the EPS object corresponds to an eigenvalue
759:    problem type that requires a positive (semi-) definite matrix B.

761:    Not collective

763:    Input Parameter:
764: .  eps - the eigenproblem solver context

766:    Output Parameter:
767: .  is - the answer

769:    Level: intermediate

771: .seealso: EPSIsGeneralized(), EPSIsHermitian()
772: @*/
773: PetscErrorCode EPSIsPositive(EPS eps,PetscBool* is)
774: {
778:   *is = eps->ispositive;
779:   return(0);
780: }