Actual source code: pepbasic.c
slepc-3.17.2 2022-08-09
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 = 0;
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 = 0;
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: PetscNewLog(pep,&pep->sc);
124: *outpep = pep;
125: PetscFunctionReturn(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) PetscFunctionReturn(0);
168: PetscFunctionListFind(PEPList,type,&r);
171: if (pep->ops->destroy) (*pep->ops->destroy)(pep);
172: PetscMemzero(pep->ops,sizeof(struct _PEPOps));
174: pep->state = PEP_STATE_INITIAL;
175: PetscObjectChangeTypeName((PetscObject)pep,type);
176: (*r)(pep);
177: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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: PetscFunctionReturn(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) PetscFunctionReturn(0);
295: if (pep->ops->reset) (pep->ops->reset)(pep);
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: PetscFunctionReturn(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) PetscFunctionReturn(0);
329: if (--((PetscObject)(*pep))->refct > 0) { *pep = 0; PetscFunctionReturn(0); }
330: PEPReset(*pep);
331: if ((*pep)->ops->destroy) (*(*pep)->ops->destroy)(*pep);
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: PetscFunctionReturn(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: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->V);
373: PetscFunctionReturn(0);
374: }
376: /*@
377: PEPGetBV - Obtain the basis vectors object associated to the polynomial
378: eigensolver object.
380: Not Collective
382: Input Parameters:
383: . pep - eigensolver context obtained from PEPCreate()
385: Output Parameter:
386: . bv - basis vectors context
388: Level: advanced
390: .seealso: PEPSetBV()
391: @*/
392: PetscErrorCode PEPGetBV(PEP pep,BV *bv)
393: {
396: if (!pep->V) {
397: BVCreate(PetscObjectComm((PetscObject)pep),&pep->V);
398: PetscObjectIncrementTabLevel((PetscObject)pep->V,(PetscObject)pep,0);
399: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->V);
400: PetscObjectSetOptions((PetscObject)pep->V,((PetscObject)pep)->options);
401: }
402: *bv = pep->V;
403: PetscFunctionReturn(0);
404: }
406: /*@
407: PEPSetRG - Associates a region object to the polynomial eigensolver.
409: Collective on pep
411: Input Parameters:
412: + pep - eigensolver context obtained from PEPCreate()
413: - rg - the region object
415: Note:
416: Use PEPGetRG() to retrieve the region context (for example,
417: to free it at the end of the computations).
419: Level: advanced
421: .seealso: PEPGetRG()
422: @*/
423: PetscErrorCode PEPSetRG(PEP pep,RG rg)
424: {
426: if (rg) {
429: }
430: PetscObjectReference((PetscObject)rg);
431: RGDestroy(&pep->rg);
432: pep->rg = rg;
433: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
434: PetscFunctionReturn(0);
435: }
437: /*@
438: PEPGetRG - Obtain the region object associated to the
439: polynomial eigensolver object.
441: Not Collective
443: Input Parameters:
444: . pep - eigensolver context obtained from PEPCreate()
446: Output Parameter:
447: . rg - region context
449: Level: advanced
451: .seealso: PEPSetRG()
452: @*/
453: PetscErrorCode PEPGetRG(PEP pep,RG *rg)
454: {
457: if (!pep->rg) {
458: RGCreate(PetscObjectComm((PetscObject)pep),&pep->rg);
459: PetscObjectIncrementTabLevel((PetscObject)pep->rg,(PetscObject)pep,0);
460: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
461: PetscObjectSetOptions((PetscObject)pep->rg,((PetscObject)pep)->options);
462: }
463: *rg = pep->rg;
464: PetscFunctionReturn(0);
465: }
467: /*@
468: PEPSetDS - Associates a direct solver object to the polynomial eigensolver.
470: Collective on pep
472: Input Parameters:
473: + pep - eigensolver context obtained from PEPCreate()
474: - ds - the direct solver object
476: Note:
477: Use PEPGetDS() to retrieve the direct solver context (for example,
478: to free it at the end of the computations).
480: Level: advanced
482: .seealso: PEPGetDS()
483: @*/
484: PetscErrorCode PEPSetDS(PEP pep,DS ds)
485: {
489: PetscObjectReference((PetscObject)ds);
490: DSDestroy(&pep->ds);
491: pep->ds = ds;
492: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
493: PetscFunctionReturn(0);
494: }
496: /*@
497: PEPGetDS - Obtain the direct solver object associated to the
498: polynomial eigensolver object.
500: Not Collective
502: Input Parameters:
503: . pep - eigensolver context obtained from PEPCreate()
505: Output Parameter:
506: . ds - direct solver context
508: Level: advanced
510: .seealso: PEPSetDS()
511: @*/
512: PetscErrorCode PEPGetDS(PEP pep,DS *ds)
513: {
516: if (!pep->ds) {
517: DSCreate(PetscObjectComm((PetscObject)pep),&pep->ds);
518: PetscObjectIncrementTabLevel((PetscObject)pep->ds,(PetscObject)pep,0);
519: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
520: PetscObjectSetOptions((PetscObject)pep->ds,((PetscObject)pep)->options);
521: }
522: *ds = pep->ds;
523: PetscFunctionReturn(0);
524: }
526: /*@
527: PEPSetST - Associates a spectral transformation object to the eigensolver.
529: Collective on pep
531: Input Parameters:
532: + pep - eigensolver context obtained from PEPCreate()
533: - st - the spectral transformation object
535: Note:
536: Use PEPGetST() to retrieve the spectral transformation context (for example,
537: to free it at the end of the computations).
539: Level: advanced
541: .seealso: PEPGetST()
542: @*/
543: PetscErrorCode PEPSetST(PEP pep,ST st)
544: {
548: PetscObjectReference((PetscObject)st);
549: STDestroy(&pep->st);
550: pep->st = st;
551: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
552: PetscFunctionReturn(0);
553: }
555: /*@
556: PEPGetST - Obtain the spectral transformation (ST) object associated
557: to the eigensolver object.
559: Not Collective
561: Input Parameters:
562: . pep - eigensolver context obtained from PEPCreate()
564: Output Parameter:
565: . st - spectral transformation context
567: Level: intermediate
569: .seealso: PEPSetST()
570: @*/
571: PetscErrorCode PEPGetST(PEP pep,ST *st)
572: {
575: if (!pep->st) {
576: STCreate(PetscObjectComm((PetscObject)pep),&pep->st);
577: PetscObjectIncrementTabLevel((PetscObject)pep->st,(PetscObject)pep,0);
578: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
579: PetscObjectSetOptions((PetscObject)pep->st,((PetscObject)pep)->options);
580: }
581: *st = pep->st;
582: PetscFunctionReturn(0);
583: }
585: /*@
586: PEPRefineGetKSP - Obtain the ksp object used by the eigensolver
587: object in the refinement phase.
589: Not Collective
591: Input Parameters:
592: . pep - eigensolver context obtained from PEPCreate()
594: Output Parameter:
595: . ksp - ksp context
597: Level: advanced
599: .seealso: PEPSetRefine()
600: @*/
601: PetscErrorCode PEPRefineGetKSP(PEP pep,KSP *ksp)
602: {
603: MPI_Comm comm;
607: if (!pep->refineksp) {
608: if (pep->npart>1) {
609: /* Split in subcomunicators */
610: PetscSubcommCreate(PetscObjectComm((PetscObject)pep),&pep->refinesubc);
611: PetscSubcommSetNumber(pep->refinesubc,pep->npart);
612: PetscSubcommSetType(pep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
613: PetscLogObjectMemory((PetscObject)pep,sizeof(PetscSubcomm));
614: PetscSubcommGetChild(pep->refinesubc,&comm);
615: } else PetscObjectGetComm((PetscObject)pep,&comm);
616: KSPCreate(comm,&pep->refineksp);
617: PetscObjectIncrementTabLevel((PetscObject)pep->refineksp,(PetscObject)pep,0);
618: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->refineksp);
619: PetscObjectSetOptions((PetscObject)pep->refineksp,((PetscObject)pep)->options);
620: KSPSetOptionsPrefix(*ksp,((PetscObject)pep)->prefix);
621: KSPAppendOptionsPrefix(*ksp,"pep_refine_");
622: KSPSetTolerances(pep->refineksp,SlepcDefaultTol(pep->rtol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
623: }
624: *ksp = pep->refineksp;
625: PetscFunctionReturn(0);
626: }
628: /*@
629: PEPSetTarget - Sets the value of the target.
631: Logically Collective on pep
633: Input Parameters:
634: + pep - eigensolver context
635: - target - the value of the target
637: Options Database Key:
638: . -pep_target <scalar> - the value of the target
640: Notes:
641: The target is a scalar value used to determine the portion of the spectrum
642: of interest. It is used in combination with PEPSetWhichEigenpairs().
644: In the case of complex scalars, a complex value can be provided in the
645: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
646: -pep_target 1.0+2.0i
648: Level: intermediate
650: .seealso: PEPGetTarget(), PEPSetWhichEigenpairs()
651: @*/
652: PetscErrorCode PEPSetTarget(PEP pep,PetscScalar target)
653: {
656: pep->target = target;
657: if (!pep->st) PEPGetST(pep,&pep->st);
658: STSetDefaultShift(pep->st,target);
659: PetscFunctionReturn(0);
660: }
662: /*@
663: PEPGetTarget - Gets the value of the target.
665: Not Collective
667: Input Parameter:
668: . pep - eigensolver context
670: Output Parameter:
671: . target - the value of the target
673: Note:
674: If the target was not set by the user, then zero is returned.
676: Level: intermediate
678: .seealso: PEPSetTarget()
679: @*/
680: PetscErrorCode PEPGetTarget(PEP pep,PetscScalar* target)
681: {
684: *target = pep->target;
685: PetscFunctionReturn(0);
686: }
688: /*@
689: PEPSetInterval - Defines the computational interval for spectrum slicing.
691: Logically Collective on pep
693: Input Parameters:
694: + pep - eigensolver context
695: . inta - left end of the interval
696: - intb - right end of the interval
698: Options Database Key:
699: . -pep_interval <a,b> - set [a,b] as the interval of interest
701: Notes:
702: Spectrum slicing is a technique employed for computing all eigenvalues of
703: symmetric eigenproblems in a given interval. This function provides the
704: interval to be considered. It must be used in combination with PEP_ALL, see
705: PEPSetWhichEigenpairs().
707: In the command-line option, two values must be provided. For an open interval,
708: one can give an infinite, e.g., -pep_interval 1.0,inf or -pep_interval -inf,1.0.
709: An open interval in the programmatic interface can be specified with
710: PETSC_MAX_REAL and -PETSC_MAX_REAL.
712: Level: intermediate
714: .seealso: PEPGetInterval(), PEPSetWhichEigenpairs()
715: @*/
716: PetscErrorCode PEPSetInterval(PEP pep,PetscReal inta,PetscReal intb)
717: {
722: if (pep->inta != inta || pep->intb != intb) {
723: pep->inta = inta;
724: pep->intb = intb;
725: pep->state = PEP_STATE_INITIAL;
726: }
727: PetscFunctionReturn(0);
728: }
730: /*@
731: PEPGetInterval - Gets the computational interval for spectrum slicing.
733: Not Collective
735: Input Parameter:
736: . pep - eigensolver context
738: Output Parameters:
739: + inta - left end of the interval
740: - intb - right end of the interval
742: Level: intermediate
744: Note:
745: If the interval was not set by the user, then zeros are returned.
747: .seealso: PEPSetInterval()
748: @*/
749: PetscErrorCode PEPGetInterval(PEP pep,PetscReal* inta,PetscReal* intb)
750: {
752: if (inta) *inta = pep->inta;
753: if (intb) *intb = pep->intb;
754: PetscFunctionReturn(0);
755: }