Actual source code: pepbasic.c

slepc-3.18.3 2023-03-24
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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>

 16: /* Logging support */
 17: PetscClassId      PEP_CLASSID = 0;
 18: PetscLogEvent     PEP_SetUp = 0,PEP_Solve = 0,PEP_Refine = 0,PEP_CISS_SVD = 0;

 20: /* List of registered PEP routines */
 21: PetscFunctionList PEPList = NULL;
 22: PetscBool         PEPRegisterAllCalled = PETSC_FALSE;

 24: /* List of registered PEP monitors */
 25: PetscFunctionList PEPMonitorList              = NULL;
 26: PetscFunctionList PEPMonitorCreateList        = NULL;
 27: PetscFunctionList PEPMonitorDestroyList       = NULL;
 28: PetscBool         PEPMonitorRegisterAllCalled = PETSC_FALSE;

 30: /*@
 31:    PEPCreate - Creates the default PEP context.

 33:    Collective

 35:    Input Parameter:
 36: .  comm - MPI communicator

 38:    Output Parameter:
 39: .  outpep - location to put the PEP context

 41:    Note:
 42:    The default PEP type is PEPTOAR

 44:    Level: beginner

 46: .seealso: PEPSetUp(), PEPSolve(), PEPDestroy(), PEP
 47: @*/
 48: PetscErrorCode PEPCreate(MPI_Comm comm,PEP *outpep)
 49: {
 50:   PEP            pep;

 53:   *outpep = NULL;
 54:   PEPInitializePackage();
 55:   SlepcHeaderCreate(pep,PEP_CLASSID,"PEP","Polynomial Eigenvalue Problem","PEP",comm,PEPDestroy,PEPView);

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

 82:   pep->converged       = PEPConvergedRelative;
 83:   pep->convergeduser   = NULL;
 84:   pep->convergeddestroy= NULL;
 85:   pep->stopping        = PEPStoppingBasic;
 86:   pep->stoppinguser    = NULL;
 87:   pep->stoppingdestroy = NULL;
 88:   pep->convergedctx    = NULL;
 89:   pep->stoppingctx     = NULL;
 90:   pep->numbermonitors  = 0;

 92:   pep->st              = NULL;
 93:   pep->ds              = NULL;
 94:   pep->V               = NULL;
 95:   pep->rg              = NULL;
 96:   pep->A               = NULL;
 97:   pep->nmat            = 0;
 98:   pep->Dl              = NULL;
 99:   pep->Dr              = NULL;
100:   pep->IS              = NULL;
101:   pep->eigr            = NULL;
102:   pep->eigi            = NULL;
103:   pep->errest          = NULL;
104:   pep->perm            = NULL;
105:   pep->pbc             = NULL;
106:   pep->solvematcoeffs  = NULL;
107:   pep->nwork           = 0;
108:   pep->work            = NULL;
109:   pep->refineksp       = NULL;
110:   pep->refinesubc      = NULL;
111:   pep->data            = NULL;

113:   pep->state           = PEP_STATE_INITIAL;
114:   pep->nconv           = 0;
115:   pep->its             = 0;
116:   pep->n               = 0;
117:   pep->nloc            = 0;
118:   pep->nrma            = NULL;
119:   pep->sfactor_set     = PETSC_FALSE;
120:   pep->lineariz        = PETSC_FALSE;
121:   pep->reason          = PEP_CONVERGED_ITERATING;

123:   PetscNew(&pep->sc);
124:   *outpep = pep;
125:   return 0;
126: }

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

131:    Logically Collective on pep

133:    Input Parameters:
134: +  pep      - the polynomial eigensolver context
135: -  type     - a known method

137:    Options Database Key:
138: .  -pep_type <method> - Sets the method; use -help for a list
139:     of available methods

141:    Notes:
142:    See "slepc/include/slepcpep.h" for available methods. The default
143:    is PEPTOAR.

