Actual source code: rginterval.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcrg.h" I*/
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;
28: if (!a && !b && !c && !d) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"At least one argument must be nonzero");
29: if (a==b && a) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, endpoints must be distinct (or both zero)");
30: if (a>b) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be a<b");
31: if (c==d && c) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, endpoints must be distinct (or both zero)");
32: if (c>d) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be c<d");
33: #if !defined(PETSC_USE_COMPLEX)
34: if (c!=-d) SETERRQ(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: return(0);
41: }
43: /*@
44: RGIntervalSetEndpoints - Sets the parameters defining the interval region.
46: Logically Collective on RG
48: Input Parameters:
49: + rg - the region context
50: . a,b - endpoints of the interval in the real axis
51: - c,d - endpoints of the interval in the imaginary axis
53: Options Database Keys:
54: . -rg_interval_endpoints - the four endpoints
56: Note:
57: The region is defined as [a,b]x[c,d]. Particular cases are an interval on
58: the real axis (c=d=0), similar for the imaginary axis (a=b=0), the whole
59: complex plane (a=-inf,b=inf,c=-inf,d=inf), and so on.
61: When PETSc is built with real scalars, the region must be symmetric with
62: respect to the real axis.
64: Level: advanced
66: .seealso: RGIntervalGetEndpoints()
67: @*/
68: PetscErrorCode RGIntervalSetEndpoints(RG rg,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
69: {
78: PetscTryMethod(rg,"RGIntervalSetEndpoints_C",(RG,PetscReal,PetscReal,PetscReal,PetscReal),(rg,a,b,c,d));
79: return(0);
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;
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: return(0);
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,b - endpoints of the interval in the real axis
104: - c,d - endpoints of the interval in the imaginary axis
106: Level: advanced
108: .seealso: RGIntervalSetEndpoints()
109: @*/
110: PetscErrorCode RGIntervalGetEndpoints(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
111: {
116: PetscUseMethod(rg,"RGIntervalGetEndpoints_C",(RG,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(rg,a,b,c,d));
117: return(0);
118: }
120: PetscErrorCode RGView_Interval(RG rg,PetscViewer viewer)
121: {
123: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
124: PetscBool isdraw,isascii;
125: int winw,winh;
126: PetscDraw draw;
127: PetscDrawAxis axis;
128: PetscReal a,b,c,d,ab,cd,lx,ly,w,scale=1.2;
131: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
132: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
133: if (isascii) {
134: PetscViewerASCIIPrintf(viewer," region: [%g,%g]x[%g,%g]\n",RGShowReal(ctx->a),RGShowReal(ctx->b),RGShowReal(ctx->c),RGShowReal(ctx->d));
135: } else if (isdraw) {
136: PetscViewerDrawGetDraw(viewer,0,&draw);
137: PetscDrawCheckResizedWindow(draw);
138: PetscDrawGetWindowSize(draw,&winw,&winh);
139: winw = PetscMax(winw,1); winh = PetscMax(winh,1);
140: PetscDrawClear(draw);
141: PetscDrawSetTitle(draw,"Interval region");
142: PetscDrawAxisCreate(draw,&axis);
143: a = ctx->a*rg->sfactor;
144: b = ctx->b*rg->sfactor;
145: c = ctx->c*rg->sfactor;
146: d = ctx->d*rg->sfactor;
147: lx = b-a;
148: ly = d-c;
149: ab = (a+b)/2;
150: cd = (c+d)/2;
151: w = scale*PetscMax(lx/winw,ly/winh)/2;
152: PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh);
153: PetscDrawAxisDraw(axis);
154: PetscDrawAxisDestroy(&axis);
155: PetscDrawRectangle(draw,a,c,b,d,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE);
156: PetscDrawFlush(draw);
157: PetscDrawSave(draw);
158: PetscDrawPause(draw);
159: }
160: return(0);
161: }
163: PetscErrorCode RGIsTrivial_Interval(RG rg,PetscBool *trivial)
164: {
165: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
168: if (rg->complement) *trivial = (ctx->a==ctx->b && ctx->c==ctx->d)? PETSC_TRUE: PETSC_FALSE;
169: 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;
170: return(0);
171: }
173: PetscErrorCode RGComputeContour_Interval(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
174: {
175: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
176: PetscInt i,pt,idx,j;
177: PetscReal hr[4],hi[4],h,off,d[4];
178: PetscScalar vr[4],vi[4];
181: if (!(ctx->a>-PETSC_MAX_REAL && ctx->b<PETSC_MAX_REAL && ctx->c>-PETSC_MAX_REAL && ctx->d<PETSC_MAX_REAL)) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"Contour not defined in unbounded regions");
182: if (ctx->a==ctx->b || ctx->c==ctx->d) {
183: if (ctx->a==ctx->b) {hi[0] = (ctx->d-ctx->c)/(n-1); hr[0] = 0.0;}
184: else {hr[0] = (ctx->b-ctx->a)/(n-1); hi[0] = 0.0;}
185: for (i=0;i<n;i++) {
186: #if defined(PETSC_USE_COMPLEX)
187: cr[i] = ctx->a+hr[0]*i + (ctx->c+hi[0]*i)*PETSC_i;
188: #else
189: cr[i] = ctx->a+hr[0]*i; ci[i] = ctx->c+hi[0]*i;
190: #endif
191: }
192: } else {
193: d[1] = d[3] = ctx->d-ctx->c; d[0] = d[2] = ctx->b-ctx->a;
194: h = (2*(d[0]+d[1]))/n;
195: vr[0] = ctx->a; vr[1] = ctx->b; vr[2] = ctx->b; vr[3] = ctx->a;
196: vi[0] = ctx->c; vi[1] = ctx->c; vi[2] = ctx->d; vi[3] = ctx->d;
197: hr[0] = h; hr[1] = 0.0; hr[2] = -h; hr[3] = 0.0;
198: hi[0] = 0.0; hi[1] = h; hi[2] = 0.0; hi[3] = -h;
199: off = 0.0; idx = 0;
200: for (i=0;i<4;i++) {
201: #if defined(PETSC_USE_COMPLEX)
202: cr[idx] = vr[i]+off*(hr[i]/h)+ (vi[i]+off*(hi[i]/h))*PETSC_i;
203: #else
204: cr[idx] = vr[i]+off*(hr[i]/h); ci[idx]=vi[i]+off*(hi[i]/h);
205: #endif
206: idx++;
207: pt = (d[i]-off)/h+1;
208: for (j=1;j<pt && idx<n;j++) {
209: #if defined(PETSC_USE_COMPLEX)
210: cr[idx] = cr[idx-1]+(hr[i]+hi[i]*PETSC_i);
211: #else
212: cr[idx] = cr[idx-1]+hr[i]; ci[idx] = ci[idx-1]+hi[i];
213: #endif
214: idx++;
215: }
216: off = off+pt*h-d[i];
217: if (off>=d[i+1]) {off -= d[i+1]; i++;}
218: }
219: }
220: return(0);
221: }
223: PetscErrorCode RGComputeBoundingBox_Interval(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
224: {
225: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
228: if (a) *a = ctx->a;
229: if (b) *b = ctx->b;
230: if (c) *c = ctx->c;
231: if (d) *d = ctx->d;
232: return(0);
233: }
235: PetscErrorCode RGCheckInside_Interval(RG rg,PetscReal dx,PetscReal dy,PetscInt *inside)
236: {
237: RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
240: if (dx>ctx->a && dx<ctx->b) *inside = 1;
241: else if (dx==ctx->a || dx==ctx->b) *inside = 0;
242: else *inside = -1;
243: if (*inside>=0) {
244: if (dy>ctx->c && dy<ctx->d) ;
245: else if (dy==ctx->c || dy==ctx->d) *inside = 0;
246: else *inside = -1;
247: }
248: return(0);
249: }
251: PetscErrorCode RGSetFromOptions_Interval(PetscOptionItems *PetscOptionsObject,RG rg)
252: {
254: PetscBool flg;
255: PetscInt k;
256: PetscReal array[4]={0,0,0,0};
259: PetscOptionsHead(PetscOptionsObject,"RG Interval Options");
261: k = 4;
262: PetscOptionsRealArray("-rg_interval_endpoints","Interval endpoints (two or four real values separated with a comma without spaces)","RGIntervalSetEndpoints",array,&k,&flg);
263: if (flg) {
264: if (k<2) SETERRQ(PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_SIZ,"Must pass at least two values in -rg_interval_endpoints (comma-separated without spaces)");
265: RGIntervalSetEndpoints(rg,array[0],array[1],array[2],array[3]);
266: }
268: PetscOptionsTail();
269: return(0);
270: }
272: PetscErrorCode RGDestroy_Interval(RG rg)
273: {
277: PetscFree(rg->data);
278: PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",NULL);
279: PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",NULL);
280: return(0);
281: }
283: SLEPC_EXTERN PetscErrorCode RGCreate_Interval(RG rg)
284: {
285: RG_INTERVAL *interval;
289: PetscNewLog(rg,&interval);
290: interval->a = -PETSC_MAX_REAL;
291: interval->b = PETSC_MAX_REAL;
292: interval->c = -PETSC_MAX_REAL;
293: interval->d = PETSC_MAX_REAL;
294: rg->data = (void*)interval;
296: rg->ops->istrivial = RGIsTrivial_Interval;
297: rg->ops->computecontour = RGComputeContour_Interval;
298: rg->ops->computebbox = RGComputeBoundingBox_Interval;
299: rg->ops->checkinside = RGCheckInside_Interval;
300: rg->ops->setfromoptions = RGSetFromOptions_Interval;
301: rg->ops->view = RGView_Interval;
302: rg->ops->destroy = RGDestroy_Interval;
303: PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",RGIntervalSetEndpoints_Interval);
304: PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",RGIntervalGetEndpoints_Interval);
305: return(0);
306: }