Actual source code: pepbasic.c
slepc-3.10.1 2018-10-23
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2018, 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: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->V);
371: }
372: *bv = pep->V;
373: return(0);
374: }
376: /*@
377: PEPSetRG - Associates a region object to the polynomial eigensolver.
379: Collective on PEP
381: Input Parameters:
382: + pep - eigensolver context obtained from PEPCreate()
383: - rg - the region object
385: Note:
386: Use PEPGetRG() to retrieve the region context (for example,
387: to free it at the end of the computations).
389: Level: advanced
391: .seealso: PEPGetRG()
392: @*/
393: PetscErrorCode PEPSetRG(PEP pep,RG rg)
394: {
401: PetscObjectReference((PetscObject)rg);
402: RGDestroy(&pep->rg);
403: pep->rg = rg;
404: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
405: return(0);
406: }
408: /*@
409: PEPGetRG - Obtain the region object associated to the
410: polynomial eigensolver object.
412: Not Collective
414: Input Parameters:
415: . pep - eigensolver context obtained from PEPCreate()
417: Output Parameter:
418: . rg - region context
420: Level: advanced
422: .seealso: PEPSetRG()
423: @*/
424: PetscErrorCode PEPGetRG(PEP pep,RG *rg)
425: {
431: if (!pep->rg) {
432: RGCreate(PetscObjectComm((PetscObject)pep),&pep->rg);
433: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->rg);
434: }
435: *rg = pep->rg;
436: return(0);
437: }
439: /*@
440: PEPSetDS - Associates a direct solver object to the polynomial eigensolver.
442: Collective on PEP
444: Input Parameters:
445: + pep - eigensolver context obtained from PEPCreate()
446: - ds - the direct solver object
448: Note:
449: Use PEPGetDS() to retrieve the direct solver context (for example,
450: to free it at the end of the computations).
452: Level: advanced
454: .seealso: PEPGetDS()
455: @*/
456: PetscErrorCode PEPSetDS(PEP pep,DS ds)
457: {
464: PetscObjectReference((PetscObject)ds);
465: DSDestroy(&pep->ds);
466: pep->ds = ds;
467: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
468: return(0);
469: }
471: /*@
472: PEPGetDS - Obtain the direct solver object associated to the
473: polynomial eigensolver object.
475: Not Collective
477: Input Parameters:
478: . pep - eigensolver context obtained from PEPCreate()
480: Output Parameter:
481: . ds - direct solver context
483: Level: advanced
485: .seealso: PEPSetDS()
486: @*/
487: PetscErrorCode PEPGetDS(PEP pep,DS *ds)
488: {
494: if (!pep->ds) {
495: DSCreate(PetscObjectComm((PetscObject)pep),&pep->ds);
496: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->ds);
497: }
498: *ds = pep->ds;
499: return(0);
500: }
502: /*@
503: PEPSetST - Associates a spectral transformation object to the eigensolver.
505: Collective on PEP
507: Input Parameters:
508: + pep - eigensolver context obtained from PEPCreate()
509: - st - the spectral transformation object
511: Note:
512: Use PEPGetST() to retrieve the spectral transformation context (for example,
513: to free it at the end of the computations).
515: Level: advanced
517: .seealso: PEPGetST()
518: @*/
519: PetscErrorCode PEPSetST(PEP pep,ST st)
520: {
527: PetscObjectReference((PetscObject)st);
528: STDestroy(&pep->st);
529: pep->st = st;
530: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
531: return(0);
532: }
534: /*@
535: PEPGetST - Obtain the spectral transformation (ST) object associated
536: to the eigensolver object.
538: Not Collective
540: Input Parameters:
541: . pep - eigensolver context obtained from PEPCreate()
543: Output Parameter:
544: . st - spectral transformation context
546: Level: intermediate
548: .seealso: PEPSetST()
549: @*/
550: PetscErrorCode PEPGetST(PEP pep,ST *st)
551: {
557: if (!pep->st) {
558: STCreate(PetscObjectComm((PetscObject)pep),&pep->st);
559: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->st);
560: }
561: *st = pep->st;
562: return(0);
563: }
565: /*@
566: PEPRefineGetKSP - Obtain the ksp object used by the eigensolver
567: object in the refinement phase.
569: Not Collective
571: Input Parameters:
572: . pep - eigensolver context obtained from PEPCreate()
574: Output Parameter:
575: . ksp - ksp context
577: Level: advanced
579: .seealso: PEPSetRefine()
580: @*/
581: PetscErrorCode PEPRefineGetKSP(PEP pep,KSP *ksp)
582: {
588: if (!pep->refineksp) {
589: if (pep->npart>1) {
590: /* Split in subcomunicators */
591: PetscSubcommCreate(PetscObjectComm((PetscObject)pep),&pep->refinesubc);
592: PetscSubcommSetNumber(pep->refinesubc,pep->npart);
593: PetscSubcommSetType(pep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
594: PetscLogObjectMemory((PetscObject)pep,sizeof(PetscSubcomm));
595: }
596: KSPCreate((pep->npart==1)?PetscObjectComm((PetscObject)pep):PetscSubcommChild(pep->refinesubc),&pep->refineksp);
597: PetscLogObjectParent((PetscObject)pep,(PetscObject)pep->refineksp);
598: KSPSetOptionsPrefix(*ksp,((PetscObject)pep)->prefix);
599: KSPAppendOptionsPrefix(*ksp,"pep_refine_");
600: KSPSetTolerances(pep->refineksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
601: }
602: *ksp = pep->refineksp;
603: return(0);
604: }
606: /*@
607: PEPSetTarget - Sets the value of the target.
609: Logically Collective on PEP
611: Input Parameters:
612: + pep - eigensolver context
613: - target - the value of the target
615: Options Database Key:
616: . -pep_target <scalar> - the value of the target
618: Notes:
619: The target is a scalar value used to determine the portion of the spectrum
620: of interest. It is used in combination with PEPSetWhichEigenpairs().
622: In the case of complex scalars, a complex value can be provided in the
623: command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
624: -pep_target 1.0+2.0i
626: Level: intermediate
628: .seealso: PEPGetTarget(), PEPSetWhichEigenpairs()
629: @*/
630: PetscErrorCode PEPSetTarget(PEP pep,PetscScalar target)
631: {
637: pep->target = target;
638: if (!pep->st) { PEPGetST(pep,&pep->st); }
639: STSetDefaultShift(pep->st,target);
640: return(0);
641: }
643: /*@
644: PEPGetTarget - Gets the value of the target.
646: Not Collective
648: Input Parameter:
649: . pep - eigensolver context
651: Output Parameter:
652: . target - the value of the target
654: Note:
655: If the target was not set by the user, then zero is returned.
657: Level: intermediate
659: .seealso: PEPSetTarget()
660: @*/
661: PetscErrorCode PEPGetTarget(PEP pep,PetscScalar* target)
662: {
666: *target = pep->target;
667: return(0);
668: }
670: /*@
671: PEPSetInterval - Defines the computational interval for spectrum slicing.
673: Logically Collective on PEP
675: Input Parameters:
676: + pep - eigensolver context
677: . inta - left end of the interval
678: - intb - right end of the interval
680: Options Database Key:
681: . -pep_interval <a,b> - set [a,b] as the interval of interest
683: Notes:
684: Spectrum slicing is a technique employed for computing all eigenvalues of
685: symmetric eigenproblems in a given interval. This function provides the
686: interval to be considered. It must be used in combination with PEP_ALL, see
687: PEPSetWhichEigenpairs().
689: In the command-line option, two values must be provided. For an open interval,
690: one can give an infinite, e.g., -pep_interval 1.0,inf or -pep_interval -inf,1.0.
691: An open interval in the programmatic interface can be specified with
692: PETSC_MAX_REAL and -PETSC_MAX_REAL.
694: Level: intermediate
696: .seealso: PEPGetInterval(), PEPSetWhichEigenpairs()
697: @*/
698: PetscErrorCode PEPSetInterval(PEP pep,PetscReal inta,PetscReal intb)
699: {
704: if (inta>intb) SETERRQ(PetscObjectComm((PetscObject)pep),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be inta<intb");
705: if (pep->inta != inta || pep->intb != intb) {
706: pep->inta = inta;
707: pep->intb = intb;
708: pep->state = PEP_STATE_INITIAL;
709: }
710: return(0);
711: }
713: /*@
714: PEPGetInterval - Gets the computational interval for spectrum slicing.
716: Not Collective
718: Input Parameter:
719: . pep - eigensolver context
721: Output Parameters:
722: + inta - left end of the interval
723: - intb - right end of the interval
725: Level: intermediate
727: Note:
728: If the interval was not set by the user, then zeros are returned.
730: .seealso: PEPSetInterval()
731: @*/
732: PetscErrorCode PEPGetInterval(PEP pep,PetscReal* inta,PetscReal* intb)
733: {
736: if (inta) *inta = pep->inta;
737: if (intb) *intb = pep->intb;
738: return(0);
739: }