145:    Normally, it is best to use the PEPSetFromOptions() command and
146:    then set the PEP type from the options database rather than by using
147:    this routine.  Using the options database provides the user with
148:    maximum flexibility in evaluating the different available methods.
149:    The PEPSetType() routine is provided for those situations where it
150:    is necessary to set the iterative solver independently of the command
151:    line or options database.

153:    Level: intermediate

155: .seealso: PEPType
156: @*/
157: PetscErrorCode PEPSetType(PEP pep,PEPType type)
158: {
159:   PetscErrorCode (*r)(PEP);
160:   PetscBool      match;


165:   PetscObjectTypeCompare((PetscObject)pep,type,&match);
166:   if (match) return 0;

168:   PetscFunctionListFind(PEPList,type,&r);

171:   PetscTryTypeMethod(pep,destroy);
172:   PetscMemzero(pep->ops,sizeof(struct _PEPOps));

174:   pep->state = PEP_STATE_INITIAL;
175:   PetscObjectChangeTypeName((PetscObject)pep,type);
176:   (*r)(pep);
177:   return 0;
178: }

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

183:    Not Collective

185:    Input Parameter:
186: .  pep - the eigensolver context

188:    Output Parameter:
189: .  type - name of PEP method

191:    Level: intermediate

193: .seealso: PEPSetType()
194: @*/
195: PetscErrorCode PEPGetType(PEP pep,PEPType *type)
196: {
199:   *type = ((PetscObject)pep)->type_name;
200:   return 0;
201: }

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

206:    Not Collective

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

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

215:    Sample usage:
216: .vb
217:     PEPRegister("my_solver",MySolverCreate);
218: .ve

220:    Then, your solver can be chosen with the procedural interface via
221: $     PEPSetType(pep,"my_solver")
222:    or at runtime via the option
223: $     -pep_type my_solver

225:    Level: advanced

227: .seealso: PEPRegisterAll()
228: @*/
229: PetscErrorCode PEPRegister(const char *name,PetscErrorCode (*function)(PEP))
230: {
231:   PEPInitializePackage();
232:   PetscFunctionListAdd(&PEPList,name,function);
233:   return 0;
234: }

