Actual source code: pepbasic.c
slepc-3.16.2 2022-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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: . pep - 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: {
51: PEP pep;
55: *outpep = 0;
56: PEPInitializePackage();
57: SlepcHeaderCreate(pep,PEP_CLASSID,"PEP","Polynomial Eigenvalue Problem","PEP",comm,PEPDestroy,PEPView);
59: pep->max_it = PETSC_DEFAULT;
60: pep->nev = 1;
61: pep->ncv = PETSC_DEFAULT;
62: pep->mpd = PETSC_DEFAULT;
63: pep->nini = 0;
64: pep->target = 0.0;
65: pep->tol = PETSC_DEFAULT;
66: pep->conv = PEP_CONV_REL;
67: pep->stop = PEP_STOP_BASIC;
68: pep->which = (PEPWhich)0;
69: pep->basis = PEP_BASIS_MONOMIAL;
70: pep->problem_type = (PEPProblemType)0;
71: pep->scale = PEP_SCALE_NONE;
72: pep->sfactor = 1.0;
73: pep->dsfactor = 1.0;
74: pep->sits = 5;
75: pep->slambda = 1.0;
76: pep->refine = PEP_REFINE_NONE;
77: pep->npart = 1;
78: pep->rtol = PETSC_DEFAULT;
79: pep->rits = PETSC_DEFAULT;
80: pep->scheme = (PEPRefineScheme)0;
81: pep->extract = (PEPExtract)0;
82: pep->trackall = PETSC_FALSE;
84: pep->converged = PEPConvergedRelative;
85: pep->convergeduser = NULL;
86: pep->convergeddestroy= NULL;
87: pep->stopping = PEPStoppingBasic;
88: pep->stoppinguser = NULL;
89: pep->stoppingdestroy = NULL;
90: pep->convergedctx = NULL;
91: pep->stoppingctx = NULL;
92: pep->numbermonitors = 0;
94: pep->st = NULL;
95: pep->ds = NULL;
96: pep->V = NULL;
97: pep->rg = NULL;
98: pep->A = NULL;
99: pep->nmat = 0;
100: pep->Dl = NULL;
101: pep->Dr = NULL;
102: pep->IS = NULL;
103: pep->eigr = NULL;
104: pep->eigi = NULL;
105: pep->errest = NULL;
106: pep->perm = NULL;
107: pep->pbc = NULL;
108: pep->solvematcoeffs = NULL;
109: pep->nwork = 0;
110: pep->work = NULL;
111: pep->refineksp = NULL;
112: pep->refinesubc = NULL;
113: pep->data = NULL;
115: pep->state = PEP_STATE_INITIAL;
116: pep->nconv = 0;
117: pep->its = 0;
118: pep->n = 0;
119: pep->nloc = 0;
120: pep->nrma = NULL;
121: pep->sfactor_set = PETSC_FALSE;
122: pep->lineariz = PETSC_FALSE;
123: pep->reason = PEP_CONVERGED_ITERATING;
125: PetscNewLog(pep,&pep->sc);
126: *outpep = pep;
127: return(0);
128: }
130: /*@C
131: PEPSetType - Selects the particular solver to be used in the PEP object.
133: Logically Collective on pep
135: Input Parameters:
136: + pep - the polynomial eigensolver context
137: - type - a known method
139: Options Database Key:
140: . -pep_type <method> - Sets the method; use -help for a list
141: of available methods
143: Notes:
144: See "slepc/include/slepcpep.h" for available methods. The default
145: is PEPTOAR.
147: Normally, it is best to use the PEPSetFromOptions() command and
148: then set the PEP type from the options database rather than by using
149: this routine. Using the options database provides the user with
150: maximum flexibility in evaluating the different available methods.
151: The PEPSetType() routine is provided for those situations where it
152: is necessary to set the iterative solver independently of the command
153: line or options database.
155: Level: intermediate
157: .seealso: PEPType
158: @*/
159: PetscErrorCode PEPSetType(PEP pep,PEPType type)
160: {
161: PetscErrorCode ierr,(*r)(PEP);
162: PetscBool match;
168: PetscObjectTypeCompare((PetscObject)pep,type,&match);
169: if (match) return(0);
171: PetscFunctionListFind(PEPList,type,&r);
172: if (!r) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown PEP type given: %s",type);
174: if (pep->ops->destroy) { (*pep->ops->destroy)(pep); }
175: PetscMemzero(pep->ops,sizeof(struct _PEPOps));
177: pep->state = PEP_STATE_INITIAL;
178: PetscObjectChangeTypeName((PetscObject)pep,type);
179: (*r)(pep);
180: return(0);
181: }
183: /*@C
184: PEPGetType - Gets the PEP type as a string from the PEP object.
186: Not Collective
188: Input Parameter:
189: . pep - the eigensolver context
191: Output Parameter:
192: . name - name of PEP method
194: Level: intermediate
196: .seealso: PEPSetType()
197: @*/
198: PetscErrorCode PEPGetType(PEP pep,PEPType *type)
199: {
203: *type = ((PetscObject)pep)->type_name;
204: return(0);
205: }
207: /*@C
208: PEPRegister - Adds a method to the polynomial eigenproblem solver package.
210: Not Collective
212: Input Parameters:
213: + name - name of a new user-defined solver
214: - function - routine to create the solver context
216: Notes:
217: PEPRegister() may be called multiple times to add several user-defined solvers.
219: Sample usage:
220: .vb
221: PEPRegister("my_solver",MySolverCreate);
222: .ve
224: Then, your solver can be chosen with the procedural interface via
225: $ PEPSetType(pep,"my_solver")
226: or at runtime via the option
227: $ -pep_type my_solver
229: Level: advanced
231: .seealso: PEPRegisterAll()
232: @*/
233: PetscErrorCode PEPRegister(const char *name,PetscErrorCode (*function)(PEP))
234: {
238: PEPInitializePackage();
239: PetscFunctionListAdd(&PEPList,name,function);
240: return(0);
241: }
243: /*@C
244: PEPMonitorRegister - Adds PEP monitor routine.
246: Not Collective
248: Input Parameters:
249: + name - name of a new monitor routine
250: . vtype - a PetscViewerType for the output
251: . format - a PetscViewerFormat for the output
252: . monitor - monitor routine
253: . create - creation routine, or NULL
254: - destroy - destruction routine, or NULL
256: Notes:
257: PEPMonitorRegister() may be called multiple times to add several user-defined monitors.
259: Sample usage:
260: .vb
261: PEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
262: .ve
264: Then, your monitor can be chosen with the procedural interface via
265: $ PEPMonitorSetFromOptions(pep,"-pep_monitor_my_monitor","my_monitor",NULL)
266: or at runtime via the option
267: $ -pep_monitor_my_monitor
269: Level: advanced
271: .seealso: PEPMonitorRegisterAll()
272: @*/
273: 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**))
274: {
275: char key[PETSC_MAX_PATH_LEN];
279: PEPInitializePackage();
280: SlepcMonitorMakeKey_Internal(name,vtype,format,key);
281: PetscFunctionListAdd(&PEPMonitorList,key,monitor);
282: if (create) { PetscFunctionListAdd(&PEPMonitorCreateList,key,create); }
283: if (destroy) { PetscFunctionListAdd(&PEPMonitorDestroyList,key,destroy); }
284: return(0);
285: }
287: /*@
288: PEPReset - Resets the PEP context to the initial state (prior to setup)
289: and destroys any allocated Vecs and Mats.
291: Collective on pep
293: Input Parameter:
294: . pep - eigensolver context obtained from PEPCreate()
296: Level: advanced
298: .seealso: PEPDestroy()
299: @*/
300: PetscErrorCode PEPReset(PEP pep)
301: {
306: if (!pep) return(0);
307: if (pep->ops->reset) { (pep->ops->reset)(pep); }
308: if (pep->st) { STReset(pep->st); }
309: if (pep->refineksp) { KSPReset(pep->refineksp); }
310: if (pep->nmat) {
311: MatDestroyMatrices(pep->nmat,&pep->A);
312: PetscFree2(pep->pbc,pep->nrma);
313: PetscFree(pep->solvematcoeffs);
314: pep->nmat = 0;
315: }
316: VecDestroy(&pep->Dl);
317: VecDestroy(&pep->Dr);
318: BVDestroy(&pep->V);
319: VecDestroyVecs(pep->nwork,&pep->work);
320: pep->nwork = 0;
321: pep->state = PEP_STATE_INITIAL;
322: return(0);
323: }
325: /*@C
326: PEPDestroy - Destroys the PEP context.
328: Collective on pep
330: Input Parameter:
331: . pep - eigensolver context obtained from PEPCreate()
333: Level: beginner
335: .seealso: PEPCreate(), PEPSetUp(), PEPSolve()
336: @*/
337: PetscErrorCode PEPDestroy(PEP *pep)
338: {
342: if (!*pep) return(0);
344: if (--((PetscObject)(*pep))->refct > 0) { *pep = 0; return(0); }
345: PEPReset(*pep);
346: if ((*pep)->ops->destroy) { (*(*pep)->ops->destroy)(*pep); }
347: if ((*pep)->eigr) {
348: PetscFree4((*pep)->eigr,(*pep)->eigi,(*pep)->errest,(*pep)->perm);
349: }
350: STDestroy(&(*pep)->st);
351: RGDestroy(&(*pep)->rg);
352: DSDestroy(&(*pep)->ds);
353: KSPDestroy(&(*pep)->refineksp);
354: PetscSubcommDestroy(&(*pep)->refinesubc);
355: PetscFree((*pep)->sc);
356: /* just in case the initial vectors have not been used */
357: SlepcBasisDestroy_Private(&(*pep)->nini,&(*pep)->IS);
358: if ((*pep)->convergeddestroy) {
359: (*(*pep)->convergeddestroy)((*pep)->convergedctx);
360: }
361: PEPMonitorCancel(*pep);
362: PetscHeaderDestroy(pep);
363: return(0);
364: }
366: /*@
367: PEPSetBV - Associates a basis vectors object to the polynomial eigensolver.
369: Collective on pep
371: Input Parameters:
372: + pep - eigensolver context obtained from PEPCreate()
373: - bv - the basis vectors object
375: Note:
376: Use PEPGetBV() to retrieve the basis vectors context (for example,
377: to free it at the end of the computations).
379: Level: advanced
381: .seealso: PEPGetBV()
382: @*/
383: PetscErrorCode PEPSetBV(PEP pep,BV bv)
384: {
391: PetscObjectReference((PetscObject)bv);
392: BVDestroy(&pep->V);
393: pep->V = bv;
394: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->V);
395: return(0);
396: }
398: /*@
399: PEPGetBV - Obtain the basis vectors object associated to the polynomial
400: eigensolver object.
402: Not Collective
404: Input Parameters:
405: . pep - eigensolver context obtained from PEPCreate()
407: Output Parameter:
408: . bv - basis vectors context
410: Level: advanced
412: .seealso: PEPSetBV()
413: @*/
414: PetscErrorCode PEPGetBV(PEP pep,BV *bv)
415: {
421: if (!pep->V) {
422: BVCreate(PetscObjectComm((PetscObject)pep),&pep->V);
423: PetscObjectIncrementTabLevel((PetscObject)pep->V,(PetscObject)pep,0);
424: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->V);
425: PetscObjectSetOptions((PetscObject)pep->V,((PetscObject)pep)->options);
426: }
427: *bv = pep->V;
428: return(0);
429: }
431: /*@
432: PEPSetRG - Associates a region object to the polynomial eigensolver.
434: Collective on pep
436: Input Parameters:
437: + pep - eigensolver context obtained from PEPCreate()
438: - rg - the region object
440: Note:
441: Use PEPGetRG() to retrieve the region context (for example,
442: to free it at the end of the computations).
444: Level: advanced
446: .seealso: PEPGetRG()
447: @*/
448: PetscErrorCode PEPSetRG(PEP pep,RG rg)
449: {
454: if (rg) {
457: }
458: PetscObjectReference((PetscObject)rg);
459: RGDestroy(&pep->rg);
460: pep->rg = rg;
461: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
462: return(0);
463: }
465: /*@
466: PEPGetRG - Obtain the region object associated to the
467: polynomial eigensolver object.
469: Not Collective
471: Input Parameters:
472: . pep - eigensolver context obtained from PEPCreate()
474: Output Parameter:
475: . rg - region context
477: Level: advanced
479: .seealso: PEPSetRG()
480: @*/
481: PetscErrorCode PEPGetRG(PEP pep,RG *rg)
482: {
488: if (!pep->rg) {
489: RGCreate(PetscObjectComm((PetscObject)pep),&pep->rg);
490: PetscObjectIncrementTabLevel((PetscObject)pep->rg,(PetscObject)pep,0);
491: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
492: PetscObjectSetOptions((PetscObject)pep->rg,((PetscObject)pep)->options);
493: }
494: *rg = pep->rg;
495: return(0);
496: }
498: /*@
499: PEPSetDS - Associates a direct solver object to the polynomial eigensolver.
501: Collective on pep
503: Input Parameters:
504: + pep - eigensolver context obtained from PEPCreate()
505: - ds - the direct solver object
507: Note:
508: Use PEPGetDS() to retrieve the direct solver context (for example,
509: to free it at the end of the computations).
511: Level: advanced
513: .seealso: PEPGetDS()
514: @*/
515: PetscErrorCode PEPSetDS(PEP pep,DS ds)
516: {
523: PetscObjectReference((PetscObject)ds);
524: DSDestroy(&pep->ds);
525: pep->ds = ds;
526: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
527: return(0);
528: }
530: /*@
531: PEPGetDS - Obtain the direct solver object associated to the
532: polynomial eigensolver object.
534: Not Collective
536: Input Parameters:
537: . pep - eigensolver context obtained from PEPCreate()
539: Output Parameter:
540: . ds - direct solver context
542: Level: advanced
544: .seealso: PEPSetDS()
545: @*/
546: PetscErrorCode PEPGetDS(PEP pep,DS *ds)
547: {
553: if (!pep->ds) {
554: DSCreate(PetscObjectComm((PetscObject)pep),&pep->ds);
555: PetscObjectIncrementTabLevel((PetscObject)pep->ds,(PetscObject)pep,0);
556: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
557: PetscObjectSetOptions((PetscObject)pep->ds,((PetscObject)pep)->options);
558: }
559: *ds = pep->ds;
560: return(0);
561: }
563: /*@
564: PEPSetST - Associates a spectral transformation object to the eigensolver.
566: Collective on pep
568: Input Parameters:
569: + pep - eigensolver context obtained from PEPCreate()
570: - st - the spectral transformation object
572: Note:
573: Use PEPGetST() to retrieve the spectral transformation context (for example,
574: to free it at the end of the computations).
576: Level: advanced
578: .seealso: PEPGetST()
579: @*/
580: PetscErrorCode PEPSetST(PEP pep,ST st)
581: {
588: PetscObjectReference((PetscObject)st);
589: STDestroy(&pep->st);
590: pep->st = st;
591: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
592: return(0);
593: }
595: /*@
596: PEPGetST - Obtain the spectral transformation (ST) object associated
597: to the eigensolver object.
599: Not Collective
601: Input Parameters:
602: . pep - eigensolver context obtained from PEPCreate()
604: Output Parameter:
605: . st - spectral transformation context
607: Level: intermediate
609: .seealso: PEPSetST()
610: @*/
611: PetscErrorCode PEPGetST(PEP pep,ST *st)
612: {
618: if (!pep->st) {
619: STCreate(PetscObjectComm((PetscObject)pep),&pep->st);
620: PetscObjectIncrementTabLevel((PetscObject)pep->st,(PetscObject)pep,0);
621: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
622: PetscObjectSetOptions((PetscObject)pep->st,((PetscObject)pep)->options);
623: }
624: *st = pep->st;
625: return(0);
626: }
628: /*@
629: PEPRefineGetKSP - Obtain the ksp object used by the eigensolver
630: object in the refinement phase.
632: Not Collective
634: Input Parameters:
635: . pep - eigensolver context obtained from PEPCreate()
637: Output Parameter:
638: . ksp - ksp context
640: Level: advanced
642: .seealso: PEPSetRefine()
643: @*/
644: PetscErrorCode PEPRefineGetKSP(PEP pep,KSP *ksp)
645: {
651: if (!pep->refineksp) {
652: if (pep->npart>1) {
653: /* Split in subcomunicators */
654: PetscSubcommCreate(PetscObjectComm((PetscObject)pep),&pep->refinesubc);
655: PetscSubcommSetNumber(pep->refinesubc,pep->npart);
656: PetscSubcommSetType(pep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
657: PetscLogObjectMemory((PetscObject)pep,sizeof(PetscSubcomm));
658: }
659: KSPCreate((pep->npart==1)?PetscObjectComm((PetscObject)pep):PetscSubcommChild(pep->refinesubc),&pep->refineksp);
660: PetscObjectIncrementTabLevel((PetscObject)pep->refineksp,(PetscObject)pep,0);
661: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->refineksp);
662: PetscObjectSetOptions((PetscObject)pep->refineksp,((PetscObject)pep)->options);
663: KSPSetOptionsPrefix(*ksp,((PetscObject)pep)->prefix);
664: KSPAppendOptionsPrefix(*ksp,"pep_refine_");
665: KSPSetTolerances(pep->refineksp,SlepcDefaultTol(pep->rtol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
666: }
667: *ksp = pep->refineksp;
668: return(0);
669: }
671: /*@
672: PEPSetTarget - Sets the value of the target.
674: Logically Collective on pep
676: Input Parameters:
677: + pep - eigensolver context
678: - target - the value of the target
680: Options Database Key:
681: . -pep_target <scalar> - the value of the target
683: Notes:
684: The target is a scalar value used to determine the portion of the spectrum
685: of interest. It is used in combination with PEPSetWhichEigenpairs().
687: In the case of complex scalars, a complex value can be provided in the
688: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
689: -pep_target 1.0+2.0i
691: Level: intermediate
693: .seealso: PEPGetTarget(), PEPSetWhichEigenpairs()
694: @*/
695: PetscErrorCode PEPSetTarget(PEP pep,PetscScalar target)
696: {
702: pep->target = target;
703: if (!pep->st) { PEPGetST(pep,&pep->st); }
704: STSetDefaultShift(pep->st,target);
705: return(0);
706: }
708: /*@
709: PEPGetTarget - Gets the value of the target.
711: Not Collective
713: Input Parameter:
714: . pep - eigensolver context
716: Output Parameter:
717: . target - the value of the target
719: Note:
720: If the target was not set by the user, then zero is returned.
722: Level: intermediate
724: .seealso: PEPSetTarget()
725: @*/
726: PetscErrorCode PEPGetTarget(PEP pep,PetscScalar* target)
727: {
731: *target = pep->target;
732: return(0);
733: }
735: /*@
736: PEPSetInterval - Defines the computational interval for spectrum slicing.
738: Logically Collective on pep
740: Input Parameters:
741: + pep - eigensolver context
742: . inta - left end of the interval
743: - intb - right end of the interval
745: Options Database Key:
746: . -pep_interval <a,b> - set [a,b] as the interval of interest
748: Notes:
749: Spectrum slicing is a technique employed for computing all eigenvalues of
750: symmetric eigenproblems in a given interval. This function provides the
751: interval to be considered. It must be used in combination with PEP_ALL, see
752: PEPSetWhichEigenpairs().
754: In the command-line option, two values must be provided. For an open interval,
755: one can give an infinite, e.g., -pep_interval 1.0,inf or -pep_interval -inf,1.0.
756: An open interval in the programmatic interface can be specified with
757: PETSC_MAX_REAL and -PETSC_MAX_REAL.
759: Level: intermediate
761: .seealso: PEPGetInterval(), PEPSetWhichEigenpairs()
762: @*/
763: PetscErrorCode PEPSetInterval(PEP pep,PetscReal inta,PetscReal intb)
764: {
769: if (inta>=intb) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
770: if (pep->inta != inta || pep->intb != intb) {
771: pep->inta = inta;
772: pep->intb = intb;
773: pep->state = PEP_STATE_INITIAL;
774: }
775: return(0);
776: }
778: /*@
779: PEPGetInterval - Gets the computational interval for spectrum slicing.
781: Not Collective
783: Input Parameter:
784: . pep - eigensolver context
786: Output Parameters:
787: + inta - left end of the interval
788: - intb - right end of the interval
790: Level: intermediate
792: Note:
793: If the interval was not set by the user, then zeros are returned.
795: .seealso: PEPSetInterval()
796: @*/
797: PetscErrorCode PEPGetInterval(PEP pep,PetscReal* inta,PetscReal* intb)
798: {
801: if (inta) *inta = pep->inta;
802: if (intb) *intb = pep->intb;
803: return(0);
804: }