1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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 RG routines
12: */
14: #include <slepc/private/rgimpl.h> 16: PetscFunctionList RGList = 0;
17: PetscBool RGRegisterAllCalled = PETSC_FALSE;
18: PetscClassId RG_CLASSID = 0;
19: static PetscBool RGPackageInitialized = PETSC_FALSE;
21: /*@C
22: RGFinalizePackage - This function destroys everything in the Slepc interface
23: to the RG package. It is called from SlepcFinalize().
25: Level: developer
27: .seealso: SlepcFinalize()
28: @*/
29: PetscErrorCode RGFinalizePackage(void) 30: {
31: PetscFunctionListDestroy(&RGList);
32: RGPackageInitialized = PETSC_FALSE;
33: RGRegisterAllCalled = PETSC_FALSE;
34: PetscFunctionReturn(0);
35: }
37: /*@C
38: RGInitializePackage - This function initializes everything in the RG package.
39: It is called from PetscDLLibraryRegister() when using dynamic libraries, and
40: on the first call to RGCreate() when using static libraries.
42: Level: developer
44: .seealso: SlepcInitialize()
45: @*/
46: PetscErrorCode RGInitializePackage(void) 47: {
48: char logList[256];
49: PetscBool opt,pkg;
50: PetscClassId classids[1];
52: if (RGPackageInitialized) PetscFunctionReturn(0);
53: RGPackageInitialized = PETSC_TRUE;
54: /* Register Classes */
55: PetscClassIdRegister("Region",&RG_CLASSID);
56: /* Register Constructors */
57: RGRegisterAll();
58: /* Process Info */
59: classids[0] = RG_CLASSID;
60: PetscInfoProcessClass("rg",1,&classids[0]);
61: /* Process summary exclusions */
62: PetscOptionsGetString(NULL,NULL,"-log_exclude",logList,sizeof(logList),&opt);
63: if (opt) {
64: PetscStrInList("rg",logList,',',&pkg);
65: if (pkg) PetscLogEventDeactivateClass(RG_CLASSID);
66: }
67: /* Register package finalizer */
68: PetscRegisterFinalize(RGFinalizePackage);
69: PetscFunctionReturn(0);
70: }
72: /*@
73: RGCreate - Creates an RG context.
75: Collective
77: Input Parameter:
78: . comm - MPI communicator
80: Output Parameter:
81: . newrg - location to put the RG context
83: Level: beginner
85: .seealso: RGDestroy(), RG 86: @*/
87: PetscErrorCode RGCreate(MPI_Comm comm,RG *newrg) 88: {
89: RG rg;
92: *newrg = 0;
93: RGInitializePackage();
94: SlepcHeaderCreate(rg,RG_CLASSID,"RG","Region","RG",comm,RGDestroy,RGView);
95: rg->complement = PETSC_FALSE;
96: rg->sfactor = 1.0;
97: rg->osfactor = 0.0;
98: rg->data = NULL;
100: *newrg = rg;
101: PetscFunctionReturn(0);
102: }
104: /*@C
105: RGSetOptionsPrefix - Sets the prefix used for searching for all
106: RG options in the database.
108: Logically Collective on rg
110: Input Parameters:
111: + rg - the region context
112: - prefix - the prefix string to prepend to all RG option requests
114: Notes:
115: A hyphen (-) must NOT be given at the beginning of the prefix name.
116: The first character of all runtime options is AUTOMATICALLY the
117: hyphen.
119: Level: advanced
121: .seealso: RGAppendOptionsPrefix()
122: @*/
123: PetscErrorCode RGSetOptionsPrefix(RG rg,const char *prefix)124: {
126: PetscObjectSetOptionsPrefix((PetscObject)rg,prefix);
127: PetscFunctionReturn(0);
128: }
130: /*@C
131: RGAppendOptionsPrefix - Appends to the prefix used for searching for all
132: RG options in the database.
134: Logically Collective on rg
136: Input Parameters:
137: + rg - the region context
138: - prefix - the prefix string to prepend to all RG option requests
140: Notes:
141: A hyphen (-) must NOT be given at the beginning of the prefix name.
142: The first character of all runtime options is AUTOMATICALLY the hyphen.
144: Level: advanced
146: .seealso: RGSetOptionsPrefix()
147: @*/
148: PetscErrorCode RGAppendOptionsPrefix(RG rg,const char *prefix)149: {
151: PetscObjectAppendOptionsPrefix((PetscObject)rg,prefix);
152: PetscFunctionReturn(0);
153: }
155: /*@C
156: RGGetOptionsPrefix - Gets the prefix used for searching for all
157: RG options in the database.
159: Not Collective
161: Input Parameters:
162: . rg - the region context
164: Output Parameters:
165: . prefix - pointer to the prefix string used is returned
167: Note:
168: On the Fortran side, the user should pass in a string 'prefix' of
169: sufficient length to hold the prefix.
171: Level: advanced
173: .seealso: RGSetOptionsPrefix(), RGAppendOptionsPrefix()
174: @*/
175: PetscErrorCode RGGetOptionsPrefix(RG rg,const char *prefix[])176: {
179: PetscObjectGetOptionsPrefix((PetscObject)rg,prefix);
180: PetscFunctionReturn(0);
181: }
183: /*@C
184: RGSetType - Selects the type for the RG object.
186: Logically Collective on rg
188: Input Parameters:
189: + rg - the region context
190: - type - a known type
192: Level: intermediate
194: .seealso: RGGetType()
195: @*/
196: PetscErrorCode RGSetType(RG rg,RGType type)197: {
198: PetscErrorCode (*r)(RG);
199: PetscBool match;
204: PetscObjectTypeCompare((PetscObject)rg,type,&match);
205: if (match) PetscFunctionReturn(0);
207: PetscFunctionListFind(RGList,type,&r);
210: if (rg->ops->destroy) (*rg->ops->destroy)(rg);
211: PetscMemzero(rg->ops,sizeof(struct _RGOps));
213: PetscObjectChangeTypeName((PetscObject)rg,type);
214: (*r)(rg);
215: PetscFunctionReturn(0);
216: }
218: /*@C
219: RGGetType - Gets the RG type name (as a string) from the RG context.
221: Not Collective
223: Input Parameter:
224: . rg - the region context
226: Output Parameter:
227: . type - name of the region
229: Level: intermediate
231: .seealso: RGSetType()
232: @*/
233: PetscErrorCode RGGetType(RG rg,RGType *type)234: {
237: *type = ((PetscObject)rg)->type_name;
238: PetscFunctionReturn(0);
239: }
241: /*@
242: RGSetFromOptions - Sets RG options from the options database.
244: Collective on rg
246: Input Parameters:
247: . rg - the region context
249: Notes:
250: To see all options, run your program with the -help option.
252: Level: beginner
254: .seealso: RGSetOptionsPrefix()
255: @*/
256: PetscErrorCode RGSetFromOptions(RG rg)257: {
259: char type[256];
260: PetscBool flg;
261: PetscReal sfactor;
264: RGRegisterAll();
265: ierr = PetscObjectOptionsBegin((PetscObject)rg);
266: PetscOptionsFList("-rg_type","Region type","RGSetType",RGList,(char*)(((PetscObject)rg)->type_name?((PetscObject)rg)->type_name:RGINTERVAL),type,sizeof(type),&flg);
267: if (flg) RGSetType(rg,type);
268: else if (!((PetscObject)rg)->type_name) RGSetType(rg,RGINTERVAL);
270: PetscOptionsBool("-rg_complement","Whether region is complemented or not","RGSetComplement",rg->complement,&rg->complement,NULL);
272: PetscOptionsReal("-rg_scale","Scaling factor","RGSetScale",1.0,&sfactor,&flg);
273: if (flg) RGSetScale(rg,sfactor);
275: if (rg->ops->setfromoptions) (*rg->ops->setfromoptions)(PetscOptionsObject,rg);
276: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)rg);
277: ierr = PetscOptionsEnd();
278: PetscFunctionReturn(0);
279: }
281: /*@C
282: RGView - Prints the RG data structure.
284: Collective on rg
286: Input Parameters:
287: + rg - the region context
288: - viewer - optional visualization context
290: Note:
291: The available visualization contexts include
292: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
293: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
294: output where only the first processor opens
295: the file. All other processors send their
296: data to the first processor to print.
298: The user can open an alternative visualization context with
299: PetscViewerASCIIOpen() - output to a specified file.
301: Level: beginner
303: .seealso: RGCreate()
304: @*/
305: PetscErrorCode RGView(RG rg,PetscViewer viewer)306: {
307: PetscBool isdraw,isascii;
310: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rg),&viewer);
313: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
314: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
315: if (isascii) {
316: PetscObjectPrintClassNamePrefixType((PetscObject)rg,viewer);
317: if (rg->ops->view) {
318: PetscViewerASCIIPushTab(viewer);
319: (*rg->ops->view)(rg,viewer);
320: PetscViewerASCIIPopTab(viewer);
321: }
322: if (rg->complement) PetscViewerASCIIPrintf(viewer," selected region is the complement of the specified one\n");
323: if (rg->sfactor!=1.0) PetscViewerASCIIPrintf(viewer," scaling factor = %g\n",(double)rg->sfactor);
324: } else if (isdraw) {
325: if (rg->ops->view) (*rg->ops->view)(rg,viewer);
326: }
327: PetscFunctionReturn(0);
328: }
330: /*@C
331: RGViewFromOptions - View from options
333: Collective on RG335: Input Parameters:
336: + rg - the region context
337: . obj - optional object
338: - name - command line option
340: Level: intermediate
342: .seealso: RGView(), RGCreate()
343: @*/
344: PetscErrorCode RGViewFromOptions(RG rg,PetscObject obj,const char name[])345: {
347: PetscObjectViewFromOptions((PetscObject)rg,obj,name);
348: PetscFunctionReturn(0);
349: }
351: /*@
352: RGIsTrivial - Whether it is the trivial region (whole complex plane).
354: Not Collective
356: Input Parameter:
357: . rg - the region context
359: Output Parameter:
360: . trivial - true if the region is equal to the whole complex plane, e.g.,
361: an interval region with all four endpoints unbounded or an
362: ellipse with infinite radius.
364: Level: beginner
366: .seealso: RGCheckInside()
367: @*/
368: PetscErrorCode RGIsTrivial(RG rg,PetscBool *trivial)369: {
373: if (rg->ops->istrivial) (*rg->ops->istrivial)(rg,trivial);
374: else *trivial = PETSC_FALSE;
375: PetscFunctionReturn(0);
376: }
378: /*@
379: RGCheckInside - Determines if a set of given points are inside the region or not.
381: Not Collective
383: Input Parameters:
384: + rg - the region context
385: . n - number of points to check
386: . ar - array of real parts
387: - ai - array of imaginary parts
389: Output Parameter:
390: . inside - array of results (1=inside, 0=on the contour, -1=outside)
392: Note:
393: The point a is expressed as a couple of PetscScalar variables ar,ai.
394: If built with complex scalars, the point is supposed to be stored in ar,
395: otherwise ar,ai contain the real and imaginary parts, respectively.
397: If a scaling factor was set, the points are scaled before checking.
399: Level: intermediate
401: .seealso: RGSetScale(), RGSetComplement()
402: @*/
403: PetscErrorCode RGCheckInside(RG rg,PetscInt n,PetscScalar *ar,PetscScalar *ai,PetscInt *inside)404: {
405: PetscReal px,py;
406: PetscInt i;
411: #if !defined(PETSC_USE_COMPLEX)
413: #endif
416: for (i=0;i<n;i++) {
417: #if defined(PETSC_USE_COMPLEX)
418: px = PetscRealPart(ar[i]);
419: py = PetscImaginaryPart(ar[i]);
420: #else
421: px = ar[i];
422: py = ai[i];
423: #endif
424: if (PetscUnlikely(rg->sfactor != 1.0)) {
425: px /= rg->sfactor;
426: py /= rg->sfactor;
427: }
428: (*rg->ops->checkinside)(rg,px,py,inside+i);
429: if (PetscUnlikely(rg->complement)) inside[i] = -inside[i];
430: }
431: PetscFunctionReturn(0);
432: }
434: /*@
435: RGIsAxisymmetric - Determines if the region is symmetric with respect
436: to the real or imaginary axis.
438: Not Collective
440: Input Parameters:
441: + rg - the region context
442: - vertical - true if symmetry must be checked against the vertical axis
444: Output Parameter:
445: . symm - true if the region is axisymmetric
447: Note:
448: If the vertical argument is true, symmetry is checked with respect to
449: the vertical axis, otherwise with respect to the horizontal axis.
451: Level: intermediate
453: .seealso: RGCanUseConjugates()
454: @*/
455: PetscErrorCode RGIsAxisymmetric(RG rg,PetscBool vertical,PetscBool *symm)456: {
461: if (rg->ops->isaxisymmetric) (*rg->ops->isaxisymmetric)(rg,vertical,symm);
462: else *symm = PETSC_FALSE;
463: PetscFunctionReturn(0);
464: }
466: /*@
467: RGCanUseConjugates - Used in contour integral methods to determine whether
468: half of integration points can be avoided (use their conjugates).
470: Not Collective
472: Input Parameters:
473: + rg - the region context
474: - realmats - true if the problem matrices are real
476: Output Parameter:
477: . useconj - whether it is possible to use conjugates
479: Notes:
480: If some integration points are the conjugates of other points, then the
481: associated computational cost can be saved. This depends on the problem
482: matrices being real and also the region being symmetric with respect to
483: the horizontal axis. The result is false if using real arithmetic or
484: in the case of a flat region (height equal to zero).
486: Level: developer
488: .seealso: RGIsAxisymmetric()
489: @*/
490: PetscErrorCode RGCanUseConjugates(RG rg,PetscBool realmats,PetscBool *useconj)491: {
492: #if defined(PETSC_USE_COMPLEX)
493: PetscReal c,d;
494: PetscBool isaxisymm;
495: #endif
500: *useconj = PETSC_FALSE;
501: #if defined(PETSC_USE_COMPLEX)
502: if (realmats) {
503: RGIsAxisymmetric(rg,PETSC_FALSE,&isaxisymm);
504: if (isaxisymm) {
505: RGComputeBoundingBox(rg,NULL,NULL,&c,&d);
506: if (c!=d) *useconj = PETSC_TRUE;
507: }
508: }
509: #endif
510: PetscFunctionReturn(0);
511: }
513: /*@
514: RGComputeContour - Computes the coordinates of several points lying on the
515: contour of the region.
517: Not Collective
519: Input Parameters:
520: + rg - the region context
521: - n - number of points to compute
523: Output Parameters:
524: + cr - location to store real parts
525: - ci - location to store imaginary parts
527: Notes:
528: In real scalars, either cr or ci can be NULL (but not both). In complex
529: scalars, the coordinates are stored in cr, which cannot be NULL (ci is
530: not referenced).
532: Level: intermediate
534: .seealso: RGComputeBoundingBox()
535: @*/
536: PetscErrorCode RGComputeContour(RG rg,PetscInt n,PetscScalar cr[],PetscScalar ci[])537: {
538: PetscInt i;
542: #if defined(PETSC_USE_COMPLEX)
544: #else
546: #endif
548: (*rg->ops->computecontour)(rg,n,cr,ci);
549: for (i=0;i<n;i++) {
550: if (cr) cr[i] *= rg->sfactor;
551: if (ci) ci[i] *= rg->sfactor;
552: }
553: PetscFunctionReturn(0);
554: }
556: /*@
557: RGComputeBoundingBox - Determines the endpoints of a rectangle in the complex plane that
558: contains the region.
560: Not Collective
562: Input Parameter:
563: . rg - the region context
565: Output Parameters:
566: + a - left endpoint of the bounding box in the real axis
567: . b - right endpoint of the bounding box in the real axis
568: . c - bottom endpoint of the bounding box in the imaginary axis
569: - d - top endpoint of the bounding box in the imaginary axis
571: Notes:
572: The bounding box is defined as [a,b]x[c,d]. In regions that are not bounded (e.g. an
573: open interval) or with the complement flag set, it makes no sense to compute a bounding
574: box, so the return values are infinite.
576: Level: intermediate
578: .seealso: RGComputeContour()
579: @*/
580: PetscErrorCode RGComputeBoundingBox(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)581: {
585: if (rg->complement) { /* cannot compute bounding box */
586: if (a) *a = -PETSC_MAX_REAL;
587: if (b) *b = PETSC_MAX_REAL;
588: if (c) *c = -PETSC_MAX_REAL;
589: if (d) *d = PETSC_MAX_REAL;
590: } else {
591: (*rg->ops->computebbox)(rg,a,b,c,d);
592: if (a && *a!=-PETSC_MAX_REAL) *a *= rg->sfactor;
593: if (b && *b!= PETSC_MAX_REAL) *b *= rg->sfactor;
594: if (c && *c!=-PETSC_MAX_REAL) *c *= rg->sfactor;
595: if (d && *d!= PETSC_MAX_REAL) *d *= rg->sfactor;
596: }
597: PetscFunctionReturn(0);
598: }
600: /*@
601: RGComputeQuadrature - Computes the values of the parameters used in a
602: quadrature rule for a contour integral around the boundary of the region.
604: Not Collective
606: Input Parameters:
607: + rg - the region context
608: . quad - the type of quadrature
609: - n - number of quadrature points to compute
611: Output Parameters:
612: + z - quadrature points
613: . zn - normalized quadrature points
614: - w - quadrature weights
616: Notes:
617: In complex scalars, the values returned in z are often the same as those
618: computed by RGComputeContour(), but this is not the case in real scalars
619: where all output arguments are real.
621: The computed values change for different quadrature rules (trapezoidal
622: or Chebyshev).
624: Level: intermediate
626: .seealso: RGComputeContour()
627: @*/
628: PetscErrorCode RGComputeQuadrature(RG rg,RGQuadRule quad,PetscInt n,PetscScalar z[],PetscScalar zn[],PetscScalar w[])629: {
636: RGComputeContour(rg,n,z,NULL);
638: (*rg->ops->computequadrature)(rg,quad,n,z,zn,w);
639: PetscFunctionReturn(0);
640: }
642: /*@
643: RGSetComplement - Sets a flag to indicate that the region is the complement
644: of the specified one.
646: Logically Collective on rg
648: Input Parameters:
649: + rg - the region context
650: - flg - the boolean flag
652: Options Database Key:
653: . -rg_complement <bool> - Activate/deactivate the complementation of the region
655: Level: intermediate
657: .seealso: RGGetComplement()
658: @*/
659: PetscErrorCode RGSetComplement(RG rg,PetscBool flg)660: {
663: rg->complement = flg;
664: PetscFunctionReturn(0);
665: }
667: /*@
668: RGGetComplement - Gets a flag that that indicates whether the region
669: is complemented or not.
671: Not Collective
673: Input Parameter:
674: . rg - the region context
676: Output Parameter:
677: . flg - the flag
679: Level: intermediate
681: .seealso: RGSetComplement()
682: @*/
683: PetscErrorCode RGGetComplement(RG rg,PetscBool *flg)684: {
687: *flg = rg->complement;
688: PetscFunctionReturn(0);
689: }
691: /*@
692: RGSetScale - Sets the scaling factor to be used when checking that a
693: point is inside the region and when computing the contour.
695: Logically Collective on rg
697: Input Parameters:
698: + rg - the region context
699: - sfactor - the scaling factor
701: Options Database Key:
702: . -rg_scale <real> - Sets the scaling factor
704: Level: advanced
706: .seealso: RGGetScale(), RGCheckInside()
707: @*/
708: PetscErrorCode RGSetScale(RG rg,PetscReal sfactor)709: {
712: if (sfactor == PETSC_DEFAULT || sfactor == PETSC_DECIDE) rg->sfactor = 1.0;
713: else {
715: rg->sfactor = sfactor;
716: }
717: PetscFunctionReturn(0);
718: }
720: /*@
721: RGGetScale - Gets the scaling factor.
723: Not Collective
725: Input Parameter:
726: . rg - the region context
728: Output Parameter:
729: . sfactor - the scaling factor
731: Level: advanced
733: .seealso: RGSetScale()
734: @*/
735: PetscErrorCode RGGetScale(RG rg,PetscReal *sfactor)736: {
739: *sfactor = rg->sfactor;
740: PetscFunctionReturn(0);
741: }
743: /*@
744: RGPushScale - Sets an additional scaling factor, that will multiply the
745: user-defined scaling factor.
747: Logically Collective on rg
749: Input Parameters:
750: + rg - the region context
751: - sfactor - the scaling factor
753: Notes:
754: The current implementation does not allow pushing several scaling factors.
756: This is intended for internal use, for instance in polynomial eigensolvers
757: that use parameter scaling.
759: Level: developer
761: .seealso: RGPopScale(), RGSetScale()
762: @*/
763: PetscErrorCode RGPushScale(RG rg,PetscReal sfactor)764: {
769: rg->osfactor = rg->sfactor;
770: rg->sfactor *= sfactor;
771: PetscFunctionReturn(0);
772: }
774: /*@
775: RGPopScale - Pops the scaling factor set with RGPushScale().
777: Not Collective
779: Input Parameter:
780: . rg - the region context
782: Level: developer
784: .seealso: RGPushScale()
785: @*/
786: PetscErrorCode RGPopScale(RG rg)787: {
790: rg->sfactor = rg->osfactor;
791: rg->osfactor = 0.0;
792: PetscFunctionReturn(0);
793: }
795: /*@C
796: RGDestroy - Destroys RG context that was created with RGCreate().
798: Collective on rg
800: Input Parameter:
801: . rg - the region context
803: Level: beginner
805: .seealso: RGCreate()
806: @*/
807: PetscErrorCode RGDestroy(RG *rg)808: {
809: if (!*rg) PetscFunctionReturn(0);
811: if (--((PetscObject)(*rg))->refct > 0) { *rg = 0; PetscFunctionReturn(0); }
812: if ((*rg)->ops->destroy) (*(*rg)->ops->destroy)(*rg);
813: PetscHeaderDestroy(rg);
814: PetscFunctionReturn(0);
815: }
817: /*@C
818: RGRegister - Adds a region to the RG package.
820: Not collective
822: Input Parameters:
823: + name - name of a new user-defined RG824: - function - routine to create context
826: Notes:
827: RGRegister() may be called multiple times to add several user-defined regions.
829: Level: advanced
831: .seealso: RGRegisterAll()
832: @*/
833: PetscErrorCode RGRegister(const char *name,PetscErrorCode (*function)(RG))834: {
835: RGInitializePackage();
836: PetscFunctionListAdd(&RGList,name,function);
837: PetscFunctionReturn(0);
838: }