236: /*@C
237:    PEPMonitorRegister - Adds PEP monitor routine.

239:    Not Collective

241:    Input Parameters:
242: +  name    - name of a new monitor routine
243: .  vtype   - a PetscViewerType for the output
244: .  format  - a PetscViewerFormat for the output
245: .  monitor - monitor routine
246: .  create  - creation routine, or NULL
247: -  destroy - destruction routine, or NULL

249:    Notes:
250:    PEPMonitorRegister() may be called multiple times to add several user-defined monitors.

252:    Sample usage:
253: .vb
254:    PEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
255: .ve

257:    Then, your monitor can be chosen with the procedural interface via
258: $      PEPMonitorSetFromOptions(pep,"-pep_monitor_my_monitor","my_monitor",NULL)
259:    or at runtime via the option
260: $      -pep_monitor_my_monitor

262:    Level: advanced

264: .seealso: PEPMonitorRegisterAll()
265: @*/
266: PetscErrorCode PEPMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(PEP,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscReal*,PetscInt,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
267: {
268:   char           key[PETSC_MAX_PATH_LEN];

270:   PEPInitializePackage();
271:   SlepcMonitorMakeKey_Internal(name,vtype,format,key);
272:   PetscFunctionListAdd(&PEPMonitorList,key,monitor);
273:   if (create)  PetscFunctionListAdd(&PEPMonitorCreateList,key,create);
274:   if (destroy) PetscFunctionListAdd(&PEPMonitorDestroyList,key,destroy);
275:   return 0;
276: }

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

282:    Collective on pep

284:    Input Parameter:
285: .  pep - eigensolver context obtained from PEPCreate()

287:    Level: advanced

289: .seealso: PEPDestroy()
290: @*/
291: PetscErrorCode PEPReset(PEP pep)
292: {
294:   if (!pep) return 0;
295:   PetscTryTypeMethod(pep,reset);
296:   if (pep->st) STReset(pep->st);
297:   if (pep->refineksp) KSPReset(pep->refineksp);
298:   if (pep->nmat) {
299:     MatDestroyMatrices(pep->nmat,&pep->A);
300:     PetscFree2(pep->pbc,pep->nrma);
301:     PetscFree(pep->solvematcoeffs);
302:     pep->nmat = 0;
303:   }
304:   VecDestroy(&pep->Dl);
305:   VecDestroy(&pep->Dr);
306:   BVDestroy(&pep->V);
307:   VecDestroyVecs(pep->nwork,&pep->work);
308:   pep->nwork = 0;
309:   pep->state = PEP_STATE_INITIAL;
310:   return 0;
311: }

313: /*@C
314:    PEPDestroy - Destroys the PEP context.

316:    Collective on pep

318:    Input Parameter:
319: .  pep - eigensolver context obtained from PEPCreate()

321:    Level: beginner

323: .seealso: PEPCreate(), PEPSetUp(), PEPSolve()
324: @*/
325: PetscErrorCode PEPDestroy(PEP *pep)
326: {
327:   if (!*pep) return 0;
329:   if (--((PetscObject)(*pep))->refct > 0) { *pep = NULL; return 0; }
330:   PEPReset(*pep);
331:   PetscTryTypeMethod(*pep,destroy);
332:   if ((*pep)->eigr) PetscFree4((*pep)->eigr,(*pep)->eigi,(*pep)->errest,(*pep)->perm);
333:   STDestroy(&(*pep)->st);
334:   RGDestroy(&(*pep)->rg);
335:   DSDestroy(&(*pep)->ds);
336:   KSPDestroy(&(*pep)->refineksp);
337:   PetscSubcommDestroy(&(*pep)->refinesubc);
338:   PetscFree((*pep)->sc);
339:   /* just in case the initial vectors have not been used */
340:   SlepcBasisDestroy_Private(&(*pep)->nini,&(*pep)->IS);
341:   if ((*pep)->convergeddestroy) (*(*pep)->convergeddestroy)((*pep)->convergedctx);
342:   PEPMonitorCancel(*pep);
343:   PetscHeaderDestroy(pep);
344:   return 0;
345: }

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

350:    Collective on pep

352:    Input Parameters:
353: +  pep - eigensolver context obtained from PEPCreate()
354: -  bv  - the basis vectors object

356:    Note:
357:    Use PEPGetBV() to retrieve the basis vectors context (for example,
358:    to free it at the end of the computations).

360:    Level: advanced

362: .seealso: PEPGetBV()
363: @*/
364: PetscErrorCode PEPSetBV(PEP pep,BV bv)
365: {
369:   PetscObjectReference((PetscObject)bv);
370:   BVDestroy(&pep->V);
371:   pep->V = bv;
372:   return 0;
373: }

375: /*@
376:    PEPGetBV - Obtain the basis vectors object associated to the polynomial
377:    eigensolver object.

379:    Not Collective

381:    Input Parameters:
382: .  pep - eigensolver context obtained from PEPCreate()

384:    Output Parameter:
385: .  bv - basis vectors context

387:    Level: advanced

389: .seealso: PEPSetBV()
390: @*/
391: PetscErrorCode PEPGetBV(PEP pep,BV *bv)
392: {
395:   if (!pep->V) {
396:     BVCreate(PetscObjectComm((PetscObject)pep),&pep->V);
397:     PetscObjectIncrementTabLevel((PetscObject)pep->V,(PetscObject)pep,0);
398:     PetscObjectSetOptions((PetscObject)pep->V,((PetscObject)pep)->options);
399:   }
400:   *bv = pep->V;
401:   return 0;
402: }

404: /*@
405:    PEPSetRG - Associates a region object to the polynomial eigensolver.

407:    Collective on pep

409:    Input Parameters:
410: +  pep - eigensolver context obtained from PEPCreate()
411: -  rg  - the region object

413:    Note:
414:    Use PEPGetRG() to retrieve the region context (for example,
415:    to free it at the end of the computations).

417:    Level: advanced

419: .seealso: PEPGetRG()
420: @*/
421: PetscErrorCode PEPSetRG(PEP pep,RG rg)
422: {
424:   if (rg) {
427:   }
428:   PetscObjectReference((PetscObject)rg);
429:   RGDestroy(&pep->rg);
430:   pep->rg = rg;
431:   return 0;
432: }

434: /*@
435:    PEPGetRG - Obtain the region object associated to the
436:    polynomial eigensolver object.

438:    Not Collective

440:    Input Parameters:
441: .  pep - eigensolver context obtained from PEPCreate()

443:    Output Parameter:
444: .  rg - region context

446:    Level: advanced

448: .seealso: PEPSetRG()
449: @*/
450: PetscErrorCode PEPGetRG(PEP pep,RG *rg)
451: {
454:   if (!pep->rg) {
455:     RGCreate(PetscObjectComm((PetscObject)pep),&pep->rg);
456:     PetscObjectIncrementTabLevel((PetscObject)pep->rg,(PetscObject)pep,0);
457:     PetscObjectSetOptions((PetscObject)pep->rg,((PetscObject)pep)->options);
458:   }
459:   *rg = pep->rg;
460:   return 0;
461: }

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

466:    Collective on pep

468:    Input Parameters:
469: +  pep - eigensolver context obtained from PEPCreate()
470: -  ds  - the direct solver object

472:    Note:
473:    Use PEPGetDS() to retrieve the direct solver context (for example,
474:    to free it at the end of the computations).

476:    Level: advanced

478: .seealso: PEPGetDS()
479: @*/
480: PetscErrorCode PEPSetDS(PEP pep,DS ds)
481: {
485:   PetscObjectReference((PetscObject)ds);
486:   DSDestroy(&pep->ds);
487:   pep->ds = ds;
488:   return 0;
489: }

491: /*@
492:    PEPGetDS - Obtain the direct solver object associated to the
493:    polynomial eigensolver object.

495:    Not Collective

497:    Input Parameters:
498: .  pep - eigensolver context obtained from PEPCreate()

500:    Output Parameter:
501: .  ds - direct solver context

503:    Level: advanced

505: .seealso: PEPSetDS()
506: @*/
507: PetscErrorCode PEPGetDS(PEP pep,DS *ds)
508: {
511:   if (!pep->ds) {
512:     DSCreate(PetscObjectComm((PetscObject)pep),&pep->ds);
513:     PetscObjectIncrementTabLevel((PetscObject)pep->ds,(PetscObject)pep,0);
514:     PetscObjectSetOptions((PetscObject)pep->ds,((PetscObject)pep)->options);
515:   }
516:   *ds = pep->ds;
517:   return 0;
518: }

520: /*@
521:    PEPSetST - Associates a spectral transformation object to the eigensolver.

523:    Collective on pep

525:    Input Parameters:
526: +  pep - eigensolver context obtained from PEPCreate()
527: -  st   - the spectral transformation object

529:    Note:
530:    Use PEPGetST() to retrieve the spectral transformation context (for example,
531:    to free it at the end of the computations).

533:    Level: advanced

535: .seealso: PEPGetST()
536: @*/
537: PetscErrorCode PEPSetST(PEP pep,ST st)
538: {
542:   PetscObjectReference((PetscObject)st);
543:   STDestroy(&pep->st);
544:   pep->st = st;
545:   return 0;
546: }

548: /*@
549:    PEPGetST - Obtain the spectral transformation (ST) object associated
550:    to the eigensolver object.

552:    Not Collective

554:    Input Parameters:
555: .  pep - eigensolver context obtained from PEPCreate()

557:    Output Parameter:
558: .  st - spectral transformation context

560:    Level: intermediate

562: .seealso: PEPSetST()
563: @*/
564: PetscErrorCode PEPGetST(PEP pep,ST *st)
565: {
568:   if (!pep->st) {
569:     STCreate(PetscObjectComm((PetscObject)pep),&pep->st);
570:     PetscObjectIncrementTabLevel((PetscObject)pep->st,(PetscObject)pep,0);
571:     PetscObjectSetOptions((PetscObject)pep->st,((PetscObject)pep)->options);
572:   }
573:   *st = pep->st;
574:   return 0;
575: }

577: /*@
578:    PEPRefineGetKSP - Obtain the ksp object used by the eigensolver
579:    object in the refinement phase.

581:    Not Collective

583:    Input Parameters:
584: .  pep - eigensolver context obtained from PEPCreate()

586:    Output Parameter:
587: .  ksp - ksp context

589:    Level: advanced

591: .seealso: PEPSetRefine()
592: @*/
593: PetscErrorCode PEPRefineGetKSP(PEP pep,KSP *ksp)
594: {
595:   MPI_Comm       comm;

599:   if (!pep->refineksp) {
600:     if (pep->npart>1) {
601:       /* Split in subcomunicators */
602:       PetscSubcommCreate(PetscObjectComm((PetscObject)pep),&pep->refinesubc);
603:       PetscSubcommSetNumber(pep->refinesubc,pep->npart);
604:       PetscSubcommSetType(pep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
605:       PetscSubcommGetChild(pep->refinesubc,&comm);
606:     } else PetscObjectGetComm((PetscObject)pep,&comm);
607:     KSPCreate(comm,&pep->refineksp);
608:     PetscObjectIncrementTabLevel((PetscObject)pep->refineksp,(PetscObject)pep,0);
609:     PetscObjectSetOptions((PetscObject)pep->refineksp,((PetscObject)pep)->options);
610:     KSPSetOptionsPrefix(*ksp,((PetscObject)pep)->prefix);
611:     KSPAppendOptionsPrefix(*ksp,"pep_refine_");
612:     KSPSetTolerances(pep->refineksp,SlepcDefaultTol(pep->rtol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
613:   }
614:   *ksp = pep->refineksp;
615:   return 0;
616: }

618: /*@
619:    PEPSetTarget - Sets the value of the target.

621:    Logically Collective on pep

623:    Input Parameters:
624: +  pep    - eigensolver context
625: -  target - the value of the target

627:    Options Database Key:
628: .  -pep_target <scalar> - the value of the target

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

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

638:    Level: intermediate

640: .seealso: PEPGetTarget(), PEPSetWhichEigenpairs()
641: @*/
642: PetscErrorCode PEPSetTarget(PEP pep,PetscScalar target)
643: {
646:   pep->target = target;
647:   if (!pep->st) PEPGetST(pep,&pep->st);
648:   STSetDefaultShift(pep->st,target);
649:   return 0;
650: }

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

655:    Not Collective

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

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

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

666:    Level: intermediate

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

678: /*@
679:    PEPSetInterval - Defines the computational interval for spectrum slicing.

681:    Logically Collective on pep

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

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

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

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

702:    Level: intermediate

704: .seealso: PEPGetInterval(), PEPSetWhichEigenpairs()
705: @*/
706: PetscErrorCode PEPSetInterval(PEP pep,PetscReal inta,PetscReal intb)
707: {
712:   if (pep->inta != inta || pep->intb != intb) {
713:     pep->inta = inta;
714:     pep->intb = intb;
715:     pep->state = PEP_STATE_INITIAL;
716:   }
717:   return 0;
718: }

720: /*@
721:    PEPGetInterval - Gets the computational interval for spectrum slicing.

723:    Not Collective

725:    Input Parameter:
726: .  pep - eigensolver context

728:    Output Parameters:
729: +  inta - left end of the interval
730: -  intb - right end of the interval

732:    Level: intermediate

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

737: .seealso: PEPSetInterval()
738: @*/
739: PetscErrorCode PEPGetInterval(PEP pep,PetscReal* inta,PetscReal* intb)
740: {
742:   if (inta) *inta = pep->inta;
743:   if (intb) *intb = pep->intb;
744:   return 0;
745: }