Actual source code: pepbasic.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 PEP routines
12: */
14: #include <slepc/private/pepimpl.h> /*I "slepcpep.h" I*/
16: PetscFunctionList PEPList = 0;
17: PetscBool PEPRegisterAllCalled = PETSC_FALSE;
18: PetscClassId PEP_CLASSID = 0;
19: PetscLogEvent PEP_SetUp = 0,PEP_Solve = 0,PEP_Refine = 0;
21: /*@
22: PEPCreate - Creates the default PEP context.
24: Collective on MPI_Comm
26: Input Parameter:
27: . comm - MPI communicator
29: Output Parameter:
30: . pep - location to put the PEP context
32: Note:
33: The default PEP type is PEPTOAR
35: Level: beginner
37: .seealso: PEPSetUp(), PEPSolve(), PEPDestroy(), PEP
38: @*/
39: PetscErrorCode PEPCreate(MPI_Comm comm,PEP *outpep)
40: {
42: PEP pep;
46: *outpep = 0;
47: PEPInitializePackage();
48: SlepcHeaderCreate(pep,PEP_CLASSID,"PEP","Polynomial Eigenvalue Problem","PEP",comm,PEPDestroy,PEPView);
50: pep->max_it = 0;
51: pep->nev = 1;
52: pep->ncv = 0;
53: pep->mpd = 0;
54: pep->nini = 0;
55: pep->target = 0.0;
56: pep->tol = PETSC_DEFAULT;
57: pep->conv = PEP_CONV_REL;
58: pep->stop = PEP_STOP_BASIC;
59: pep->which = (PEPWhich)0;
60: pep->basis = PEP_BASIS_MONOMIAL;
61: pep->problem_type = (PEPProblemType)0;
62: pep->scale = PEP_SCALE_NONE;
63: pep->sfactor = 1.0;
64: pep->dsfactor = 1.0;
65: pep->sits = 5;
66: pep->slambda = 1.0;
67: pep->refine = PEP_REFINE_NONE;
68: pep->npart = 1;
69: pep->rtol = PETSC_DEFAULT;
70: pep->rits = PETSC_DEFAULT;
71: pep->scheme = (PEPRefineScheme)0;
72: pep->extract = (PEPExtract)0;
73: pep->trackall = PETSC_FALSE;
75: pep->converged = PEPConvergedRelative;
76: pep->convergeduser = NULL;
77: pep->convergeddestroy= NULL;
78: pep->stopping = PEPStoppingBasic;
79: pep->stoppinguser = NULL;
80: pep->stoppingdestroy = NULL;
81: pep->convergedctx = NULL;
82: pep->stoppingctx = NULL;
83: pep->numbermonitors = 0;
85: pep->st = NULL;
86: pep->ds = NULL;
87: pep->V = NULL;
88: pep->rg = NULL;
89: pep->A = NULL;
90: pep->nmat = 0;
91: pep->Dl = NULL;
92: pep->Dr = NULL;
93: pep->IS = NULL;
94: pep->eigr = NULL;
95: pep->eigi = NULL;
96: pep->errest = NULL;
97: pep->perm = NULL;
98: pep->pbc = NULL;
99: pep->solvematcoeffs = NULL;
100: pep->nwork = 0;
101: pep->work = NULL;
102: pep->refineksp = NULL;
103: pep->refinesubc = NULL;
104: pep->data = NULL;
106: pep->state = PEP_STATE_INITIAL;
107: pep->nconv = 0;
108: pep->its = 0;
109: pep->n = 0;
110: pep->nloc = 0;
111: pep->nrma = NULL;
112: pep->sfactor_set = PETSC_FALSE;
113: pep->lineariz = PETSC_FALSE;
114: pep->reason = PEP_CONVERGED_ITERATING;
116: PetscNewLog(pep,&pep->sc);
117: *outpep = pep;
118: return(0);
119: }
121: /*@C
122: PEPSetType - Selects the particular solver to be used in the PEP object.
124: Logically Collective on PEP
126: Input Parameters:
127: + pep - the polynomial eigensolver context
128: - type - a known method
130: Options Database Key:
131: . -pep_type <method> - Sets the method; use -help for a list
132: of available methods
134: Notes:
135: See "slepc/include/slepcpep.h" for available methods. The default
136: is PEPTOAR.
138: Normally, it is best to use the PEPSetFromOptions() command and
139: then set the PEP type from the options database rather than by using
140: this routine. Using the options database provides the user with
141: maximum flexibility in evaluating the different available methods.
142: The PEPSetType() routine is provided for those situations where it
143: is necessary to set the iterative solver independently of the command
144: line or options database.
146: Level: intermediate
148: .seealso: PEPType
149: @*/
150: PetscErrorCode PEPSetType(PEP pep,PEPType type)
151: {
152: PetscErrorCode ierr,(*r)(PEP);
153: PetscBool match;
159: PetscObjectTypeCompare((PetscObject)pep,type,&match);
160: if (match) return(0);
162: PetscFunctionListFind(PEPList,type,&r);
163: if (!r) SETERRQ1(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown PEP type given: %s",type);
165: if (pep->ops->destroy) { (*pep->ops->destroy)(pep); }
166: PetscMemzero(pep->ops,sizeof(struct _PEPOps));
168: pep->state = PEP_STATE_INITIAL;
169: PetscObjectChangeTypeName((PetscObject)pep,type);
170: (*r)(pep);
171: return(0);
172: }
174: /*@C
175: PEPGetType - Gets the PEP type as a string from the PEP object.
177: Not Collective
179: Input Parameter:
180: . pep - the eigensolver context
182: Output Parameter:
183: . name - name of PEP method
185: Level: intermediate
187: .seealso: PEPSetType()
188: @*/
189: PetscErrorCode PEPGetType(PEP pep,PEPType *type)
190: {
194: *type = ((PetscObject)pep)->type_name;
195: return(0);
196: }
198: /*@C
199: PEPRegister - Adds a method to the polynomial eigenproblem solver package.
201: Not Collective
203: Input Parameters:
204: + name - name of a new user-defined solver
205: - function - routine to create the solver context
207: Notes:
208: PEPRegister() may be called multiple times to add several user-defined solvers.
210: Sample usage:
211: .vb
212: PEPRegister("my_solver",MySolverCreate);
213: .ve
215: Then, your solver can be chosen with the procedural interface via
216: $ PEPSetType(pep,"my_solver")
217: or at runtime via the option
218: $ -pep_type my_solver
220: Level: advanced
222: .seealso: PEPRegisterAll()
223: @*/
224: PetscErrorCode PEPRegister(const char *name,PetscErrorCode (*function)(PEP))
225: {
229: PEPInitializePackage();
230: PetscFunctionListAdd(&PEPList,name,function);
231: return(0);
232: }
234: /*@
235: PEPReset - Resets the PEP context to the initial state (prior to setup)
236: and destroys any allocated Vecs and Mats.
238: Collective on PEP
240: Input Parameter:
241: . pep - eigensolver context obtained from PEPCreate()
243: Level: advanced
245: .seealso: PEPDestroy()
246: @*/
247: PetscErrorCode PEPReset(PEP pep)
248: {
253: if (!pep) return(0);
254: if (pep->ops->reset) { (pep->ops->reset)(pep); }
255: if (pep->st) { STReset(pep->st); }
256: if (pep->refineksp) { KSPReset(pep->refineksp); }
257: if (pep->nmat) {
258: MatDestroyMatrices(pep->nmat,&pep->A);
259: PetscFree2(pep->pbc,pep->nrma);
260: PetscFree(pep->solvematcoeffs);
261: pep->nmat = 0;
262: }
263: VecDestroy(&pep->Dl);
264: VecDestroy(&pep->Dr);
265: BVDestroy(&pep->V);
266: VecDestroyVecs(pep->nwork,&pep->work);
267: pep->nwork = 0;
268: pep->state = PEP_STATE_INITIAL;
269: return(0);
270: }
272: /*@
273: PEPDestroy - Destroys the PEP context.
275: Collective on PEP
277: Input Parameter:
278: . pep - eigensolver context obtained from PEPCreate()
280: Level: beginner
282: .seealso: PEPCreate(), PEPSetUp(), PEPSolve()
283: @*/
284: PetscErrorCode PEPDestroy(PEP *pep)
285: {
289: if (!*pep) return(0);
291: if (--((PetscObject)(*pep))->refct > 0) { *pep = 0; return(0); }
292: PEPReset(*pep);
293: if ((*pep)->ops->destroy) { (*(*pep)->ops->destroy)(*pep); }
294: if ((*pep)->eigr) {
295: PetscFree4((*pep)->eigr,(*pep)->eigi,(*pep)->errest,(*pep)->perm);
296: }
297: STDestroy(&(*pep)->st);
298: RGDestroy(&(*pep)->rg);
299: DSDestroy(&(*pep)->ds);
300: KSPDestroy(&(*pep)->refineksp);
301: PetscSubcommDestroy(&(*pep)->refinesubc);
302: PetscFree((*pep)->sc);
303: /* just in case the initial vectors have not been used */
304: SlepcBasisDestroy_Private(&(*pep)->nini,&(*pep)->IS);
305: if ((*pep)->convergeddestroy) {
306: (*(*pep)->convergeddestroy)((*pep)->convergedctx);
307: }
308: PEPMonitorCancel(*pep);
309: PetscHeaderDestroy(pep);
310: return(0);
311: }
313: /*@
314: PEPSetBV - Associates a basis vectors object to the polynomial eigensolver.
316: Collective on PEP
318: Input Parameters:
319: + pep - eigensolver context obtained from PEPCreate()
320: - bv - the basis vectors object
322: Note:
323: Use PEPGetBV() to retrieve the basis vectors context (for example,
324: to free it at the end of the computations).
326: Level: advanced
328: .seealso: PEPGetBV()
329: @*/
330: PetscErrorCode PEPSetBV(PEP pep,BV bv)
331: {
338: PetscObjectReference((PetscObject)bv);
339: BVDestroy(&pep->V);
340: pep->V = bv;
341: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->V);
342: return(0);
343: }
345: /*@
346: PEPGetBV - Obtain the basis vectors object associated to the polynomial
347: eigensolver object.
349: Not Collective
351: Input Parameters:
352: . pep - eigensolver context obtained from PEPCreate()
354: Output Parameter:
355: . bv - basis vectors context
357: Level: advanced
359: .seealso: PEPSetBV()
360: @*/
361: PetscErrorCode PEPGetBV(PEP pep,BV *bv)
362: {
368: if (!pep->V) {
369: BVCreate(PetscObjectComm((PetscObject)pep),&pep->V);
370: PetscObjectIncrementTabLevel((PetscObject)pep->V,(PetscObject)pep,0);
371: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->V);
372: PetscObjectSetOptions((PetscObject)pep->V,((PetscObject)pep)->options);
373: }
374: *bv = pep->V;
375: return(0);
376: }
378: /*@
379: PEPSetRG - Associates a region object to the polynomial eigensolver.
381: Collective on PEP
383: Input Parameters:
384: + pep - eigensolver context obtained from PEPCreate()
385: - rg - the region object
387: Note:
388: Use PEPGetRG() to retrieve the region context (for example,
389: to free it at the end of the computations).
391: Level: advanced
393: .seealso: PEPGetRG()
394: @*/
395: PetscErrorCode PEPSetRG(PEP pep,RG rg)
396: {
403: PetscObjectReference((PetscObject)rg);
404: RGDestroy(&pep->rg);
405: pep->rg = rg;
406: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
407: return(0);
408: }
410: /*@
411: PEPGetRG - Obtain the region object associated to the
412: polynomial eigensolver object.
414: Not Collective
416: Input Parameters:
417: . pep - eigensolver context obtained from PEPCreate()
419: Output Parameter:
420: . rg - region context
422: Level: advanced
424: .seealso: PEPSetRG()
425: @*/
426: PetscErrorCode PEPGetRG(PEP pep,RG *rg)
427: {
433: if (!pep->rg) {
434: RGCreate(PetscObjectComm((PetscObject)pep),&pep->rg);
435: PetscObjectIncrementTabLevel((PetscObject)pep->rg,(PetscObject)pep,0);
436: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
437: PetscObjectSetOptions((PetscObject)pep->rg,((PetscObject)pep)->options);
438: }
439: *rg = pep->rg;
440: return(0);
441: }
443: /*@
444: PEPSetDS - Associates a direct solver object to the polynomial eigensolver.
446: Collective on PEP
448: Input Parameters:
449: + pep - eigensolver context obtained from PEPCreate()
450: - ds - the direct solver object
452: Note:
453: Use PEPGetDS() to retrieve the direct solver context (for example,
454: to free it at the end of the computations).
456: Level: advanced
458: .seealso: PEPGetDS()
459: @*/
460: PetscErrorCode PEPSetDS(PEP pep,DS ds)
461: {
468: PetscObjectReference((PetscObject)ds);
469: DSDestroy(&pep->ds);
470: pep->ds = ds;
471: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
472: return(0);
473: }
475: /*@
476: PEPGetDS - Obtain the direct solver object associated to the
477: polynomial eigensolver object.
479: Not Collective
481: Input Parameters:
482: . pep - eigensolver context obtained from PEPCreate()
484: Output Parameter:
485: . ds - direct solver context
487: Level: advanced
489: .seealso: PEPSetDS()
490: @*/
491: PetscErrorCode PEPGetDS(PEP pep,DS *ds)
492: {
498: if (!pep->ds) {
499: DSCreate(PetscObjectComm((PetscObject)pep),&pep->ds);
500: PetscObjectIncrementTabLevel((PetscObject)pep->ds,(PetscObject)pep,0);
501: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
502: PetscObjectSetOptions((PetscObject)pep->ds,((PetscObject)pep)->options);
503: }
504: *ds = pep->ds;
505: return(0);
506: }
508: /*@
509: PEPSetST - Associates a spectral transformation object to the eigensolver.
511: Collective on PEP
513: Input Parameters:
514: + pep - eigensolver context obtained from PEPCreate()
515: - st - the spectral transformation object
517: Note:
518: Use PEPGetST() to retrieve the spectral transformation context (for example,
519: to free it at the end of the computations).
521: Level: advanced
523: .seealso: PEPGetST()
524: @*/
525: PetscErrorCode PEPSetST(PEP pep,ST st)
526: {
533: PetscObjectReference((PetscObject)st);
534: STDestroy(&pep->st);
535: pep->st = st;
536: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
537: return(0);
538: }
540: /*@
541: PEPGetST - Obtain the spectral transformation (ST) object associated
542: to the eigensolver object.
544: Not Collective
546: Input Parameters:
547: . pep - eigensolver context obtained from PEPCreate()
549: Output Parameter:
550: . st - spectral transformation context
552: Level: intermediate
554: .seealso: PEPSetST()
555: @*/
556: PetscErrorCode PEPGetST(PEP pep,ST *st)
557: {
563: if (!pep->st) {
564: STCreate(PetscObjectComm((PetscObject)pep),&pep->st);
565: PetscObjectIncrementTabLevel((PetscObject)pep->st,(PetscObject)pep,0);
566: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
567: PetscObjectSetOptions((PetscObject)pep->st,((PetscObject)pep)->options);
568: }
569: *st = pep->st;
570: return(0);
571: }
573: /*@
574: PEPRefineGetKSP - Obtain the ksp object used by the eigensolver
575: object in the refinement phase.
577: Not Collective
579: Input Parameters:
580: . pep - eigensolver context obtained from PEPCreate()
582: Output Parameter:
583: . ksp - ksp context
585: Level: advanced
587: .seealso: PEPSetRefine()
588: @*/
589: PetscErrorCode PEPRefineGetKSP(PEP pep,KSP *ksp)
590: {
596: if (!pep->refineksp) {
597: if (pep->npart>1) {
598: /* Split in subcomunicators */
599: PetscSubcommCreate(PetscObjectComm((PetscObject)pep),&pep->refinesubc);
600: PetscSubcommSetNumber(pep->refinesubc,pep->npart);
601: PetscSubcommSetType(pep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
602: PetscLogObjectMemory((PetscObject)pep,sizeof(PetscSubcomm));
603: }
604: KSPCreate((pep->npart==1)?PetscObjectComm((PetscObject)pep):PetscSubcommChild(pep->refinesubc),&pep->refineksp);
605: PetscObjectIncrementTabLevel((PetscObject)pep->refineksp,(PetscObject)pep,0);
606: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->refineksp);
607: PetscObjectSetOptions((PetscObject)pep->refineksp,((PetscObject)pep)->options);
608: KSPSetOptionsPrefix(*ksp,((PetscObject)pep)->prefix);
609: KSPAppendOptionsPrefix(*ksp,"pep_refine_");
610: KSPSetTolerances(pep->refineksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
611: }
612: *ksp = pep->refineksp;
613: return(0);
614: }
616: /*@
617: PEPSetTarget - Sets the value of the target.
619: Logically Collective on PEP
621: Input Parameters:
622: + pep - eigensolver context
623: - target - the value of the target
625: Options Database Key:
626: . -pep_target <scalar> - the value of the target
628: Notes:
629: The target is a scalar value used to determine the portion of the spectrum
630: of interest. It is used in combination with PEPSetWhichEigenpairs().
632: In the case of complex scalars, a complex value can be provided in the
633: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
634: -pep_target 1.0+2.0i
636: Level: intermediate
638: .seealso: PEPGetTarget(), PEPSetWhichEigenpairs()
639: @*/
640: PetscErrorCode PEPSetTarget(PEP pep,PetscScalar target)
641: {
647: pep->target = target;
648: if (!pep->st) { PEPGetST(pep,&pep->st); }
649: STSetDefaultShift(pep->st,target);
650: return(0);
651: }
653: /*@
654: PEPGetTarget - Gets the value of the target.
656: Not Collective
658: Input Parameter:
659: . pep - eigensolver context
661: Output Parameter:
662: . target - the value of the target
664: Note:
665: If the target was not set by the user, then zero is returned.
667: Level: intermediate
669: .seealso: PEPSetTarget()
670: @*/
671: PetscErrorCode PEPGetTarget(PEP pep,PetscScalar* target)
672: {
676: *target = pep->target;
677: return(0);
678: }
680: /*@
681: PEPSetInterval - Defines the computational interval for spectrum slicing.
683: Logically Collective on PEP
685: Input Parameters:
686: + pep - eigensolver context
687: . inta - left end of the interval
688: - intb - right end of the interval
690: Options Database Key:
691: . -pep_interval <a,b> - set [a,b] as the interval of interest
693: Notes:
694: Spectrum slicing is a technique employed for computing all eigenvalues of
695: symmetric eigenproblems in a given interval. This function provides the
696: interval to be considered. It must be used in combination with PEP_ALL, see
697: PEPSetWhichEigenpairs().
699: In the command-line option, two values must be provided. For an open interval,
700: one can give an infinite, e.g., -pep_interval 1.0,inf or -pep_interval -inf,1.0.
701: An open interval in the programmatic interface can be specified with
702: PETSC_MAX_REAL and -PETSC_MAX_REAL.
704: Level: intermediate
706: .seealso: PEPGetInterval(), PEPSetWhichEigenpairs()
707: @*/
708: PetscErrorCode PEPSetInterval(PEP pep,PetscReal inta,PetscReal intb)
709: {
714: if (inta>intb) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
715: if (pep->inta != inta || pep->intb != intb) {
716: pep->inta = inta;
717: pep->intb = intb;
718: pep->state = PEP_STATE_INITIAL;
719: }
720: return(0);
721: }
723: /*@
724: PEPGetInterval - Gets the computational interval for spectrum slicing.
726: Not Collective
728: Input Parameter:
729: . pep - eigensolver context
731: Output Parameters:
732: + inta - left end of the interval
733: - intb - right end of the interval
735: Level: intermediate
737: Note:
738: If the interval was not set by the user, then zeros are returned.
740: .seealso: PEPSetInterval()
741: @*/
742: PetscErrorCode PEPGetInterval(PEP pep,PetscReal* inta,PetscReal* intb)
743: {
746: if (inta) *inta = pep->inta;
747: if (intb) *intb = pep->intb;
748: return(0);
749: }