Actual source code: pepbasic.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:    Basic PEP routines
 12: */

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

 16: PetscFunctionList PEPList = 0;
 17: PetscBool         PEPRegisterAllCalled = PETSC_FALSE;
 18: PetscClassId      PEP_CLASSID = 0;
 19: PetscLogEvent     PEP_SetUp = 0,PEP_Solve = 0,PEP_Refine = 0;

 21: /*@
 22:    PEPCreate - Creates the default PEP context.

 24:    Collective on MPI_Comm

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

 29:    Output Parameter:
 30: .  pep - location to put the PEP context

 32:    Note:
 33:    The default PEP type is PEPTOAR

 35:    Level: beginner

 37: .seealso: PEPSetUp(), PEPSolve(), PEPDestroy(), PEP
 38: @*/
 39: PetscErrorCode PEPCreate(MPI_Comm comm,PEP *outpep)
 40: {
 42:   PEP            pep;

 46:   *outpep = 0;
 47:   PEPInitializePackage();
 48:   SlepcHeaderCreate(pep,PEP_CLASSID,"PEP","Polynomial Eigenvalue Problem","PEP",comm,PEPDestroy,PEPView);

 50:   pep->max_it          = 0;
 51:   pep->nev             = 1;
 52:   pep->ncv             = 0;
 53:   pep->mpd             = 0;
 54:   pep->nini            = 0;
 55:   pep->target          = 0.0;
 56:   pep->tol             = PETSC_DEFAULT;
 57:   pep->conv            = PEP_CONV_REL;
 58:   pep->stop            = PEP_STOP_BASIC;
 59:   pep->which           = (PEPWhich)0;
 60:   pep->basis           = PEP_BASIS_MONOMIAL;
 61:   pep->problem_type    = (PEPProblemType)0;
 62:   pep->scale           = PEP_SCALE_NONE;
 63:   pep->sfactor         = 1.0;
 64:   pep->dsfactor        = 1.0;
 65:   pep->sits            = 5;
 66:   pep->slambda         = 1.0;
 67:   pep->refine          = PEP_REFINE_NONE;
 68:   pep->npart           = 1;
 69:   pep->rtol            = PETSC_DEFAULT;
 70:   pep->rits            = PETSC_DEFAULT;
 71:   pep->scheme          = (PEPRefineScheme)0;
 72:   pep->extract         = (PEPExtract)0;
 73:   pep->trackall        = PETSC_FALSE;

 75:   pep->converged       = PEPConvergedRelative;
 76:   pep->convergeduser   = NULL;
 77:   pep->convergeddestroy= NULL;
 78:   pep->stopping        = PEPStoppingBasic;
 79:   pep->stoppinguser    = NULL;
 80:   pep->stoppingdestroy = NULL;
 81:   pep->convergedctx    = NULL;
 82:   pep->stoppingctx     = NULL;
 83:   pep->numbermonitors  = 0;

 85:   pep->st              = NULL;
 86:   pep->ds              = NULL;
 87:   pep->V               = NULL;
 88:   pep->rg              = NULL;
 89:   pep->A               = NULL;
 90:   pep->nmat            = 0;
 91:   pep->Dl              = NULL;
 92:   pep->Dr              = NULL;
 93:   pep->IS              = NULL;
 94:   pep->eigr            = NULL;
 95:   pep->eigi            = NULL;
 96:   pep->errest          = NULL;
 97:   pep->perm            = NULL;
 98:   pep->pbc             = NULL;
 99:   pep->solvematcoeffs  = NULL;
100:   pep->nwork           = 0;
101:   pep->work            = NULL;
102:   pep->refineksp       = NULL;
103:   pep->refinesubc      = NULL;
104:   pep->data            = NULL;

106:   pep->state           = PEP_STATE_INITIAL;
107:   pep->nconv           = 0;
108:   pep->its             = 0;
109:   pep->n               = 0;
110:   pep->nloc            = 0;
111:   pep->nrma            = NULL;
112:   pep->sfactor_set     = PETSC_FALSE;
113:   pep->lineariz        = PETSC_FALSE;
114:   pep->reason          = PEP_CONVERGED_ITERATING;

116:   PetscNewLog(pep,&pep->sc);
117:   *outpep = pep;
118:   return(0);
119: }

