Actual source code: rginterval.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }