Actual source code: rginterval.c
slepc-3.18.0 2022-10-01
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: Interval of the real axis or more generally a (possibly open) rectangle
12: of the complex plane
13: */
15: #include <slepc/private/rgimpl.h>
16: #include <petscdraw.h>
18: typedef struct {
19: PetscReal a,b; /* interval in the real axis */
20: PetscReal c,d; /* interval in the imaginary axis */
21: } RG_INTERVAL;
23: static PetscErrorCode RGIntervalSetEndpoints_Interval(RG rg,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
24: {
25: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
32: #if !defined(PETSC_USE_COMPLEX)
34: #endif
35: ctx->a = a;
36: ctx->b = b;
37: ctx->c = c;
38: ctx->d = d;
39: return 0;
40: }
42: /*@
43: RGIntervalSetEndpoints - Sets the parameters defining the interval region.
45: Logically Collective on rg
47: Input Parameters:
48: + rg - the region context
49: . a - left endpoint of the interval in the real axis
50: . b - right endpoint of the interval in the real axis
51: . c - bottom endpoint of the interval in the imaginary axis
52: - d - top endpoint of the interval in the imaginary axis
54: Options Database Keys:
55: . -rg_interval_endpoints - the four endpoints
57: Note:
58: The region is defined as [a,b]x[c,d]. Particular cases are an interval on
59: the real axis (c=d=0), similar for the imaginary axis (a=b=0), the whole
60: complex plane (a=-inf,b=inf,c=-inf,d=inf), and so on.
62: When PETSc is built with real scalars, the region must be symmetric with
63: respect to the real axis.
65: Level: advanced
67: .seealso: RGIntervalGetEndpoints()
68: @*/
69: PetscErrorCode RGIntervalSetEndpoints(RG rg,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
70: {
76: PetscTryMethod(rg,"RGIntervalSetEndpoints_C",(RG,PetscReal,PetscReal,PetscReal,PetscReal),(rg,a,b,c,d));
77: return 0;
78: }
80: static PetscErrorCode RGIntervalGetEndpoints_Interval(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
81: {
82: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
84: if (a) *a = ctx->a;
85: if (b) *b = ctx->b;
86: if (c) *c = ctx->c;
87: if (d) *d = ctx->d;
88: return 0;
89: }
91: /*@C
92: RGIntervalGetEndpoints - Gets the parameters that define the interval region.
94: Not Collective
96: Input Parameter:
97: . rg - the region context
99: Output Parameters:
100: + a - left endpoint of the interval in the real axis
101: . b - right endpoint of the interval in the real axis
102: . c - bottom endpoint of the interval in the imaginary axis
103: - d - top endpoint of the interval in the imaginary axis
105: Level: advanced
107: .seealso: RGIntervalSetEndpoints()
108: @*/
109: PetscErrorCode RGIntervalGetEndpoints(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
110: {
112: PetscUseMethod(rg,"RGIntervalGetEndpoints_C",(RG,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(rg,a,b,c,d));
113: return 0;
114: }
116: PetscErrorCode RGView_Interval(RG rg,PetscViewer viewer)
117: {
118: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
119: PetscBool isdraw,isascii;
120: int winw,winh;
121: PetscDraw draw;
122: PetscDrawAxis axis;
123: PetscReal a,b,c,d,ab,cd,lx,ly,w,scale=1.2;
125: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
126: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
127: if (isascii) PetscViewerASCIIPrintf(viewer," region: [%g,%g]x[%g,%g]\n",RGShowReal(ctx->a),RGShowReal(ctx->b),RGShowReal(ctx->c),RGShowReal(ctx->d));
128: else if (isdraw) {
129: PetscViewerDrawGetDraw(viewer,0,&draw);
130: PetscDrawCheckResizedWindow(draw);
131: PetscDrawGetWindowSize(draw,&winw,&winh);
132: winw = PetscMax(winw,1); winh = PetscMax(winh,1);
133: PetscDrawClear(draw);
134: PetscDrawSetTitle(draw,"Interval region");
135: PetscDrawAxisCreate(draw,&axis);
136: a = ctx->a*rg->sfactor;
137: b = ctx->b*rg->sfactor;
138: c = ctx->c*rg->sfactor;
139: d = ctx->d*rg->sfactor;
140: lx = b-a;
141: ly = d-c;
142: ab = (a+b)/2;
143: cd = (c+d)/2;
144: w = scale*PetscMax(lx/winw,ly/winh)/2;
145: PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh);
146: PetscDrawAxisDraw(axis);
147: PetscDrawAxisDestroy(&axis);
148: PetscDrawRectangle(draw,a,c,b,d,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE);
149: PetscDrawFlush(draw);
150: PetscDrawSave(draw);
151: PetscDrawPause(draw);
152: }
153: return 0;
154: }
156: PetscErrorCode RGIsTrivial_Interval(RG rg,PetscBool *trivial)
157: {
158: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
160: if (rg->complement) *trivial = (ctx->a==ctx->b && ctx->c==ctx->d)? PETSC_TRUE: PETSC_FALSE;
161: else *trivial = (ctx->a<=-PETSC_MAX_REAL && ctx->b>=PETSC_MAX_REAL && ctx->c<=-PETSC_MAX_REAL && ctx->d>=PETSC_MAX_REAL)? PETSC_TRUE: PETSC_FALSE;
162: return 0;
163: }
165: PetscErrorCode RGComputeContour_Interval(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
166: {
167: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
168: PetscInt i,N,Nv,Nh,k1,k0;
169: PetscReal hv,hh,t;
172: if (ctx->a==ctx->b || ctx->c==ctx->d) {
173: if (ctx->a==ctx->b) {hv = (ctx->d-ctx->c)/(n-1); hh = 0.0;}
174: else {hh = (ctx->b-ctx->a)/(n-1); hv = 0.0;}
175: for (i=0;i<n;i++) {
176: #if defined(PETSC_USE_COMPLEX)
177: cr[i] = PetscCMPLX(ctx->a+hh*i,ctx->c+hv*i);
178: #else
179: if (cr) cr[i] = ctx->a+hh*i;
180: if (ci) ci[i] = ctx->c+hv*i;
181: #endif
182: }
183: } else {
185: N = n/2;
186: t = ((ctx->d-ctx->c)/(ctx->d-ctx->c+ctx->b-ctx->a))*N;
187: Nv = t-PetscFloorReal(t)>0.5?PetscCeilReal(t):PetscFloorReal(t);
188: if (Nv==0) Nv++;
189: else if (Nv==N) Nv--;
190: Nh = N-Nv;
191: hh = (ctx->b-ctx->a)/Nh;
192: hv = (ctx->d-ctx->c)/Nv;
193: /* positive imaginary part first */
194: k1 = Nv/2+1;
195: k0 = Nv-k1;
197: for (i=k1;i<Nv;i++) {
198: #if defined(PETSC_USE_COMPLEX)
199: cr[i-k1] = PetscCMPLX(ctx->b,ctx->c+i*hv);
200: cr[i-k1+N] = PetscCMPLX(ctx->a,ctx->d-i*hv);
201: #else
202: if (cr) {cr[i-k1] = ctx->b; cr[i-k1+N] = ctx->a;}
203: if (ci) {ci[i-k1] = ctx->c+i*hv; ci[i-k1+N] = ctx->d-i*hv;}
204: #endif
205: }
206: for (i=0;i<Nh;i++) {
207: #if defined(PETSC_USE_COMPLEX)
208: cr[i+k0] = PetscCMPLX(ctx->b-i*hh,ctx->d);
209: cr[i+k0+N] = PetscCMPLX(ctx->a+i*hh,ctx->c);
210: #else
211: if (cr) {cr[i+k0] = ctx->b-i*hh; cr[i+k0+N] = ctx->a+i*hh;}
212: if (ci) {ci[i+k0] = ctx->d; ci[i+k0+N] = ctx->c;}
213: #endif
214: }
215: for (i=0;i<k1;i++) {
216: #if defined(PETSC_USE_COMPLEX)
217: cr[i+k0+Nh] = PetscCMPLX(ctx->a,ctx->d-i*hv);
218: cr[i+k0+Nh+N] = PetscCMPLX(ctx->b,ctx->c+i*hv);
219: #else
220: if (cr) {cr[i+k0+Nh] = ctx->a; cr[i+k0+Nh+N] = ctx->b;}
221: if (ci) {ci[i+k0+Nh] = ctx->d+i*hv; ci[i+k0+Nh+N] = ctx->c-i*hv;}
222: #endif
223: }
224: }
225: return 0;
226: }
228: PetscErrorCode RGComputeBoundingBox_Interval(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
229: {
230: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
232: if (a) *a = ctx->a;
233: if (b) *b = ctx->b;
234: if (c) *c = ctx->c;
235: if (d) *d = ctx->d;
236: return 0;
237: }
239: PetscErrorCode RGComputeQuadrature_Interval(RG rg,RGQuadRule quad,PetscInt n,PetscScalar *z,PetscScalar *zn,PetscScalar *w)
240: {
241: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
242: PetscReal theta,max_w=0.0,radius=1.0;
243: PetscScalar tmp,tmp2,center=0.0;
244: PetscInt i,j;
246: if (quad == RG_QUADRULE_CHEBYSHEV) {
248: for (i=0;i<n;i++) {
249: theta = PETSC_PI*(i+0.5)/n;
250: zn[i] = PetscCosReal(theta);
251: w[i] = PetscCosReal((n-1)*theta)/n;
252: if (ctx->c==ctx->d) z[i] = ((ctx->b-ctx->a)*(zn[i]+1.0)/2.0+ctx->a)*rg->sfactor;
253: else if (ctx->a==ctx->b) {
254: #if defined(PETSC_USE_COMPLEX)
255: z[i] = ((ctx->d-ctx->c)*(zn[i]+1.0)/2.0+ctx->c)*rg->sfactor*PETSC_i;
256: #else
257: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
258: #endif
259: }
260: }
261: } else { /* RG_QUADRULE_TRAPEZOIDAL */
262: #if defined(PETSC_USE_COMPLEX)
263: center = rg->sfactor*PetscCMPLX(ctx->b+ctx->a,ctx->d+ctx->c)/2.0;
264: #else
265: center = rg->sfactor*(ctx->b+ctx->a)/2.0;
266: #endif
267: radius = PetscSqrtReal(PetscPowRealInt(rg->sfactor*(ctx->b-ctx->a)/2.0,2)+PetscPowRealInt(rg->sfactor*(ctx->d-ctx->c)/2.0,2));
268: for (i=0;i<n;i++) {
269: zn[i] = (z[i]-center)/radius;
270: tmp = 1.0; tmp2 = 1.0;
271: for (j=0;j<n;j++) {
272: tmp *= z[j];
273: if (i != j) tmp2 *= z[j]-z[i];
274: }
275: w[i] = tmp/tmp2;
276: max_w = PetscMax(PetscAbsScalar(w[i]),max_w);
277: }
278: for (i=0;i<n;i++) w[i] /= (PetscScalar)max_w;
279: }
280: return 0;
281: }
283: PetscErrorCode RGCheckInside_Interval(RG rg,PetscReal dx,PetscReal dy,PetscInt *inside)
284: {
285: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
287: if (dx>ctx->a && dx<ctx->b) *inside = 1;
288: else if (dx==ctx->a || dx==ctx->b) *inside = 0;
289: else *inside = -1;
290: if (*inside>=0) {
291: if (dy>ctx->c && dy<ctx->d) ;
292: else if (dy==ctx->c || dy==ctx->d) *inside = 0;
293: else *inside = -1;
294: }
295: return 0;
296: }
298: PetscErrorCode RGIsAxisymmetric_Interval(RG rg,PetscBool vertical,PetscBool *symm)
299: {
300: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
302: if (vertical) *symm = (ctx->a == -ctx->b)? PETSC_TRUE: PETSC_FALSE;
303: else *symm = (ctx->c == -ctx->d)? PETSC_TRUE: PETSC_FALSE;
304: return 0;
305: }
307: PetscErrorCode RGSetFromOptions_Interval(RG rg,PetscOptionItems *PetscOptionsObject)
308: {
309: PetscBool flg;
310: PetscInt k;
311: PetscReal array[4]={0,0,0,0};
313: PetscOptionsHeadBegin(PetscOptionsObject,"RG Interval Options");
315: k = 4;
316: PetscOptionsRealArray("-rg_interval_endpoints","Interval endpoints (two or four real values separated with a comma without spaces)","RGIntervalSetEndpoints",array,&k,&flg);
317: if (flg) {
319: RGIntervalSetEndpoints(rg,array[0],array[1],array[2],array[3]);
320: }
322: PetscOptionsHeadEnd();
323: return 0;
324: }
326: PetscErrorCode RGDestroy_Interval(RG rg)
327: {
328: PetscFree(rg->data);
329: PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",NULL);
330: PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",NULL);
331: return 0;
332: }
334: SLEPC_EXTERN PetscErrorCode RGCreate_Interval(RG rg)
335: {
336: RG_INTERVAL *interval;
338: PetscNew(&interval);
339: interval->a = -PETSC_MAX_REAL;
340: interval->b = PETSC_MAX_REAL;
341: interval->c = -PETSC_MAX_REAL;
342: interval->d = PETSC_MAX_REAL;
343: rg->data = (void*)interval;
345: rg->ops->istrivial = RGIsTrivial_Interval;
346: rg->ops->computecontour = RGComputeContour_Interval;
347: rg->ops->computebbox = RGComputeBoundingBox_Interval;
348: rg->ops->computequadrature = RGComputeQuadrature_Interval;
349: rg->ops->checkinside = RGCheckInside_Interval;
350: rg->ops->isaxisymmetric = RGIsAxisymmetric_Interval;
351: rg->ops->setfromoptions = RGSetFromOptions_Interval;
352: rg->ops->view = RGView_Interval;
353: rg->ops->destroy = RGDestroy_Interval;
354: PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",RGIntervalSetEndpoints_Interval);
355: PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",RGIntervalGetEndpoints_Interval);
356: return 0;
357: }