121: /*@C
122:    PEPSetType - Selects the particular solver to be used in the PEP object.

124:    Logically Collective on PEP

126:    Input Parameters:
127: +  pep      - the polynomial eigensolver context
128: -  type     - a known method

130:    Options Database Key:
131: .  -pep_type <method> - Sets the method; use -help for a list
132:     of available methods

134:    Notes:
135:    See "slepc/include/slepcpep.h" for available methods. The default
136:    is PEPTOAR.

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

146:    Level: intermediate

148: .seealso: PEPType
149: @*/
150: PetscErrorCode PEPSetType(PEP pep,PEPType type)
151: {
152:   PetscErrorCode ierr,(*r)(PEP);
153:   PetscBool      match;


159:   PetscObjectTypeCompare((PetscObject)pep,type,&match);
160:   if (match) return(0);

162:   PetscFunctionListFind(PEPList,type,&r);
163:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown PEP type given: %s",type);

165:   if (pep->ops->destroy) { (*pep->ops->destroy)(pep); }
166:   PetscMemzero(pep->ops,sizeof(struct _PEPOps));

168:   pep->state = PEP_STATE_INITIAL;
169:   PetscObjectChangeTypeName((PetscObject)pep,type);
170:   (*r)(pep);
171:   return(0);
172: }

174: /*@C
175:    PEPGetType - Gets the PEP type as a string from the PEP object.

177:    Not Collective

179:    Input Parameter:
180: .  pep - the eigensolver context

182:    Output Parameter:
183: .  name - name of PEP method

185:    Level: intermediate

187: .seealso: PEPSetType()
188: @*/
189: PetscErrorCode PEPGetType(PEP pep,PEPType *type)
190: {
194:   *type = ((PetscObject)pep)->type_name;
195:   return(0);
196: }

198: /*@C
199:    PEPRegister - Adds a method to the polynomial eigenproblem solver package.

201:    Not Collective

203:    Input Parameters:
204: +  name - name of a new user-defined solver
205: -  function - routine to create the solver context

207:    Notes:
208:    PEPRegister() may be called multiple times to add several user-defined solvers.

210:    Sample usage:
211: .vb
212:     PEPRegister("my_solver",MySolverCreate);
213: .ve

215:    Then, your solver can be chosen with the procedural interface via
216: $     PEPSetType(pep,"my_solver")
217:    or at runtime via the option
218: $     -pep_type my_solver

220:    Level: advanced

222: .seealso: PEPRegisterAll()
223: @*/
224: PetscErrorCode PEPRegister(const char *name,PetscErrorCode (*function)(PEP))
225: {

229:   PEPInitializePackage();
230:   PetscFunctionListAdd(&PEPList,name,function);
231:   return(0);
232: }

234: /*@
235:    PEPReset - Resets the PEP context to the initial state (prior to setup)
236:    and destroys any allocated Vecs and Mats.

238:    Collective on PEP

240:    Input Parameter:
241: .  pep - eigensolver context obtained from PEPCreate()

243:    Level: advanced

245: .seealso: PEPDestroy()
246: @*/
247: PetscErrorCode PEPReset(PEP pep)
248: {

253:   if (!pep) return(0);
254:   if (pep->ops->reset) { (pep->ops->reset)(pep); }
255:   if (pep->st) { STReset(pep->st); }
256:   if (pep->refineksp) { KSPReset(pep->refineksp); }
257:   if (pep->nmat) {
258:     MatDestroyMatrices(pep->nmat,&pep->A);
259:     PetscFree2(pep->pbc,pep->nrma);
260:     PetscFree(pep->solvematcoeffs);
261:     pep->nmat = 0;
262:   }
263:   VecDestroy(&pep->Dl);
264:   VecDestroy(&pep->Dr);
265:   BVDestroy(&pep->V);
266:   VecDestroyVecs(pep->nwork,&pep->work);
267:   pep->nwork = 0;
268:   pep->state = PEP_STATE_INITIAL;
269:   return(0);
270: }

272: /*@
273:    PEPDestroy - Destroys the PEP context.

275:    Collective on PEP

277:    Input Parameter:
278: .  pep - eigensolver context obtained from PEPCreate()

280:    Level: beginner

282: .seealso: PEPCreate(), PEPSetUp(), PEPSolve()
283: @*/
284: PetscErrorCode PEPDestroy(PEP *pep)
285: {

289:   if (!*pep) return(0);
291:   if (--((PetscObject)(*pep))->refct > 0) { *pep = 0; return(0); }
292:   PEPReset(*pep);
293:   if ((*pep)->ops->destroy) { (*(*pep)->ops->destroy)(*pep); }
294:   if ((*pep)->eigr) {
295:     PetscFree4((*pep)->eigr,(*pep)->eigi,(*pep)->errest,(*pep)->perm);
296:   }
297:   STDestroy(&(*pep)->st);
298:   RGDestroy(&(*pep)->rg);
299:   DSDestroy(&(*pep)->ds);
300:   KSPDestroy(&(*pep)->refineksp);
301:   PetscSubcommDestroy(&(*pep)->refinesubc);
302:   PetscFree((*pep)->sc);
303:   /* just in case the initial vectors have not been used */
304:   SlepcBasisDestroy_Private(&(*pep)->nini,&(*pep)->IS);
305:   if ((*pep)->convergeddestroy) {
306:     (*(*pep)->convergeddestroy)((*pep)->convergedctx);
307:   }
308:   PEPMonitorCancel(*pep);
309:   PetscHeaderDestroy(pep);
310:   return(0);
311: }

313: /*@
314:    PEPSetBV - Associates a basis vectors object to the polynomial eigensolver.

316:    Collective on PEP

318:    Input Parameters:
319: +  pep - eigensolver context obtained from PEPCreate()
320: -  bv  - the basis vectors object

322:    Note:
323:    Use PEPGetBV() to retrieve the basis vectors context (for example,
324:    to free it at the end of the computations).

326:    Level: advanced

328: .seealso: PEPGetBV()
329: @*/
330: PetscErrorCode PEPSetBV(PEP pep,BV bv)
331: {

338:   PetscObjectReference((PetscObject)bv);
339:   BVDestroy(&pep->V);
340:   pep->V = bv;
341:   PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->V);
342:   return(0);
343: }

345: /*@
346:    PEPGetBV - Obtain the basis vectors object associated to the polynomial
347:    eigensolver object.

349:    Not Collective

351:    Input Parameters:
352: .  pep - eigensolver context obtained from PEPCreate()

354:    Output Parameter:
355: .  bv - basis vectors context

357:    Level: advanced

359: .seealso: PEPSetBV()
360: @*/
361: PetscErrorCode PEPGetBV(PEP pep,BV *bv)
362: {

368:   if (!pep->V) {
369:     BVCreate(PetscObjectComm((PetscObject)pep),&pep->V);
370:     PetscObjectIncrementTabLevel((PetscObject)pep->V,(PetscObject)pep,0);
371:     PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->V);
372:     PetscObjectSetOptions((PetscObject)pep->V,((PetscObject)pep)->options);
373:   }
374:   *bv = pep->V;
375:   return(0);
376: }

378: /*@
379:    PEPSetRG - Associates a region object to the polynomial eigensolver.

381:    Collective on PEP

383:    Input Parameters:
384: +  pep - eigensolver context obtained from PEPCreate()
385: -  rg  - the region object

387:    Note:
388:    Use PEPGetRG() to retrieve the region context (for example,
389:    to free it at the end of the computations).

391:    Level: advanced

393: .seealso: PEPGetRG()
394: @*/
395: PetscErrorCode PEPSetRG(PEP pep,RG rg)
396: {

403:   PetscObjectReference((PetscObject)rg);
404:   RGDestroy(&pep->rg);
405:   pep->rg = rg;
406:   PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
407:   return(0);
408: }

410: /*@
411:    PEPGetRG - Obtain the region object associated to the
412:    polynomial eigensolver object.

414:    Not Collective

416:    Input Parameters:
417: .  pep - eigensolver context obtained from PEPCreate()

419:    Output Parameter:
420: .  rg - region context

422:    Level: advanced

424: .seealso: PEPSetRG()
425: @*/
426: PetscErrorCode PEPGetRG(PEP pep,RG *rg)
427: {

433:   if (!pep->rg) {
434:     RGCreate(PetscObjectComm((PetscObject)pep),&pep->rg);
435:     PetscObjectIncrementTabLevel((PetscObject)pep->rg,(PetscObject)pep,0);
436:     PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
437:     PetscObjectSetOptions((PetscObject)pep->rg,((PetscObject)pep)->options);
438:   }
439:   *rg = pep->rg;
440:   return(0);
441: }

443: /*@
444:    PEPSetDS - Associates a direct solver object to the polynomial eigensolver.

446:    Collective on PEP

448:    Input Parameters:
449: +  pep - eigensolver context obtained from PEPCreate()
450: -  ds  - the direct solver object

452:    Note:
453:    Use PEPGetDS() to retrieve the direct solver context (for example,
454:    to free it at the end of the computations).

456:    Level: advanced

458: .seealso: PEPGetDS()
459: @*/
460: PetscErrorCode PEPSetDS(PEP pep,DS ds)
461: {

468:   PetscObjectReference((PetscObject)ds);
469:   DSDestroy(&pep->ds);
470:   pep->ds = ds;
471:   PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
472:   return(0);
473: }

475: /*@
476:    PEPGetDS - Obtain the direct solver object associated to the
477:    polynomial eigensolver object.

479:    Not Collective

481:    Input Parameters:
482: .  pep - eigensolver context obtained from PEPCreate()

484:    Output Parameter:
485: .  ds - direct solver context

487:    Level: advanced

489: .seealso: PEPSetDS()
490: @*/
491: PetscErrorCode PEPGetDS(PEP pep,DS *ds)
492: {

498:   if (!pep->ds) {
499:     DSCreate(PetscObjectComm((PetscObject)pep),&pep->ds);
500:     PetscObjectIncrementTabLevel((PetscObject)pep->ds,(PetscObject)pep,0);
501:     PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
502:     PetscObjectSetOptions((PetscObject)pep->ds,((PetscObject)pep)->options);
503:   }
504:   *ds = pep->ds;
505:   return(0);
506: }

508: /*@
509:    PEPSetST - Associates a spectral transformation object to the eigensolver.

511:    Collective on PEP

513:    Input Parameters:
514: +  pep - eigensolver context obtained from PEPCreate()
515: -  st   - the spectral transformation object

517:    Note:
518:    Use PEPGetST() to retrieve the spectral transformation context (for example,
519:    to free it at the end of the computations).

521:    Level: advanced

523: .seealso: PEPGetST()
524: @*/
525: PetscErrorCode PEPSetST(PEP pep,ST st)
526: {

533:   PetscObjectReference((PetscObject)st);
534:   STDestroy(&pep->st);
535:   pep->st = st;
536:   PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
537:   return(0);
538: }

540: /*@
541:    PEPGetST - Obtain the spectral transformation (ST) object associated
542:    to the eigensolver object.

544:    Not Collective

546:    Input Parameters:
547: .  pep - eigensolver context obtained from PEPCreate()

549:    Output Parameter:
550: .  st - spectral transformation context

552:    Level: intermediate

554: .seealso: PEPSetST()
555: @*/
556: PetscErrorCode PEPGetST(PEP pep,ST *st)
557: {

563:   if (!pep->st) {
564:     STCreate(PetscObjectComm((PetscObject)pep),&pep->st);
565:     PetscObjectIncrementTabLevel((PetscObject)pep->st,(PetscObject)pep,0);
566:     PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
567:     PetscObjectSetOptions((PetscObject)pep->st,((PetscObject)pep)->options);
568:   }
569:   *st = pep->st;
570:   return(0);
571: }

573: /*@
574:    PEPRefineGetKSP - Obtain the ksp object used by the eigensolver
575:    object in the refinement phase.

577:    Not Collective

579:    Input Parameters:
580: .  pep - eigensolver context obtained from PEPCreate()

582:    Output Parameter:
583: .  ksp - ksp context

585:    Level: advanced

587: .seealso: PEPSetRefine()
588: @*/
589: PetscErrorCode PEPRefineGetKSP(PEP pep,KSP *ksp)
590: {

596:   if (!pep->refineksp) {
597:     if (pep->npart>1) {
598:       /* Split in subcomunicators */
599:       PetscSubcommCreate(PetscObjectComm((PetscObject)pep),&pep->refinesubc);
600:       PetscSubcommSetNumber(pep->refinesubc,pep->npart);
601:       PetscSubcommSetType(pep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
602:       PetscLogObjectMemory((PetscObject)pep,sizeof(PetscSubcomm));
603:     }
604:     KSPCreate((pep->npart==1)?PetscObjectComm((PetscObject)pep):PetscSubcommChild(pep->refinesubc),&pep->refineksp);
605:     PetscObjectIncrementTabLevel((PetscObject)pep->refineksp,(PetscObject)pep,0);
606:     PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->refineksp);
607:     PetscObjectSetOptions((PetscObject)pep->refineksp,((PetscObject)pep)->options);
608:     KSPSetOptionsPrefix(*ksp,((PetscObject)pep)->prefix);
609:     KSPAppendOptionsPrefix(*ksp,"pep_refine_");
610:     KSPSetTolerances(pep->refineksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
611:   }
612:   *ksp = pep->refineksp;
613:   return(0);
614: }

616: /*@
617:    PEPSetTarget - Sets the value of the target.

619:    Logically Collective on PEP

621:    Input Parameters:
622: +  pep    - eigensolver context
623: -  target - the value of the target

625:    Options Database Key:
626: .  -pep_target <scalar> - the value of the target

628:    Notes:
629:    The target is a scalar value used to determine the portion of the spectrum
630:    of interest. It is used in combination with PEPSetWhichEigenpairs().

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

636:    Level: intermediate

638: .seealso: PEPGetTarget(), PEPSetWhichEigenpairs()
639: @*/
640: PetscErrorCode PEPSetTarget(PEP pep,PetscScalar target)
641: {

647:   pep->target = target;
648:   if (!pep->st) { PEPGetST(pep,&pep->st); }
649:   STSetDefaultShift(pep->st,target);
650:   return(0);
651: }

653: /*@
654:    PEPGetTarget - Gets the value of the target.

656:    Not Collective

658:    Input Parameter:
659: .  pep - eigensolver context

661:    Output Parameter:
662: .  target - the value of the target

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

667:    Level: intermediate

669: .seealso: PEPSetTarget()
670: @*/
671: PetscErrorCode PEPGetTarget(PEP pep,PetscScalar* target)
672: {
676:   *target = pep->target;
677:   return(0);
678: }

680: /*@
681:    PEPSetInterval - Defines the computational interval for spectrum slicing.

683:    Logically Collective on PEP

685:    Input Parameters:
686: +  pep  - eigensolver context
687: .  inta - left end of the interval
688: -  intb - right end of the interval

690:    Options Database Key:
691: .  -pep_interval <a,b> - set [a,b] as the interval of interest

693:    Notes:
694:    Spectrum slicing is a technique employed for computing all eigenvalues of
695:    symmetric eigenproblems in a given interval. This function provides the
696:    interval to be considered. It must be used in combination with PEP_ALL, see
697:    PEPSetWhichEigenpairs().

699:    In the command-line option, two values must be provided. For an open interval,
700:    one can give an infinite, e.g., -pep_interval 1.0,inf or -pep_interval -inf,1.0.
701:    An open interval in the programmatic interface can be specified with
702:    PETSC_MAX_REAL and -PETSC_MAX_REAL.

704:    Level: intermediate

706: .seealso: PEPGetInterval(), PEPSetWhichEigenpairs()
707: @*/
708: PetscErrorCode PEPSetInterval(PEP pep,PetscReal inta,PetscReal intb)
709: {
714:   if (inta>intb) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
715:   if (pep->inta != inta || pep->intb != intb) {
716:     pep->inta = inta;
717:     pep->intb = intb;
718:     pep->state = PEP_STATE_INITIAL;
719:   }
720:   return(0);
721: }

723: /*@
724:    PEPGetInterval - Gets the computational interval for spectrum slicing.

726:    Not Collective

728:    Input Parameter:
729: .  pep - eigensolver context

731:    Output Parameters:
732: +  inta - left end of the interval
733: -  intb - right end of the interval

735:    Level: intermediate

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

740: .seealso: PEPSetInterval()
741: @*/
742: PetscErrorCode PEPGetInterval(PEP pep,PetscReal* inta,PetscReal* intb)
743: {
746:   if (inta) *inta = pep->inta;
747:   if (intb) *intb = pep->intb;
748:   return(0);
749: }