Actual source code: epsbasic.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepceps.h" I*/
16: PetscFunctionList EPSList = 0;
17: PetscBool EPSRegisterAllCalled = PETSC_FALSE;
18: PetscClassId EPS_CLASSID = 0;
19: PetscLogEvent EPS_SetUp = 0,EPS_Solve = 0;
21: /*@
22: EPSCreate - Creates the default EPS context.
24: Collective on MPI_Comm
26: Input Parameter:
27: . comm - MPI communicator
29: Output Parameter:
30: . eps - location to put the EPS context
32: Note:
33: The default EPS type is EPSKRYLOVSCHUR
35: Level: beginner
37: .seealso: EPSSetUp(), EPSSolve(), EPSDestroy(), EPS
38: @*/
39: PetscErrorCode EPSCreate(MPI_Comm comm,EPS *outeps)
40: {
42: EPS eps;
46: *outeps = 0;
47: EPSInitializePackage();
48: SlepcHeaderCreate(eps,EPS_CLASSID,"EPS","Eigenvalue Problem Solver","EPS",comm,EPSDestroy,EPSView);
50: eps->max_it = 0;
51: eps->nev = 1;
52: eps->ncv = 0;
53: eps->mpd = 0;
54: eps->nini = 0;
55: eps->nds = 0;
56: eps->target = 0.0;
57: eps->tol = PETSC_DEFAULT;
58: eps->conv = EPS_CONV_REL;
59: eps->stop = EPS_STOP_BASIC;
60: eps->which = (EPSWhich)0;
61: eps->inta = 0.0;
62: eps->intb = 0.0;
63: eps->problem_type = (EPSProblemType)0;
64: eps->extraction = EPS_RITZ;
65: eps->balance = EPS_BALANCE_NONE;
66: eps->balance_its = 5;
67: eps->balance_cutoff = 1e-8;
68: eps->trueres = PETSC_FALSE;
69: eps->trackall = PETSC_FALSE;
70: eps->purify = PETSC_TRUE;
71: eps->twosided = PETSC_FALSE;
73: eps->converged = EPSConvergedRelative;
74: eps->convergeduser = NULL;
75: eps->convergeddestroy= NULL;
76: eps->stopping = EPSStoppingBasic;
77: eps->stoppinguser = NULL;
78: eps->stoppingdestroy = NULL;
79: eps->arbitrary = NULL;
80: eps->convergedctx = NULL;
81: eps->stoppingctx = NULL;
82: eps->arbitraryctx = NULL;
83: eps->numbermonitors = 0;
85: eps->st = NULL;
86: eps->ds = NULL;
87: eps->dsts = NULL;
88: eps->V = NULL;
89: eps->W = NULL;
90: eps->rg = NULL;
91: eps->D = NULL;
92: eps->IS = NULL;
93: eps->defl = NULL;
94: eps->eigr = NULL;
95: eps->eigi = NULL;
96: eps->errest = NULL;
97: eps->rr = NULL;
98: eps->ri = NULL;
99: eps->perm = NULL;
100: eps->nwork = 0;
101: eps->work = NULL;
102: eps->data = NULL;
104: eps->state = EPS_STATE_INITIAL;
105: eps->categ = EPS_CATEGORY_KRYLOV;
106: eps->nconv = 0;
107: eps->its = 0;
108: eps->nloc = 0;
109: eps->nrma = 0.0;
110: eps->nrmb = 0.0;
111: eps->useds = PETSC_FALSE;
112: eps->hasts = PETSC_FALSE;
113: eps->isgeneralized = PETSC_FALSE;
114: eps->ispositive = PETSC_FALSE;
115: eps->ishermitian = PETSC_FALSE;
116: eps->reason = EPS_CONVERGED_ITERATING;
118: PetscNewLog(eps,&eps->sc);
119: *outeps = eps;
120: return(0);
121: }
123: /*@C
124: EPSSetType - Selects the particular solver to be used in the EPS object.
126: Logically Collective on EPS
128: Input Parameters:
129: + eps - the eigensolver context
130: - type - a known method
132: Options Database Key:
133: . -eps_type <method> - Sets the method; use -help for a list
134: of available methods
136: Notes:
137: See "slepc/include/slepceps.h" for available methods. The default
138: is EPSKRYLOVSCHUR.
140: Normally, it is best to use the EPSSetFromOptions() command and
141: then set the EPS type from the options database rather than by using
142: this routine. Using the options database provides the user with
143: maximum flexibility in evaluating the different available methods.
144: The EPSSetType() routine is provided for those situations where it
145: is necessary to set the iterative solver independently of the command
146: line or options database.
148: Level: intermediate
150: .seealso: STSetType(), EPSType
151: @*/
152: PetscErrorCode EPSSetType(EPS eps,EPSType type)
153: {
154: PetscErrorCode ierr,(*r)(EPS);
155: PetscBool match;
161: PetscObjectTypeCompare((PetscObject)eps,type,&match);
162: if (match) return(0);
164: PetscFunctionListFind(EPSList,type,&r);
165: if (!r) SETERRQ1(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown EPS type given: %s",type);
167: if (eps->ops->destroy) { (*eps->ops->destroy)(eps); }
168: PetscMemzero(eps->ops,sizeof(struct _EPSOps));
170: eps->state = EPS_STATE_INITIAL;
171: PetscObjectChangeTypeName((PetscObject)eps,type);
172: (*r)(eps);
173: return(0);
174: }
176: /*@C
177: EPSGetType - Gets the EPS type as a string from the EPS object.
179: Not Collective
181: Input Parameter:
182: . eps - the eigensolver context
184: Output Parameter:
185: . name - name of EPS method
187: Level: intermediate
189: .seealso: EPSSetType()
190: @*/
191: PetscErrorCode EPSGetType(EPS eps,EPSType *type)
192: {
196: *type = ((PetscObject)eps)->type_name;
197: return(0);
198: }
200: /*@C
201: EPSRegister - Adds a method to the eigenproblem solver package.
203: Not Collective
205: Input Parameters:
206: + name - name of a new user-defined solver
207: - function - routine to create the solver context
209: Notes:
210: EPSRegister() may be called multiple times to add several user-defined solvers.
212: Sample usage:
213: .vb
214: EPSRegister("my_solver",MySolverCreate);
215: .ve
217: Then, your solver can be chosen with the procedural interface via
218: $ EPSSetType(eps,"my_solver")
219: or at runtime via the option
220: $ -eps_type my_solver
222: Level: advanced
224: .seealso: EPSRegisterAll()
225: @*/
226: PetscErrorCode EPSRegister(const char *name,PetscErrorCode (*function)(EPS))
227: {
231: EPSInitializePackage();
232: PetscFunctionListAdd(&EPSList,name,function);
233: return(0);
234: }
236: /*@
237: EPSReset - Resets the EPS context to the initial state (prior to setup)
238: and destroys any allocated Vecs and Mats.
240: Collective on EPS
242: Input Parameter:
243: . eps - eigensolver context obtained from EPSCreate()
245: Note:
246: This can be used when a problem of different matrix size wants to be solved.
247: All options that have previously been set are preserved, so in a next use
248: the solver configuration is the same, but new sizes for matrices and vectors
249: are allowed.
251: Level: advanced
253: .seealso: EPSDestroy()
254: @*/
255: PetscErrorCode EPSReset(EPS eps)
256: {
261: if (!eps) return(0);
262: if (eps->ops->reset) { (eps->ops->reset)(eps); }
263: if (eps->st) { STReset(eps->st); }
264: VecDestroy(&eps->D);
265: BVDestroy(&eps->V);
266: BVDestroy(&eps->W);
267: VecDestroyVecs(eps->nwork,&eps->work);
268: eps->nwork = 0;
269: eps->state = EPS_STATE_INITIAL;
270: return(0);
271: }
273: /*@
274: EPSDestroy - Destroys the EPS context.
276: Collective on EPS
278: Input Parameter:
279: . eps - eigensolver context obtained from EPSCreate()
281: Level: beginner
283: .seealso: EPSCreate(), EPSSetUp(), EPSSolve()
284: @*/
285: PetscErrorCode EPSDestroy(EPS *eps)
286: {
290: if (!*eps) return(0);
292: if (--((PetscObject)(*eps))->refct > 0) { *eps = 0; return(0); }
293: EPSReset(*eps);
294: if ((*eps)->ops->destroy) { (*(*eps)->ops->destroy)(*eps); }
295: if ((*eps)->eigr) {
296: PetscFree4((*eps)->eigr,(*eps)->eigi,(*eps)->errest,(*eps)->perm);
297: }
298: if ((*eps)->rr) {
299: PetscFree2((*eps)->rr,(*eps)->ri);
300: }
301: STDestroy(&(*eps)->st);
302: RGDestroy(&(*eps)->rg);
303: DSDestroy(&(*eps)->ds);
304: DSDestroy(&(*eps)->dsts);
305: PetscFree((*eps)->sc);
306: /* just in case the initial vectors have not been used */
307: SlepcBasisDestroy_Private(&(*eps)->nds,&(*eps)->defl);
308: SlepcBasisDestroy_Private(&(*eps)->nini,&(*eps)->IS);
309: if ((*eps)->convergeddestroy) {
310: (*(*eps)->convergeddestroy)((*eps)->convergedctx);
311: }
312: EPSMonitorCancel(*eps);
313: PetscHeaderDestroy(eps);
314: return(0);
315: }
317: /*@
318: EPSSetTarget - Sets the value of the target.
320: Logically Collective on EPS
322: Input Parameters:
323: + eps - eigensolver context
324: - target - the value of the target
326: Options Database Key:
327: . -eps_target <scalar> - the value of the target
329: Notes:
330: The target is a scalar value used to determine the portion of the spectrum
331: of interest. It is used in combination with EPSSetWhichEigenpairs().
333: In the case of complex scalars, a complex value can be provided in the
334: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
335: -eps_target 1.0+2.0i
337: Level: intermediate
339: .seealso: EPSGetTarget(), EPSSetWhichEigenpairs()
340: @*/
341: PetscErrorCode EPSSetTarget(EPS eps,PetscScalar target)
342: {
348: eps->target = target;
349: if (!eps->st) { EPSGetST(eps,&eps->st); }
350: STSetDefaultShift(eps->st,target);
351: return(0);
352: }
354: /*@
355: EPSGetTarget - Gets the value of the target.
357: Not Collective
359: Input Parameter:
360: . eps - eigensolver context
362: Output Parameter:
363: . target - the value of the target
365: Note:
366: If the target was not set by the user, then zero is returned.
368: Level: intermediate
370: .seealso: EPSSetTarget()
371: @*/
372: PetscErrorCode EPSGetTarget(EPS eps,PetscScalar* target)
373: {
377: *target = eps->target;
378: return(0);
379: }
381: /*@
382: EPSSetInterval - Defines the computational interval for spectrum slicing.
384: Logically Collective on EPS
386: Input Parameters:
387: + eps - eigensolver context
388: . inta - left end of the interval
389: - intb - right end of the interval
391: Options Database Key:
392: . -eps_interval <a,b> - set [a,b] as the interval of interest
394: Notes:
395: Spectrum slicing is a technique employed for computing all eigenvalues of
396: symmetric eigenproblems in a given interval. This function provides the
397: interval to be considered. It must be used in combination with EPS_ALL, see
398: EPSSetWhichEigenpairs().
400: In the command-line option, two values must be provided. For an open interval,
401: one can give an infinite, e.g., -eps_interval 1.0,inf or -eps_interval -inf,1.0.
402: An open interval in the programmatic interface can be specified with
403: PETSC_MAX_REAL and -PETSC_MAX_REAL.
405: Level: intermediate
407: .seealso: EPSGetInterval(), EPSSetWhichEigenpairs()
408: @*/
409: PetscErrorCode EPSSetInterval(EPS eps,PetscReal inta,PetscReal intb)
410: {
415: if (inta>intb) SETERRQ(PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
416: if (eps->inta != inta || eps->intb != intb) {
417: eps->inta = inta;
418: eps->intb = intb;
419: eps->state = EPS_STATE_INITIAL;
420: }
421: return(0);
422: }
424: /*@
425: EPSGetInterval - Gets the computational interval for spectrum slicing.
427: Not Collective
429: Input Parameter:
430: . eps - eigensolver context
432: Output Parameters:
433: + inta - left end of the interval
434: - intb - right end of the interval
436: Level: intermediate
438: Note:
439: If the interval was not set by the user, then zeros are returned.
441: .seealso: EPSSetInterval()
442: @*/
443: PetscErrorCode EPSGetInterval(EPS eps,PetscReal* inta,PetscReal* intb)
444: {
447: if (inta) *inta = eps->inta;
448: if (intb) *intb = eps->intb;
449: return(0);
450: }
452: /*@
453: EPSSetST - Associates a spectral transformation object to the eigensolver.
455: Collective on EPS
457: Input Parameters:
458: + eps - eigensolver context obtained from EPSCreate()
459: - st - the spectral transformation object
461: Note:
462: Use EPSGetST() to retrieve the spectral transformation context (for example,
463: to free it at the end of the computations).
465: Level: advanced
467: .seealso: EPSGetST()
468: @*/
469: PetscErrorCode EPSSetST(EPS eps,ST st)
470: {
477: PetscObjectReference((PetscObject)st);
478: STDestroy(&eps->st);
479: eps->st = st;
480: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->st);
481: return(0);
482: }
484: /*@
485: EPSGetST - Obtain the spectral transformation (ST) object associated
486: to the eigensolver object.
488: Not Collective
490: Input Parameters:
491: . eps - eigensolver context obtained from EPSCreate()
493: Output Parameter:
494: . st - spectral transformation context
496: Level: intermediate
498: .seealso: EPSSetST()
499: @*/
500: PetscErrorCode EPSGetST(EPS eps,ST *st)
501: {
507: if (!eps->st) {
508: STCreate(PetscObjectComm((PetscObject)eps),&eps->st);
509: PetscObjectIncrementTabLevel((PetscObject)eps->st,(PetscObject)eps,0);
510: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->st);
511: PetscObjectSetOptions((PetscObject)eps->st,((PetscObject)eps)->options);
512: }
513: *st = eps->st;
514: return(0);
515: }
517: /*@
518: EPSSetBV - Associates a basis vectors object to the eigensolver.
520: Collective on EPS
522: Input Parameters:
523: + eps - eigensolver context obtained from EPSCreate()
524: - V - the basis vectors object
526: Level: advanced
528: .seealso: EPSGetBV()
529: @*/
530: PetscErrorCode EPSSetBV(EPS eps,BV V)
531: {
538: PetscObjectReference((PetscObject)V);
539: BVDestroy(&eps->V);
540: eps->V = V;
541: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->V);
542: return(0);
543: }
545: /*@
546: EPSGetBV - Obtain the basis vectors object associated to the eigensolver object.
548: Not Collective
550: Input Parameters:
551: . eps - eigensolver context obtained from EPSCreate()
553: Output Parameter:
554: . V - basis vectors context
556: Level: advanced
558: .seealso: EPSSetBV()
559: @*/
560: PetscErrorCode EPSGetBV(EPS eps,BV *V)
561: {
567: if (!eps->V) {
568: BVCreate(PetscObjectComm((PetscObject)eps),&eps->V);
569: PetscObjectIncrementTabLevel((PetscObject)eps->V,(PetscObject)eps,0);
570: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->V);
571: PetscObjectSetOptions((PetscObject)eps->V,((PetscObject)eps)->options);
572: }
573: *V = eps->V;
574: return(0);
575: }
577: /*@
578: EPSSetRG - Associates a region object to the eigensolver.
580: Collective on EPS
582: Input Parameters:
583: + eps - eigensolver context obtained from EPSCreate()
584: - rg - the region object
586: Note:
587: Use EPSGetRG() to retrieve the region context (for example,
588: to free it at the end of the computations).
590: Level: advanced
592: .seealso: EPSGetRG()
593: @*/
594: PetscErrorCode EPSSetRG(EPS eps,RG rg)
595: {
602: PetscObjectReference((PetscObject)rg);
603: RGDestroy(&eps->rg);
604: eps->rg = rg;
605: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->rg);
606: return(0);
607: }
609: /*@
610: EPSGetRG - Obtain the region object associated to the eigensolver.
612: Not Collective
614: Input Parameters:
615: . eps - eigensolver context obtained from EPSCreate()
617: Output Parameter:
618: . rg - region context
620: Level: advanced
622: .seealso: EPSSetRG()
623: @*/
624: PetscErrorCode EPSGetRG(EPS eps,RG *rg)
625: {
631: if (!eps->rg) {
632: RGCreate(PetscObjectComm((PetscObject)eps),&eps->rg);
633: PetscObjectIncrementTabLevel((PetscObject)eps->rg,(PetscObject)eps,0);
634: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->rg);
635: PetscObjectSetOptions((PetscObject)eps->rg,((PetscObject)eps)->options);
636: }
637: *rg = eps->rg;
638: return(0);
639: }
641: /*@
642: EPSSetDS - Associates a direct solver object to the eigensolver.
644: Collective on EPS
646: Input Parameters:
647: + eps - eigensolver context obtained from EPSCreate()
648: - ds - the direct solver object
650: Note:
651: Use EPSGetDS() to retrieve the direct solver context (for example,
652: to free it at the end of the computations).
654: Level: advanced
656: .seealso: EPSGetDS()
657: @*/
658: PetscErrorCode EPSSetDS(EPS eps,DS ds)
659: {
666: PetscObjectReference((PetscObject)ds);
667: DSDestroy(&eps->ds);
668: eps->ds = ds;
669: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->ds);
670: return(0);
671: }
673: /*@
674: EPSGetDS - Obtain the direct solver object associated to the eigensolver object.
676: Not Collective
678: Input Parameters:
679: . eps - eigensolver context obtained from EPSCreate()
681: Output Parameter:
682: . ds - direct solver context
684: Level: advanced
686: .seealso: EPSSetDS()
687: @*/
688: PetscErrorCode EPSGetDS(EPS eps,DS *ds)
689: {
695: if (!eps->ds) {
696: DSCreate(PetscObjectComm((PetscObject)eps),&eps->ds);
697: PetscObjectIncrementTabLevel((PetscObject)eps->ds,(PetscObject)eps,0);
698: PetscLogObjectParent((PetscObject)eps,(PetscObject)eps->ds);
699: PetscObjectSetOptions((PetscObject)eps->ds,((PetscObject)eps)->options);
700: }
701: *ds = eps->ds;
702: return(0);
703: }
705: /*@
706: EPSIsGeneralized - Ask if the EPS object corresponds to a generalized
707: eigenvalue problem.
709: Not collective
711: Input Parameter:
712: . eps - the eigenproblem solver context
714: Output Parameter:
715: . is - the answer
717: Level: intermediate
719: .seealso: EPSIsHermitian(), EPSIsPositive()
720: @*/
721: PetscErrorCode EPSIsGeneralized(EPS eps,PetscBool* is)
722: {
726: *is = eps->isgeneralized;
727: return(0);
728: }
730: /*@
731: EPSIsHermitian - Ask if the EPS object corresponds to a Hermitian
732: eigenvalue problem.
734: Not collective
736: Input Parameter:
737: . eps - the eigenproblem solver context
739: Output Parameter:
740: . is - the answer
742: Level: intermediate
744: .seealso: EPSIsGeneralized(), EPSIsPositive()
745: @*/
746: PetscErrorCode EPSIsHermitian(EPS eps,PetscBool* is)
747: {
751: *is = eps->ishermitian;
752: return(0);
753: }
755: /*@
756: EPSIsPositive - Ask if the EPS object corresponds to an eigenvalue
757: problem type that requires a positive (semi-) definite matrix B.
759: Not collective
761: Input Parameter:
762: . eps - the eigenproblem solver context
764: Output Parameter:
765: . is - the answer
767: Level: intermediate
769: .seealso: EPSIsGeneralized(), EPSIsHermitian()
770: @*/
771: PetscErrorCode EPSIsPositive(EPS eps,PetscBool* is)
772: {
776: *is = eps->ispositive;
777: return(0);
778: }