Actual source code: rgbasic.c
slepc-3.18.1 2022-11-02
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 = NULL;
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: return 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) return 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: return 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 = NULL;
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: return 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: return 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: return 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: return 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) return 0;
207: PetscFunctionListFind(RGList,type,&r);
210: PetscTryTypeMethod(rg,destroy);
211: PetscMemzero(rg->ops,sizeof(struct _RGOps));
213: PetscObjectChangeTypeName((PetscObject)rg,type);
214: (*r)(rg);
215: return 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: return 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: {
258: char type[256];
259: PetscBool flg;
260: PetscReal sfactor;
263: RGRegisterAll();
264: PetscObjectOptionsBegin((PetscObject)rg);
265: PetscOptionsFList("-rg_type","Region type","RGSetType",RGList,(char*)(((PetscObject)rg)->type_name?((PetscObject)rg)->type_name:RGINTERVAL),type,sizeof(type),&flg);
266: if (flg) RGSetType(rg,type);
267: else if (!((PetscObject)rg)->type_name) RGSetType(rg,RGINTERVAL);
269: PetscOptionsBool("-rg_complement","Whether region is complemented or not","RGSetComplement",rg->complement,&rg->complement,NULL);
271: PetscOptionsReal("-rg_scale","Scaling factor","RGSetScale",1.0,&sfactor,&flg);
272: if (flg) RGSetScale(rg,sfactor);
274: PetscTryTypeMethod(rg,setfromoptions,PetscOptionsObject);
275: PetscObjectProcessOptionsHandlers((PetscObject)rg,PetscOptionsObject);
276: PetscOptionsEnd();
277: return 0;
278: }
280: /*@C
281: RGView - Prints the RG data structure.
283: Collective on rg
285: Input Parameters:
286: + rg - the region context
287: - viewer - optional visualization context
289: Note:
290: The available visualization contexts include
291: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
292: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
293: output where only the first processor opens
294: the file. All other processors send their
295: data to the first processor to print.
297: The user can open an alternative visualization context with
298: PetscViewerASCIIOpen() - output to a specified file.
300: Level: beginner
302: .seealso: RGCreate()
303: @*/
304: PetscErrorCode RGView(RG rg,PetscViewer viewer)
305: {
306: PetscBool isdraw,isascii;
309: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)rg),&viewer);
312: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
313: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
314: if (isascii) {
315: PetscObjectPrintClassNamePrefixType((PetscObject)rg,viewer);
316: PetscViewerASCIIPushTab(viewer);
317: PetscTryTypeMethod(rg,view,viewer);
318: PetscViewerASCIIPopTab(viewer);
319: if (rg->complement) PetscViewerASCIIPrintf(viewer," selected region is the complement of the specified one\n");
320: if (rg->sfactor!=1.0) PetscViewerASCIIPrintf(viewer," scaling factor = %g\n",(double)rg->sfactor);
321: } else if (isdraw) PetscTryTypeMethod(rg,view,viewer);
322: return 0;
323: }
325: /*@C
326: RGViewFromOptions - View from options
328: Collective on RG
330: Input Parameters:
331: + rg - the region context
332: . obj - optional object
333: - name - command line option
335: Level: intermediate
337: .seealso: RGView(), RGCreate()
338: @*/
339: PetscErrorCode RGViewFromOptions(RG rg,PetscObject obj,const char name[])
340: {
342: PetscObjectViewFromOptions((PetscObject)rg,obj,name);
343: return 0;
344: }
346: /*@
347: RGIsTrivial - Whether it is the trivial region (whole complex plane).
349: Not Collective
351: Input Parameter:
352: . rg - the region context
354: Output Parameter:
355: . trivial - true if the region is equal to the whole complex plane, e.g.,
356: an interval region with all four endpoints unbounded or an
357: ellipse with infinite radius.
359: Level: beginner
361: .seealso: RGCheckInside()
362: @*/
363: PetscErrorCode RGIsTrivial(RG rg,PetscBool *trivial)
364: {
368: *trivial = PETSC_FALSE;
369: PetscTryTypeMethod(rg,istrivial,trivial);
370: return 0;
371: }
373: /*@
374: RGCheckInside - Determines if a set of given points are inside the region or not.
376: Not Collective
378: Input Parameters:
379: + rg - the region context
380: . n - number of points to check
381: . ar - array of real parts
382: - ai - array of imaginary parts
384: Output Parameter:
385: . inside - array of results (1=inside, 0=on the contour, -1=outside)
387: Note:
388: The point a is expressed as a couple of PetscScalar variables ar,ai.
389: If built with complex scalars, the point is supposed to be stored in ar,
390: otherwise ar,ai contain the real and imaginary parts, respectively.
392: If a scaling factor was set, the points are scaled before checking.
394: Level: intermediate
396: .seealso: RGSetScale(), RGSetComplement()
397: @*/
398: PetscErrorCode RGCheckInside(RG rg,PetscInt n,PetscScalar *ar,PetscScalar *ai,PetscInt *inside)
399: {
400: PetscReal px,py;
401: PetscInt i;
406: #if !defined(PETSC_USE_COMPLEX)
408: #endif
411: for (i=0;i<n;i++) {
412: #if defined(PETSC_USE_COMPLEX)
413: px = PetscRealPart(ar[i]);
414: py = PetscImaginaryPart(ar[i]);
415: #else
416: px = ar[i];
417: py = ai[i];
418: #endif
419: if (PetscUnlikely(rg->sfactor != 1.0)) {
420: px /= rg->sfactor;
421: py /= rg->sfactor;
422: }
423: PetscUseTypeMethod(rg,checkinside,px,py,inside+i);
424: if (PetscUnlikely(rg->complement)) inside[i] = -inside[i];
425: }
426: return 0;
427: }
429: /*@
430: RGIsAxisymmetric - Determines if the region is symmetric with respect
431: to the real or imaginary axis.
433: Not Collective
435: Input Parameters:
436: + rg - the region context
437: - vertical - true if symmetry must be checked against the vertical axis
439: Output Parameter:
440: . symm - true if the region is axisymmetric
442: Note:
443: If the vertical argument is true, symmetry is checked with respect to
444: the vertical axis, otherwise with respect to the horizontal axis.
446: Level: intermediate
448: .seealso: RGCanUseConjugates()
449: @*/
450: PetscErrorCode RGIsAxisymmetric(RG rg,PetscBool vertical,PetscBool *symm)
451: {
455: *symm = PETSC_FALSE;
456: PetscTryTypeMethod(rg,isaxisymmetric,vertical,symm);
457: return 0;
458: }
460: /*@
461: RGCanUseConjugates - Used in contour integral methods to determine whether
462: half of integration points can be avoided (use their conjugates).
464: Not Collective
466: Input Parameters:
467: + rg - the region context
468: - realmats - true if the problem matrices are real
470: Output Parameter:
471: . useconj - whether it is possible to use conjugates
473: Notes:
474: If some integration points are the conjugates of other points, then the
475: associated computational cost can be saved. This depends on the problem
476: matrices being real and also the region being symmetric with respect to
477: the horizontal axis. The result is false if using real arithmetic or
478: in the case of a flat region (height equal to zero).
480: Level: developer
482: .seealso: RGIsAxisymmetric()
483: @*/
484: PetscErrorCode RGCanUseConjugates(RG rg,PetscBool realmats,PetscBool *useconj)
485: {
486: #if defined(PETSC_USE_COMPLEX)
487: PetscReal c,d;
488: PetscBool isaxisymm;
489: #endif
494: *useconj = PETSC_FALSE;
495: #if defined(PETSC_USE_COMPLEX)
496: if (realmats) {
497: RGIsAxisymmetric(rg,PETSC_FALSE,&isaxisymm);
498: if (isaxisymm) {
499: RGComputeBoundingBox(rg,NULL,NULL,&c,&d);
500: if (c!=d) *useconj = PETSC_TRUE;
501: }
502: }
503: #endif
504: return 0;
505: }
507: /*@
508: RGComputeContour - Computes the coordinates of several points lying on the
509: contour of the region.
511: Not Collective
513: Input Parameters:
514: + rg - the region context
515: - n - number of points to compute
517: Output Parameters:
518: + cr - location to store real parts
519: - ci - location to store imaginary parts
521: Notes:
522: In real scalars, either cr or ci can be NULL (but not both). In complex
523: scalars, the coordinates are stored in cr, which cannot be NULL (ci is
524: not referenced).
526: Level: intermediate
528: .seealso: RGComputeBoundingBox()
529: @*/
530: PetscErrorCode RGComputeContour(RG rg,PetscInt n,PetscScalar cr[],PetscScalar ci[])
531: {
532: PetscInt i;
536: #if defined(PETSC_USE_COMPLEX)
538: #else
540: #endif
542: PetscUseTypeMethod(rg,computecontour,n,cr,ci);
543: for (i=0;i<n;i++) {
544: if (cr) cr[i] *= rg->sfactor;
545: if (ci) ci[i] *= rg->sfactor;
546: }
547: return 0;
548: }
550: /*@
551: RGComputeBoundingBox - Determines the endpoints of a rectangle in the complex plane that
552: contains the region.
554: Not Collective
556: Input Parameter:
557: . rg - the region context
559: Output Parameters:
560: + a - left endpoint of the bounding box in the real axis
561: . b - right endpoint of the bounding box in the real axis
562: . c - bottom endpoint of the bounding box in the imaginary axis
563: - d - top endpoint of the bounding box in the imaginary axis
565: Notes:
566: The bounding box is defined as [a,b]x[c,d]. In regions that are not bounded (e.g. an
567: open interval) or with the complement flag set, it makes no sense to compute a bounding
568: box, so the return values are infinite.
570: Level: intermediate
572: .seealso: RGComputeContour()
573: @*/
574: PetscErrorCode RGComputeBoundingBox(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
575: {
579: if (rg->complement) { /* cannot compute bounding box */
580: if (a) *a = -PETSC_MAX_REAL;
581: if (b) *b = PETSC_MAX_REAL;
582: if (c) *c = -PETSC_MAX_REAL;
583: if (d) *d = PETSC_MAX_REAL;
584: } else {
585: PetscUseTypeMethod(rg,computebbox,a,b,c,d);
586: if (a && *a!=-PETSC_MAX_REAL) *a *= rg->sfactor;
587: if (b && *b!= PETSC_MAX_REAL) *b *= rg->sfactor;
588: if (c && *c!=-PETSC_MAX_REAL) *c *= rg->sfactor;
589: if (d && *d!= PETSC_MAX_REAL) *d *= rg->sfactor;
590: }
591: return 0;
592: }
594: /*@
595: RGComputeQuadrature - Computes the values of the parameters used in a
596: quadrature rule for a contour integral around the boundary of the region.
598: Not Collective
600: Input Parameters:
601: + rg - the region context
602: . quad - the type of quadrature
603: - n - number of quadrature points to compute
605: Output Parameters:
606: + z - quadrature points
607: . zn - normalized quadrature points
608: - w - quadrature weights
610: Notes:
611: In complex scalars, the values returned in z are often the same as those
612: computed by RGComputeContour(), but this is not the case in real scalars
613: where all output arguments are real.
615: The computed values change for different quadrature rules (trapezoidal
616: or Chebyshev).
618: Level: intermediate
620: .seealso: RGComputeContour()
621: @*/
622: PetscErrorCode RGComputeQuadrature(RG rg,RGQuadRule quad,PetscInt n,PetscScalar z[],PetscScalar zn[],PetscScalar w[])
623: {
630: RGComputeContour(rg,n,z,NULL);
631: PetscUseTypeMethod(rg,computequadrature,quad,n,z,zn,w);
632: return 0;
633: }
635: /*@
636: RGSetComplement - Sets a flag to indicate that the region is the complement
637: of the specified one.
639: Logically Collective on rg
641: Input Parameters:
642: + rg - the region context
643: - flg - the boolean flag
645: Options Database Key:
646: . -rg_complement <bool> - Activate/deactivate the complementation of the region
648: Level: intermediate
650: .seealso: RGGetComplement()
651: @*/
652: PetscErrorCode RGSetComplement(RG rg,PetscBool flg)
653: {
656: rg->complement = flg;
657: return 0;
658: }
660: /*@
661: RGGetComplement - Gets a flag that that indicates whether the region
662: is complemented or not.
664: Not Collective
666: Input Parameter:
667: . rg - the region context
669: Output Parameter:
670: . flg - the flag
672: Level: intermediate
674: .seealso: RGSetComplement()
675: @*/
676: PetscErrorCode RGGetComplement(RG rg,PetscBool *flg)
677: {
680: *flg = rg->complement;
681: return 0;
682: }
684: /*@
685: RGSetScale - Sets the scaling factor to be used when checking that a
686: point is inside the region and when computing the contour.
688: Logically Collective on rg
690: Input Parameters:
691: + rg - the region context
692: - sfactor - the scaling factor
694: Options Database Key:
695: . -rg_scale <real> - Sets the scaling factor
697: Level: advanced
699: .seealso: RGGetScale(), RGCheckInside()
700: @*/
701: PetscErrorCode RGSetScale(RG rg,PetscReal sfactor)
702: {
705: if (sfactor == PETSC_DEFAULT || sfactor == PETSC_DECIDE) rg->sfactor = 1.0;
706: else {
708: rg->sfactor = sfactor;
709: }
710: return 0;
711: }
713: /*@
714: RGGetScale - Gets the scaling factor.
716: Not Collective
718: Input Parameter:
719: . rg - the region context
721: Output Parameter:
722: . sfactor - the scaling factor
724: Level: advanced
726: .seealso: RGSetScale()
727: @*/
728: PetscErrorCode RGGetScale(RG rg,PetscReal *sfactor)
729: {
732: *sfactor = rg->sfactor;
733: return 0;
734: }
736: /*@
737: RGPushScale - Sets an additional scaling factor, that will multiply the
738: user-defined scaling factor.
740: Logically Collective on rg
742: Input Parameters:
743: + rg - the region context
744: - sfactor - the scaling factor
746: Notes:
747: The current implementation does not allow pushing several scaling factors.
749: This is intended for internal use, for instance in polynomial eigensolvers
750: that use parameter scaling.
752: Level: developer
754: .seealso: RGPopScale(), RGSetScale()
755: @*/
756: PetscErrorCode RGPushScale(RG rg,PetscReal sfactor)
757: {
762: rg->osfactor = rg->sfactor;
763: rg->sfactor *= sfactor;
764: return 0;
765: }
767: /*@
768: RGPopScale - Pops the scaling factor set with RGPushScale().
770: Not Collective
772: Input Parameter:
773: . rg - the region context
775: Level: developer
777: .seealso: RGPushScale()
778: @*/
779: PetscErrorCode RGPopScale(RG rg)
780: {
783: rg->sfactor = rg->osfactor;
784: rg->osfactor = 0.0;
785: return 0;
786: }
788: /*@C
789: RGDestroy - Destroys RG context that was created with RGCreate().
791: Collective on rg
793: Input Parameter:
794: . rg - the region context
796: Level: beginner
798: .seealso: RGCreate()
799: @*/
800: PetscErrorCode RGDestroy(RG *rg)
801: {
802: if (!*rg) return 0;
804: if (--((PetscObject)(*rg))->refct > 0) { *rg = NULL; return 0; }
805: PetscTryTypeMethod(*rg,destroy);
806: PetscHeaderDestroy(rg);
807: return 0;
808: }
810: /*@C
811: RGRegister - Adds a region to the RG package.
813: Not collective
815: Input Parameters:
816: + name - name of a new user-defined RG
817: - function - routine to create context
819: Notes:
820: RGRegister() may be called multiple times to add several user-defined regions.
822: Level: advanced
824: .seealso: RGRegisterAll()
825: @*/
826: PetscErrorCode RGRegister(const char *name,PetscErrorCode (*function)(RG))
827: {
828: RGInitializePackage();
829: PetscFunctionListAdd(&RGList,name,function);
830: return 0;
831: }