Actual source code: rginterval.c

slepc-3.14.0 2020-09-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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;

 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: /*@C
 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],vr[4],vi[4];

180:   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");
181:   if (ctx->a==ctx->b || ctx->c==ctx->d) {
182:     if (ctx->a==ctx->b) {hi[0] = (ctx->d-ctx->c)/(n-1); hr[0] = 0.0;}
183:     else {hr[0] = (ctx->b-ctx->a)/(n-1); hi[0] = 0.0;}
184:     for (i=0;i<n;i++) {
185: #if defined(PETSC_USE_COMPLEX)
186:       cr[i] = PetscCMPLX(ctx->a+hr[0]*i,ctx->c+hi[0]*i);
187: #else
188:       cr[i] = ctx->a+hr[0]*i; ci[i] = ctx->c+hi[0]*i;
189: #endif
190:     }
191:   } else {
192:     d[1] = d[3] = ctx->d-ctx->c; d[0] = d[2] = ctx->b-ctx->a;
193:     h = 2.0*(d[0]+d[1])/n;
194:     vr[0] = ctx->a; vr[1] = ctx->b; vr[2] = ctx->b; vr[3] = ctx->a;
195:     vi[0] = ctx->c; vi[1] = ctx->c; vi[2] = ctx->d; vi[3] = ctx->d;
196:     hr[0] = h;   hr[1] = 0.0; hr[2] = -h;  hr[3] = 0.0;
197:     hi[0] = 0.0; hi[1] = h;   hi[2] = 0.0; hi[3] = -h;
198:     off = 0.0; idx = 0;
199:     for (i=0;i<4;i++) {
200: #if defined(PETSC_USE_COMPLEX)
201:       cr[idx] = PetscCMPLX(vr[i]+off*(hr[i]/h),vi[i]+off*(hi[i]/h));
202: #else
203:       cr[idx] = vr[i]+off*(hr[i]/h); ci[idx] = vi[i]+off*(hi[i]/h);
204: #endif
205:       idx++;
206:       pt = (PetscInt)((d[i]-off)/h)+1;
207:       for (j=1;j<pt && idx<n;j++) {
208: #if defined(PETSC_USE_COMPLEX)
209:         cr[idx] = cr[idx-1]+PetscCMPLX(hr[i],hi[i]);
210: #else
211:         cr[idx] = cr[idx-1]+hr[i]; ci[idx] = ci[idx-1]+hi[i];
212: #endif
213:         idx++;
214:       }
215:       off += pt*h-d[i];
216:       if (off>=d[i+1]) {off -= d[i+1]; i++;}
217:     }
218:   }
219:   return(0);
220: }

222: PetscErrorCode RGComputeBoundingBox_Interval(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
223: {
224:   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;

227:   if (a) *a = ctx->a;
228:   if (b) *b = ctx->b;
229:   if (c) *c = ctx->c;
230:   if (d) *d = ctx->d;
231:   return(0);
232: }

234: PetscErrorCode RGCheckInside_Interval(RG rg,PetscReal dx,PetscReal dy,PetscInt *inside)
235: {
236:   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;

239:   if (dx>ctx->a && dx<ctx->b) *inside = 1;
240:   else if (dx==ctx->a || dx==ctx->b) *inside = 0;
241:   else *inside = -1;
242:   if (*inside>=0) {
243:     if (dy>ctx->c && dy<ctx->d) ;
244:     else if (dy==ctx->c || dy==ctx->d) *inside = 0;
245:     else *inside = -1;
246:   }
247:   return(0);
248: }

250: PetscErrorCode RGSetFromOptions_Interval(PetscOptionItems *PetscOptionsObject,RG rg)
251: {
253:   PetscBool      flg;
254:   PetscInt       k;
255:   PetscReal      array[4]={0,0,0,0};

258:   PetscOptionsHead(PetscOptionsObject,"RG Interval Options");

260:     k = 4;
261:     PetscOptionsRealArray("-rg_interval_endpoints","Interval endpoints (two or four real values separated with a comma without spaces)","RGIntervalSetEndpoints",array,&k,&flg);
262:     if (flg) {
263:       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)");
264:       RGIntervalSetEndpoints(rg,array[0],array[1],array[2],array[3]);
265:     }

267:   PetscOptionsTail();
268:   return(0);
269: }

271: PetscErrorCode RGDestroy_Interval(RG rg)
272: {

276:   PetscFree(rg->data);
277:   PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",NULL);
278:   PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",NULL);
279:   return(0);
280: }

282: SLEPC_EXTERN PetscErrorCode RGCreate_Interval(RG rg)
283: {
284:   RG_INTERVAL    *interval;

288:   PetscNewLog(rg,&interval);
289:   interval->a = -PETSC_MAX_REAL;
290:   interval->b = PETSC_MAX_REAL;
291:   interval->c = -PETSC_MAX_REAL;
292:   interval->d = PETSC_MAX_REAL;
293:   rg->data = (void*)interval;

295:   rg->ops->istrivial      = RGIsTrivial_Interval;
296:   rg->ops->computecontour = RGComputeContour_Interval;
297:   rg->ops->computebbox    = RGComputeBoundingBox_Interval;
298:   rg->ops->checkinside    = RGCheckInside_Interval;
299:   rg->ops->setfromoptions = RGSetFromOptions_Interval;
300:   rg->ops->view           = RGView_Interval;
301:   rg->ops->destroy        = RGDestroy_Interval;
302:   PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",RGIntervalSetEndpoints_Interval);
303:   PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",RGIntervalGetEndpoints_Interval);
304:   return(0);
305: }