Actual source code: pepbasic.c
slepc-3.18.3 2023-03-24
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: }