Actual source code: epsbasic.c
slepc-3.15.0 2021-03-31
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 EPS routines
12: */
14: #include <slepc/private/epsimpl.h>
16: /* Logging support */
17: PetscClassId EPS_CLASSID = 0;
18: PetscLogEvent EPS_SetUp = 0,EPS_Solve = 0;
20: /* List of registered EPS routines */
21: PetscFunctionList EPSList = 0;
22: PetscBool EPSRegisterAllCalled = PETSC_FALSE;
24: /* List of registered EPS monitors */
25: PetscFunctionList EPSMonitorList = NULL;
26: PetscFunctionList EPSMonitorCreateList = NULL;
27: PetscFunctionList EPSMonitorDestroyList = NULL;
28: PetscBool EPSMonitorRegisterAllCalled = PETSC_FALSE;
30: /*@
31: EPSCreate - Creates the default EPS context.
33: Collective
35: Input Parameter:
36: . comm - MPI communicator
38: Output Parameter:
39: . eps - location to put the EPS context
41: Note:
42: The default EPS type is EPSKRYLOVSCHUR
44: Level: beginner
46: .seealso: EPSSetUp(), EPSSolve(), EPSDestroy(), EPS
47: @*/
48: PetscErrorCode EPSCreate(MPI_Comm comm,EPS *outeps)
49: {
51: EPS eps;
55: *outeps = 0;
56: EPSInitializePackage();
57: SlepcHeaderCreate(eps,EPS_CLASSID,"EPS","Eigenvalue Problem Solver","EPS",comm,EPSDestroy,EPSView);
59: eps->max_it = PETSC_DEFAULT;
60: eps->nev = 1;
61: eps->ncv = PETSC_DEFAULT;
62: eps->mpd = PETSC_DEFAULT;
63: eps->nini = 0;
64: eps->nds = 0;
65: eps->target = 0.0;
66: eps->tol = PETSC_DEFAULT;
67: eps->conv = EPS_CONV_REL;
68: eps->stop = EPS_STOP_BASIC;
69: eps->which = (EPSWhich)0;
70: eps->inta = 0.0;
71: eps->intb = 0.0;
72: eps->problem_type = (EPSProblemType)0;
73: eps->extraction = EPS_RITZ;
74: eps->balance = EPS_BALANCE_NONE;
75: eps->balance_its = 5;
76: eps->balance_cutoff = 1e-8;
77: eps->trueres = PETSC_FALSE;
78: eps->trackall = PETSC_FALSE;
79: eps->purify = PETSC_TRUE;
80: eps->twosided = PETSC_FALSE;
82: eps->converged = EPSConvergedRelative;
83: eps->convergeduser = NULL;
84: eps->convergeddestroy= NULL;
85: eps->stopping = EPSStoppingBasic;
86: eps->stoppinguser = NULL;
87: eps->stoppingdestroy = NULL;
88: eps->arbitrary = NULL;
89: eps->convergedctx = NULL;
90: eps->stoppingctx = NULL;
91: eps->arbitraryctx = NULL;
92: eps->numbermonitors = 0;
94: eps->st = NULL;
95: eps->ds = NULL;
96: eps->V = NULL;
97: eps->W = NULL;
98: eps->rg = NULL;
99: eps->D = NULL;
100: eps->IS = NULL;
101: eps->ISL = NULL;
102: eps->defl = NULL;
103: eps->eigr = NULL;
104: eps->eigi = NULL;
105: eps->errest = NULL;
106: eps->rr = NULL;
107: eps->ri = NULL;
108: eps->perm = NULL;
109: eps->nwork = 0;
110: eps->work = NULL;
111: eps->data = NULL;
113: eps->state = EPS_STATE_INITIAL;
114: eps->categ = EPS_CATEGORY_KRYLOV;
115: eps->nconv = 0;
116: eps->its = 0;
117: eps->nloc = 0;
118: eps->nrma = 0.0;
119: eps->nrmb = 0.0;
120: eps->useds = PETSC_FALSE;
121: eps->isgeneralized = PETSC_FALSE;
122: eps->ispositive = PETSC_FALSE;
123: eps->ishermitian = PETSC_FALSE;
124: eps->reason = EPS_CONVERGED_ITERATING;
126: PetscNewLog(eps,&eps->sc);
127: *outeps = eps;
128: return(0);
129: }
131: /*@C
132: EPSSetType - Selects the particular solver to be used in the EPS object.
134: Logically Collective on eps
136: Input Parameters:
137: + eps - the eigensolver context
138: - type - a known method
140: Options Database Key:
141: . -eps_type <method> - Sets the method; use -help for a list
142: of available methods
144: Notes:
145: See "slepc/include/slepceps.h" for available methods. The default
146: is EPSKRYLOVSCHUR.
148: Normally, it is best to use the EPSSetFromOptions() command and
149: then set the EPS type from the options database rather than by using
150: this routine. Using the options database provides the user with
151: maximum flexibility in evaluating the different available methods.
152: The EPSSetType() routine is provided for those situations where it
153: is necessary to set the iterative solver independently of the command
154: line or options database.
156: Level: intermediate
158: .seealso: STSetType(), EPSType
159: @*/
160: PetscErrorCode EPSSetType(EPS eps,EPSType type)
161: {
162: PetscErrorCode ierr,(*r)(EPS);
163: PetscBool match;
169: PetscObjectTypeCompare((PetscObject)eps,type,&match);
170: if (match) return(0);
172: PetscFunctionListFind(EPSList,type,&r);
173: if (!r) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown EPS type given: %s",type);
175: if (eps->ops->destroy) { (*eps->ops->destroy)(eps); }
176: PetscMemzero(eps->ops,sizeof(struct _EPSOps));
178: eps->state = EPS_STATE_INITIAL;
179: PetscObjectChangeTypeName((PetscObject)eps,type);
180: (*r)(eps);
181: return(0);
182: }
184: /*@C
185: EPSGetType - Gets the EPS type as a string from the EPS object.
187: Not Collective
189: Input Parameter:
190: . eps - the eigensolver context
192: Output Parameter:
193: . name - name of EPS method
195: Level: intermediate
197: .seealso: EPSSetType()
198: @*/
199: PetscErrorCode EPSGetType(EPS eps,EPSType *type)
200: {
204: *type = ((PetscObject)eps)->type_name;
205: return(0);
206: }
208: /*@C
209: EPSRegister - Adds a method to the eigenproblem solver package.
211: Not Collective
213: Input Parameters:
214: + name - name of a new user-defined solver
215: - function - routine to create the solver context
217: Notes:
218: EPSRegister() may be called multiple times to add several user-defined solvers.
220: Sample usage:
221: .vb
222: EPSRegister("my_solver",MySolverCreate);
223: .ve
225: Then, your solver can be chosen with the procedural interface via
226: $ EPSSetType(eps,"my_solver")
227: or at runtime via the option
228: $ -eps_type my_solver
230: Level: advanced
232: .seealso: EPSRegisterAll()
233: @*/
234: PetscErrorCode EPSRegister(const char *name,PetscErrorCode (*function)(EPS))
235: {
239: EPSInitializePackage();
240: PetscFunctionListAdd(&EPSList,name,function);
241: return(0);
242: }
244: /*@C
245: EPSMonitorRegister - Adds EPS monitor routine.
247: Not Collective
249: Input Parameters:
250: + name - name of a new monitor routine
251: . vtype - a PetscViewerType for the output
252: . format - a PetscViewerFormat for the output
253: . monitor - monitor routine
254: . create - creation routine, or NULL
255: - destroy - destruction routine, or NULL
257: Notes:
258: EPSMonitorRegister() may be called multiple times to add several user-defined monitors.
260: Sample usage:
261: .vb
262: EPSMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
263: .ve
265: Then, your monitor can be chosen with the procedural interface via
266: $ EPSMonitorSetFromOptions(eps,"-eps_monitor_my_monitor","my_monitor",NULL)
267: or at runtime via the option
268: $ -eps_monitor_my_monitor
270: Level: advanced
272: .seealso: EPSMonitorRegisterAll()
273: @*/
274: PetscErrorCode EPSMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(EPS,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscReal*,PetscInt,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
275: {
276: char key[PETSC_MAX_PATH_LEN];
280: EPSInitializePackage();
281: SlepcMonitorMakeKey_Internal(name,vtype,format,key);
282: PetscFunctionListAdd(&EPSMonitorList,key,monitor);
283: if (create) { PetscFunctionListAdd(&EPSMonitorCreateList,key,create); }
284: if (destroy) { PetscFunctionListAdd(&EPSMonitorDestroyList,key,destroy); }
285: return(0);
286: }
288: /*@
289: EPSReset - Resets the EPS context to the initial state (prior to setup)
290: and destroys any allocated Vecs and Mats.
292: Collective on eps
294: Input Parameter:
295: . eps - eigensolver context obtained from EPSCreate()
297: Note:
298: This can be used when a problem of different matrix size wants to be solved.
299: All options that have previously been set are preserved, so in a next use
300: the solver configuration is the same, but new sizes for matrices and vectors
301: are allowed.
303: Level: advanced
305: .seealso: EPSDestroy()
306: @*/
307: PetscErrorCode EPSReset(EPS eps)
308: {
313: if (!eps) return(0);
314: if (eps->ops->reset) { (eps->ops->reset)(eps); }
315: if (eps->st) { STReset(eps->st); }
316: VecDestroy(&eps->D);
317: BVDestroy(&eps->V);
318: BVDestroy(&eps->W);
319: VecDestroyVecs(eps->nwork,&eps->work);
320: eps->nwork = 0;
321: eps->state = EPS_STATE_INITIAL;
322: return(0);
323: }
325: /*@C
326: EPSDestroy - Destroys the EPS context.
328: Collective on eps
330: Input Parameter:
331: . eps - eigensolver context obtained from EPSCreate()
333: Level: beginner
335: .seealso: EPSCreate(), EPSSetUp(), EPSSolve()
336: @*/
337: PetscErrorCode EPSDestroy(EPS *eps)
338: {
342: if (!*eps) return(0);
344: if (--((PetscObject)(*eps))->refct > 0) { *eps = 0; return(0); }
345: EPSReset(*eps);
346: if ((*eps)->ops->destroy) { (*(*eps)->ops->destroy)(*eps); }
347: if ((*eps)->eigr) {
348: PetscFree4((*eps)->eigr,(*eps)->eigi,(*eps)->errest,(*eps)->perm);
349: }
350: if ((*eps)->rr) {
351: PetscFree2((*eps)->rr,(*eps)->ri);
352: }
353: STDestroy(&(*eps)->st);
354: RGDestroy(&(*eps)->rg);
355: DSDestroy(&(*eps)->ds);
356: PetscFree((*eps)->sc);
357: /* just in case the initial vectors have not been used */
358: SlepcBasisDestroy_Private(&(*eps)->nds,&(*eps)->defl);
359: SlepcBasisDestroy_Private(&(*eps)->nini,&(*eps)->IS);
360: SlepcBasisDestroy_Private(&(*eps)->ninil,&(*eps)->ISL);
361: if ((*eps)->convergeddestroy) {
362: (*(*eps)->convergeddestroy)((*eps)->convergedctx);
363: }
364: EPSMonitorCancel(*eps);
365: PetscHeaderDestroy(eps);
366: return(0);
367: }
369: /*@
370: EPSSetTarget - Sets the value of the target.
372: Logically Collective on eps
374: Input Parameters:
375: + eps - eigensolver context
376: - target - the value of the target
378: Options Database Key:
379: . -eps_target <scalar> - the value of the target
381: Notes:
382: The target is a scalar value used to determine the portion of the spectrum
383: of interest. It is used in combination with EPSSetWhichEigenpairs().
385: In the case of complex scalars, a complex value can be provided in the
386: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
387: -eps_target 1.0+2.0i
389: Level: intermediate
391: .seealso: EPSGetTarget(), EPSSetWhichEigenpairs()
392: @*/
393: PetscErrorCode EPSSetTarget(EPS eps,PetscScalar target)
394: {
400: eps->target = target;
401: if (!eps->st) { EPSGetST(eps,&eps->st); }
402: STSetDefaultShift(eps->st,target);
403: return(0);
404: }
406: /*@
407: EPSGetTarget - Gets the value of the target.
409: Not Collective
411: Input Parameter:
412: . eps - eigensolver context
414: Output Parameter:
415: . target - the value of the target
417: Note:
418: If the target was not set by the user, then zero is returned.
420: Level: intermediate
422: .seealso: EPSSetTarget()
423: @*/
424: PetscErrorCode EPSGetTarget(EPS eps,PetscScalar* target)
425: {
429: *target = eps->target;
430: return(0);
431: }
433: /*@
434: EPSSetInterval - Defines the computational interval for spectrum slicing.
436: Logically Collective on eps
438: Input Parameters:
439: + eps - eigensolver context
440: . inta - left end of the interval
441: - intb - right end of the interval
443: Options Database Key:
444: . -eps_interval <a,b> - set [a,b] as the interval of interest
446: Notes:
447: Spectrum slicing is a technique employed for computing all eigenvalues of
448: symmetric eigenproblems in a given interval. This function provides the
449: interval to be considered. It must be used in combination with EPS_ALL, see
450: EPSSetWhichEigenpairs().
452: In the command-line option, two values must be provided. For an open interval,
453: one can give an infinite, e.g., -eps_interval 1.0,inf or -eps_interval -inf,1.0.
454: An open interval in the programmatic interface can be specified with
455: PETSC_MAX_REAL and -PETSC_MAX_REAL.
457: Level: intermediate
459: .seealso: EPSGetInterval(), EPSSetWhichEigenpairs()
460: @*/
461: PetscErrorCode EPSSetInterval(EPS eps,PetscReal inta,PetscReal intb)
462: {
467: if (inta>=intb) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
468: if (eps->inta != inta || eps->intb != intb) {
469: eps->inta = inta;
470: eps->intb = intb;
471: eps->state = EPS_STATE_INITIAL;
472: }
473: return(0);
474: }
476: /*@
477: EPSGetInterval - Gets the computational interval for spectrum slicing.
479: Not Collective
481: Input Parameter:
482: . eps - eigensolver context
484: Output Parameters:
485: + inta - left end of the interval
486: - intb - right end of the interval
488: Level: intermediate
490: Note:
491: If the interval was not set by the user, then zeros are returned.
493: .seealso: EPSSetInterval()
494: @*/
495: PetscErrorCode EPSGetInterval(EPS eps,PetscReal* inta,PetscReal* intb)
496: {
499: if (inta) *inta = eps->inta;
500: if (intb) *intb = eps->intb;
501: return(0);
502: }
504: /*@
505: EPSSetST - Associates a spectral transformation object to the eigensolver.
507: Collective on eps
509: Input Parameters:
510: + eps - eigensolver context obtained from EPSCreate()
511: - st - the spectral transformation object
513: Note:
514: Use EPSGetST() to retrieve the spectral transformation context (for example,
515: to free it at the end of the computations).
517: Level: advanced
519: .seealso: EPSGetST()
520: @*/
521: PetscErrorCode EPSSetST(EPS eps,ST st)
522: {
529: PetscObjectReference((PetscObject)st);
530: STDestroy(&eps->st);
531: eps->st = st;
532: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->st);
533: return(0);
534: }
536: /*@
537: EPSGetST - Obtain the spectral transformation (ST) object associated
538: to the eigensolver object.
540: Not Collective
542: Input Parameters:
543: . eps - eigensolver context obtained from EPSCreate()
545: Output Parameter:
546: . st - spectral transformation context
548: Level: intermediate
550: .seealso: EPSSetST()
551: @*/
552: PetscErrorCode EPSGetST(EPS eps,ST *st)
553: {
559: if (!eps->st) {
560: STCreate(PetscObjectComm((PetscObject)eps),&eps->st);
561: PetscObjectIncrementTabLevel((PetscObject)eps->st,(PetscObject)eps,0);
562: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->st);
563: PetscObjectSetOptions((PetscObject)eps->st,((PetscObject)eps)->options);
564: }
565: *st = eps->st;
566: return(0);
567: }
569: /*@
570: EPSSetBV - Associates a basis vectors object to the eigensolver.
572: Collective on eps
574: Input Parameters:
575: + eps - eigensolver context obtained from EPSCreate()
576: - V - the basis vectors object
578: Level: advanced
580: .seealso: EPSGetBV()
581: @*/
582: PetscErrorCode EPSSetBV(EPS eps,BV V)
583: {
590: PetscObjectReference((PetscObject)V);
591: BVDestroy(&eps->V);
592: eps->V = V;
593: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->V);
594: return(0);
595: }
597: /*@
598: EPSGetBV - Obtain the basis vectors object associated to the eigensolver object.
600: Not Collective
602: Input Parameters:
603: . eps - eigensolver context obtained from EPSCreate()
605: Output Parameter:
606: . V - basis vectors context
608: Level: advanced
610: .seealso: EPSSetBV()
611: @*/
612: PetscErrorCode EPSGetBV(EPS eps,BV *V)
613: {
619: if (!eps->V) {
620: BVCreate(PetscObjectComm((PetscObject)eps),&eps->V);
621: PetscObjectIncrementTabLevel((PetscObject)eps->V,(PetscObject)eps,0);
622: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->V);
623: PetscObjectSetOptions((PetscObject)eps->V,((PetscObject)eps)->options);
624: }
625: *V = eps->V;
626: return(0);
627: }
629: /*@
630: EPSSetRG - Associates a region object to the eigensolver.
632: Collective on eps
634: Input Parameters:
635: + eps - eigensolver context obtained from EPSCreate()
636: - rg - the region object
638: Note:
639: Use EPSGetRG() to retrieve the region context (for example,
640: to free it at the end of the computations).
642: Level: advanced
644: .seealso: EPSGetRG()
645: @*/
646: PetscErrorCode EPSSetRG(EPS eps,RG rg)
647: {
654: PetscObjectReference((PetscObject)rg);
655: RGDestroy(&eps->rg);
656: eps->rg = rg;
657: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->rg);
658: return(0);
659: }
661: /*@
662: EPSGetRG - Obtain the region object associated to the eigensolver.
664: Not Collective
666: Input Parameters:
667: . eps - eigensolver context obtained from EPSCreate()
669: Output Parameter:
670: . rg - region context
672: Level: advanced
674: .seealso: EPSSetRG()
675: @*/
676: PetscErrorCode EPSGetRG(EPS eps,RG *rg)
677: {
683: if (!eps->rg) {
684: RGCreate(PetscObjectComm((PetscObject)eps),&eps->rg);
685: PetscObjectIncrementTabLevel((PetscObject)eps->rg,(PetscObject)eps,0);
686: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->rg);
687: PetscObjectSetOptions((PetscObject)eps->rg,((PetscObject)eps)->options);
688: }
689: *rg = eps->rg;
690: return(0);
691: }
693: /*@
694: EPSSetDS - Associates a direct solver object to the eigensolver.
696: Collective on eps
698: Input Parameters:
699: + eps - eigensolver context obtained from EPSCreate()
700: - ds - the direct solver object
702: Note:
703: Use EPSGetDS() to retrieve the direct solver context (for example,
704: to free it at the end of the computations).
706: Level: advanced
708: .seealso: EPSGetDS()
709: @*/
710: PetscErrorCode EPSSetDS(EPS eps,DS ds)
711: {
718: PetscObjectReference((PetscObject)ds);
719: DSDestroy(&eps->ds);
720: eps->ds = ds;
721: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->ds);
722: return(0);
723: }
725: /*@
726: EPSGetDS - Obtain the direct solver object associated to the eigensolver object.
728: Not Collective
730: Input Parameters:
731: . eps - eigensolver context obtained from EPSCreate()
733: Output Parameter:
734: . ds - direct solver context
736: Level: advanced
738: .seealso: EPSSetDS()
739: @*/
740: PetscErrorCode EPSGetDS(EPS eps,DS *ds)
741: {
747: if (!eps->ds) {
748: DSCreate(PetscObjectComm((PetscObject)eps),&eps->ds);
749: PetscObjectIncrementTabLevel((PetscObject)eps->ds,(PetscObject)eps,0);
750: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->ds);
751: PetscObjectSetOptions((PetscObject)eps->ds,((PetscObject)eps)->options);
752: }
753: *ds = eps->ds;
754: return(0);
755: }
757: /*@
758: EPSIsGeneralized - Ask if the EPS object corresponds to a generalized
759: eigenvalue problem.
761: Not collective
763: Input Parameter:
764: . eps - the eigenproblem solver context
766: Output Parameter:
767: . is - the answer
769: Level: intermediate
771: .seealso: EPSIsHermitian(), EPSIsPositive()
772: @*/
773: PetscErrorCode EPSIsGeneralized(EPS eps,PetscBool* is)
774: {
778: *is = eps->isgeneralized;
779: return(0);
780: }
782: /*@
783: EPSIsHermitian - Ask if the EPS object corresponds to a Hermitian
784: eigenvalue problem.
786: Not collective
788: Input Parameter:
789: . eps - the eigenproblem solver context
791: Output Parameter:
792: . is - the answer
794: Level: intermediate
796: .seealso: EPSIsGeneralized(), EPSIsPositive()
797: @*/
798: PetscErrorCode EPSIsHermitian(EPS eps,PetscBool* is)
799: {
803: *is = eps->ishermitian;
804: return(0);
805: }
807: /*@
808: EPSIsPositive - Ask if the EPS object corresponds to an eigenvalue
809: problem type that requires a positive (semi-) definite matrix B.
811: Not collective
813: Input Parameter:
814: . eps - the eigenproblem solver context
816: Output Parameter:
817: . is - the answer
819: Level: intermediate
821: .seealso: EPSIsGeneralized(), EPSIsHermitian()
822: @*/
823: PetscErrorCode EPSIsPositive(EPS eps,PetscBool* is)
824: {
828: *is = eps->ispositive;
829: return(0);
830: }