Actual source code: rginterval.c
slepc-3.22.1 2024-10-28
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;
27: PetscFunctionBegin;
28: PetscCheck(a || b || c || d,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"At least one argument must be nonzero");
29: PetscCheck(a!=b || !a,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, endpoints must be distinct (or both zero)");
30: PetscCheck(a<=b,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be a<b");
31: PetscCheck(c!=d || !c,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, endpoints must be distinct (or both zero)");
32: PetscCheck(c<=d,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be c<d");
33: #if !defined(PETSC_USE_COMPLEX)
34: PetscCheck(c==-d,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"In real scalars the region must be symmetric wrt real axis");
35: #endif
36: ctx->a = a;
37: ctx->b = b;
38: ctx->c = c;
39: ctx->d = d;
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: /*@
44: RGIntervalSetEndpoints - Sets the parameters defining the interval region.
46: Logically Collective
48: Input Parameters:
49: + rg - the region context
50: . a - left endpoint of the interval in the real axis
51: . b - right endpoint of the interval in the real axis
52: . c - bottom endpoint of the interval in the imaginary axis
53: - d - top endpoint of the interval in the imaginary axis
55: Options Database Keys:
56: . -rg_interval_endpoints - the four endpoints
58: Note:
59: The region is defined as [a,b]x[c,d]. Particular cases are an interval on
60: the real axis (c=d=0), similar for the imaginary axis (a=b=0), the whole
61: complex plane (a=-inf,b=inf,c=-inf,d=inf), and so on.
63: When PETSc is built with real scalars, the region must be symmetric with
64: respect to the real axis.
66: Level: advanced
68: .seealso: RGIntervalGetEndpoints()
69: @*/
70: PetscErrorCode RGIntervalSetEndpoints(RG rg,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
71: {
72: PetscFunctionBegin;
78: PetscTryMethod(rg,"RGIntervalSetEndpoints_C",(RG,PetscReal,PetscReal,PetscReal,PetscReal),(rg,a,b,c,d));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: static PetscErrorCode RGIntervalGetEndpoints_Interval(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
83: {
84: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
86: PetscFunctionBegin;
87: if (a) *a = ctx->a;
88: if (b) *b = ctx->b;
89: if (c) *c = ctx->c;
90: if (d) *d = ctx->d;
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /*@
95: RGIntervalGetEndpoints - Gets the parameters that define the interval region.
97: Not Collective
99: Input Parameter:
100: . rg - the region context
102: Output Parameters:
103: + a - left endpoint of the interval in the real axis
104: . b - right endpoint of the interval in the real axis
105: . c - bottom endpoint of the interval in the imaginary axis
106: - d - top endpoint of the interval in the imaginary axis
108: Level: advanced
110: .seealso: RGIntervalSetEndpoints()
111: @*/
112: PetscErrorCode RGIntervalGetEndpoints(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
113: {
114: PetscFunctionBegin;
116: PetscUseMethod(rg,"RGIntervalGetEndpoints_C",(RG,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(rg,a,b,c,d));
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: static PetscErrorCode RGView_Interval(RG rg,PetscViewer viewer)
121: {
122: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
123: PetscBool isdraw,isascii;
124: int winw,winh;
125: PetscDraw draw;
126: PetscDrawAxis axis;
127: PetscReal a,b,c,d,ab,cd,lx,ly,w,scale=1.2;
129: PetscFunctionBegin;
130: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
131: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
132: if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer," region: [%g,%g]x[%g,%g]\n",RGShowReal(ctx->a),RGShowReal(ctx->b),RGShowReal(ctx->c),RGShowReal(ctx->d)));
133: else if (isdraw) {
134: PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
135: PetscCall(PetscDrawCheckResizedWindow(draw));
136: PetscCall(PetscDrawGetWindowSize(draw,&winw,&winh));
137: winw = PetscMax(winw,1); winh = PetscMax(winh,1);
138: PetscCall(PetscDrawClear(draw));
139: PetscCall(PetscDrawSetTitle(draw,"Interval region"));
140: PetscCall(PetscDrawAxisCreate(draw,&axis));
141: a = ctx->a*rg->sfactor;
142: b = ctx->b*rg->sfactor;
143: c = ctx->c*rg->sfactor;
144: d = ctx->d*rg->sfactor;
145: lx = b-a;
146: ly = d-c;
147: ab = (a+b)/2;
148: cd = (c+d)/2;
149: w = scale*PetscMax(lx/winw,ly/winh)/2;
150: PetscCall(PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh));
151: PetscCall(PetscDrawAxisDraw(axis));
152: PetscCall(PetscDrawAxisDestroy(&axis));
153: PetscCall(PetscDrawRectangle(draw,a,c,b,d,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE));
154: PetscCall(PetscDrawFlush(draw));
155: PetscCall(PetscDrawSave(draw));
156: PetscCall(PetscDrawPause(draw));
157: }
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: static PetscErrorCode RGIsTrivial_Interval(RG rg,PetscBool *trivial)
162: {
163: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
165: PetscFunctionBegin;
166: if (rg->complement) *trivial = (ctx->a==ctx->b && ctx->c==ctx->d)? PETSC_TRUE: PETSC_FALSE;
167: 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;
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: static PetscErrorCode RGComputeContour_Interval(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
172: {
173: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
174: PetscInt i,N,Nv,Nh,k1,k0;
175: PetscReal hv,hh,t;
177: PetscFunctionBegin;
178: PetscCheck(ctx->a>-PETSC_MAX_REAL && ctx->b<PETSC_MAX_REAL && ctx->c>-PETSC_MAX_REAL && ctx->d<PETSC_MAX_REAL,PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"Contour not defined in unbounded regions");
179: if (ctx->a==ctx->b || ctx->c==ctx->d) {
180: if (ctx->a==ctx->b) {hv = (ctx->d-ctx->c)/(n-1); hh = 0.0;}
181: else {hh = (ctx->b-ctx->a)/(n-1); hv = 0.0;}
182: for (i=0;i<n;i++) {
183: #if defined(PETSC_USE_COMPLEX)
184: cr[i] = PetscCMPLX(ctx->a+hh*i,ctx->c+hv*i);
185: #else
186: if (cr) cr[i] = ctx->a+hh*i;
187: if (ci) ci[i] = ctx->c+hv*i;
188: #endif
189: }
190: } else {
191: PetscCheck(n>3,PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"Minimum number of contour points: 4");
192: N = n/2;
193: t = ((ctx->d-ctx->c)/(ctx->d-ctx->c+ctx->b-ctx->a))*N;
194: Nv = t-PetscFloorReal(t)>0.5?PetscCeilReal(t):PetscFloorReal(t);
195: if (Nv==0) Nv++;
196: else if (Nv==N) Nv--;
197: Nh = N-Nv;
198: hh = (ctx->b-ctx->a)/Nh;
199: hv = (ctx->d-ctx->c)/Nv;
200: /* positive imaginary part first */
201: k1 = Nv/2+1;
202: k0 = Nv-k1;
204: for (i=k1;i<Nv;i++) {
205: #if defined(PETSC_USE_COMPLEX)
206: cr[i-k1] = PetscCMPLX(ctx->b,ctx->c+i*hv);
207: cr[i-k1+N] = PetscCMPLX(ctx->a,ctx->d-i*hv);
208: #else
209: if (cr) {cr[i-k1] = ctx->b; cr[i-k1+N] = ctx->a;}
210: if (ci) {ci[i-k1] = ctx->c+i*hv; ci[i-k1+N] = ctx->d-i*hv;}
211: #endif
212: }
213: for (i=0;i<Nh;i++) {
214: #if defined(PETSC_USE_COMPLEX)
215: cr[i+k0] = PetscCMPLX(ctx->b-i*hh,ctx->d);
216: cr[i+k0+N] = PetscCMPLX(ctx->a+i*hh,ctx->c);
217: #else
218: if (cr) {cr[i+k0] = ctx->b-i*hh; cr[i+k0+N] = ctx->a+i*hh;}
219: if (ci) {ci[i+k0] = ctx->d; ci[i+k0+N] = ctx->c;}
220: #endif
221: }
222: for (i=0;i<k1;i++) {
223: #if defined(PETSC_USE_COMPLEX)
224: cr[i+k0+Nh] = PetscCMPLX(ctx->a,ctx->d-i*hv);
225: cr[i+k0+Nh+N] = PetscCMPLX(ctx->b,ctx->c+i*hv);
226: #else
227: if (cr) {cr[i+k0+Nh] = ctx->a; cr[i+k0+Nh+N] = ctx->b;}
228: if (ci) {ci[i+k0+Nh] = ctx->d+i*hv; ci[i+k0+Nh+N] = ctx->c-i*hv;}
229: #endif
230: }
231: }
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: static PetscErrorCode RGComputeBoundingBox_Interval(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
236: {
237: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
239: PetscFunctionBegin;
240: if (a) *a = ctx->a;
241: if (b) *b = ctx->b;
242: if (c) *c = ctx->c;
243: if (d) *d = ctx->d;
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: static PetscErrorCode RGComputeQuadrature_Interval(RG rg,RGQuadRule quad,PetscInt n,PetscScalar *z,PetscScalar *zn,PetscScalar *w)
248: {
249: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
250: PetscReal theta,max_w=0.0,radius=1.0;
251: PetscScalar tmp,tmp2,center=0.0;
252: PetscInt i,j;
254: PetscFunctionBegin;
255: if (quad == RG_QUADRULE_CHEBYSHEV) {
256: PetscCheck((ctx->c==ctx->d && ctx->c==0.0) || (ctx->a==ctx->b && ctx->a==0.0),PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"Endpoints of the imaginary axis or the real axis must be both zero");
257: for (i=0;i<n;i++) {
258: theta = PETSC_PI*(i+0.5)/n;
259: zn[i] = PetscCosReal(theta);
260: w[i] = PetscCosReal((n-1)*theta)/n;
261: if (ctx->c==ctx->d) z[i] = ((ctx->b-ctx->a)*(zn[i]+1.0)/2.0+ctx->a)*rg->sfactor;
262: else if (ctx->a==ctx->b) {
263: #if defined(PETSC_USE_COMPLEX)
264: z[i] = ((ctx->d-ctx->c)*(zn[i]+1.0)/2.0+ctx->c)*rg->sfactor*PETSC_i;
265: #else
266: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Integration points on a vertical line require complex arithmetic");
267: #endif
268: }
269: }
270: } else { /* RG_QUADRULE_TRAPEZOIDAL */
271: #if defined(PETSC_USE_COMPLEX)
272: center = rg->sfactor*PetscCMPLX(ctx->b+ctx->a,ctx->d+ctx->c)/2.0;
273: #else
274: center = rg->sfactor*(ctx->b+ctx->a)/2.0;
275: #endif
276: radius = PetscSqrtReal(PetscPowRealInt(rg->sfactor*(ctx->b-ctx->a)/2.0,2)+PetscPowRealInt(rg->sfactor*(ctx->d-ctx->c)/2.0,2));
277: for (i=0;i<n;i++) {
278: zn[i] = (z[i]-center)/radius;
279: tmp = 1.0; tmp2 = 1.0;
280: for (j=0;j<n;j++) {
281: tmp *= z[j];
282: if (i != j) tmp2 *= z[j]-z[i];
283: }
284: w[i] = tmp/tmp2;
285: max_w = PetscMax(PetscAbsScalar(w[i]),max_w);
286: }
287: for (i=0;i<n;i++) w[i] /= (PetscScalar)max_w;
288: }
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: static PetscErrorCode RGCheckInside_Interval(RG rg,PetscReal dx,PetscReal dy,PetscInt *inside)
293: {
294: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
296: PetscFunctionBegin;
297: if (dx>ctx->a && dx<ctx->b) *inside = 1;
298: else if (dx==ctx->a || dx==ctx->b) *inside = 0;
299: else *inside = -1;
300: if (*inside>=0) {
301: if (dy>ctx->c && dy<ctx->d) ;
302: else if (dy==ctx->c || dy==ctx->d) *inside = 0;
303: else *inside = -1;
304: }
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: static PetscErrorCode RGIsAxisymmetric_Interval(RG rg,PetscBool vertical,PetscBool *symm)
309: {
310: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
312: PetscFunctionBegin;
313: if (vertical) *symm = (ctx->a == -ctx->b)? PETSC_TRUE: PETSC_FALSE;
314: else *symm = (ctx->c == -ctx->d)? PETSC_TRUE: PETSC_FALSE;
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: static PetscErrorCode RGSetFromOptions_Interval(RG rg,PetscOptionItems *PetscOptionsObject)
319: {
320: PetscBool flg;
321: PetscInt k;
322: PetscReal array[4]={0,0,0,0};
324: PetscFunctionBegin;
325: PetscOptionsHeadBegin(PetscOptionsObject,"RG Interval Options");
327: k = 4;
328: PetscCall(PetscOptionsRealArray("-rg_interval_endpoints","Interval endpoints (two or four real values separated with a comma without spaces)","RGIntervalSetEndpoints",array,&k,&flg));
329: if (flg) {
330: PetscCheck(k>1,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_SIZ,"Must pass at least two values in -rg_interval_endpoints (comma-separated without spaces)");
331: PetscCall(RGIntervalSetEndpoints(rg,array[0],array[1],array[2],array[3]));
332: }
334: PetscOptionsHeadEnd();
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }
338: static PetscErrorCode RGDestroy_Interval(RG rg)
339: {
340: PetscFunctionBegin;
341: PetscCall(PetscFree(rg->data));
342: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",NULL));
343: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",NULL));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: SLEPC_EXTERN PetscErrorCode RGCreate_Interval(RG rg)
348: {
349: RG_INTERVAL *interval;
351: PetscFunctionBegin;
352: PetscCall(PetscNew(&interval));
353: interval->a = -PETSC_MAX_REAL;
354: interval->b = PETSC_MAX_REAL;
355: interval->c = -PETSC_MAX_REAL;
356: interval->d = PETSC_MAX_REAL;
357: rg->data = (void*)interval;
359: rg->ops->istrivial = RGIsTrivial_Interval;
360: rg->ops->computecontour = RGComputeContour_Interval;
361: rg->ops->computebbox = RGComputeBoundingBox_Interval;
362: rg->ops->computequadrature = RGComputeQuadrature_Interval;
363: rg->ops->checkinside = RGCheckInside_Interval;
364: rg->ops->isaxisymmetric = RGIsAxisymmetric_Interval;
365: rg->ops->setfromoptions = RGSetFromOptions_Interval;
366: rg->ops->view = RGView_Interval;
367: rg->ops->destroy = RGDestroy_Interval;
368: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",RGIntervalSetEndpoints_Interval));
369: PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",RGIntervalGetEndpoints_Interval));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }