Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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             : */
      14             : 
      15             : #include <slepc/private/rgimpl.h>      /*I "slepcrg.h" I*/
      16             : #include <petscdraw.h>
      17             : 
      18             : typedef struct {
      19             :   PetscReal   a,b;     /* interval in the real axis */
      20             :   PetscReal   c,d;     /* interval in the imaginary axis */
      21             : } RG_INTERVAL;
      22             : 
      23         100 : static PetscErrorCode RGIntervalSetEndpoints_Interval(RG rg,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
      24             : {
      25         100 :   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
      26             : 
      27         100 :   PetscFunctionBegin;
      28         100 :   PetscCheck(a || b || c || d,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"At least one argument must be nonzero");
      29         100 :   PetscCheck(a!=b || !a,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, endpoints must be distinct (or both zero)");
      30         100 :   PetscCheck(a<=b,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be a<b");
      31         100 :   PetscCheck(c!=d || !c,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, endpoints must be distinct (or both zero)");
      32         100 :   PetscCheck(c<=d,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"Badly defined interval, must be c<d");
      33             : #if !defined(PETSC_USE_COMPLEX)
      34         100 :   PetscCheck(c==-d,PetscObjectComm((PetscObject)rg),PETSC_ERR_ARG_WRONG,"In real scalars the region must be symmetric wrt real axis");
      35             : #endif
      36         100 :   ctx->a = a;
      37         100 :   ctx->b = b;
      38         100 :   ctx->c = c;
      39         100 :   ctx->d = d;
      40         100 :   PetscFunctionReturn(PETSC_SUCCESS);
      41             : }
      42             : 
      43             : /*@
      44             :    RGIntervalSetEndpoints - Sets the parameters defining the interval region.
      45             : 
      46             :    Logically Collective
      47             : 
      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
      54             : 
      55             :    Options Database Keys:
      56             : .  -rg_interval_endpoints - the four endpoints
      57             : 
      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.
      62             : 
      63             :    When PETSc is built with real scalars, the region must be symmetric with
      64             :    respect to the real axis.
      65             : 
      66             :    Level: advanced
      67             : 
      68             : .seealso: RGIntervalGetEndpoints()
      69             : @*/
      70         100 : PetscErrorCode RGIntervalSetEndpoints(RG rg,PetscReal a,PetscReal b,PetscReal c,PetscReal d)
      71             : {
      72         100 :   PetscFunctionBegin;
      73         100 :   PetscValidHeaderSpecific(rg,RG_CLASSID,1);
      74         400 :   PetscValidLogicalCollectiveReal(rg,a,2);
      75         400 :   PetscValidLogicalCollectiveReal(rg,b,3);
      76         400 :   PetscValidLogicalCollectiveReal(rg,c,4);
      77         400 :   PetscValidLogicalCollectiveReal(rg,d,5);
      78         100 :   PetscTryMethod(rg,"RGIntervalSetEndpoints_C",(RG,PetscReal,PetscReal,PetscReal,PetscReal),(rg,a,b,c,d));
      79         100 :   PetscFunctionReturn(PETSC_SUCCESS);
      80             : }
      81             : 
      82         346 : static PetscErrorCode RGIntervalGetEndpoints_Interval(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
      83             : {
      84         346 :   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
      85             : 
      86         346 :   PetscFunctionBegin;
      87         346 :   if (a) *a = ctx->a;
      88         346 :   if (b) *b = ctx->b;
      89         346 :   if (c) *c = ctx->c;
      90         346 :   if (d) *d = ctx->d;
      91         346 :   PetscFunctionReturn(PETSC_SUCCESS);
      92             : }
      93             : 
      94             : /*@C
      95             :    RGIntervalGetEndpoints - Gets the parameters that define the interval region.
      96             : 
      97             :    Not Collective
      98             : 
      99             :    Input Parameter:
     100             : .  rg - the region context
     101             : 
     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
     107             : 
     108             :    Level: advanced
     109             : 
     110             : .seealso: RGIntervalSetEndpoints()
     111             : @*/
     112         346 : PetscErrorCode RGIntervalGetEndpoints(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
     113             : {
     114         346 :   PetscFunctionBegin;
     115         346 :   PetscValidHeaderSpecific(rg,RG_CLASSID,1);
     116         346 :   PetscUseMethod(rg,"RGIntervalGetEndpoints_C",(RG,PetscReal*,PetscReal*,PetscReal*,PetscReal*),(rg,a,b,c,d));
     117         346 :   PetscFunctionReturn(PETSC_SUCCESS);
     118             : }
     119             : 
     120           5 : static PetscErrorCode RGView_Interval(RG rg,PetscViewer viewer)
     121             : {
     122           5 :   RG_INTERVAL    *ctx = (RG_INTERVAL*)rg->data;
     123           5 :   PetscBool      isdraw,isascii;
     124           5 :   int            winw,winh;
     125           5 :   PetscDraw      draw;
     126           5 :   PetscDrawAxis  axis;
     127           5 :   PetscReal      a,b,c,d,ab,cd,lx,ly,w,scale=1.2;
     128             : 
     129           5 :   PetscFunctionBegin;
     130           5 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw));
     131           5 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     132           5 :   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           1 :   else if (isdraw) {
     134           1 :     PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw));
     135           1 :     PetscCall(PetscDrawCheckResizedWindow(draw));
     136           1 :     PetscCall(PetscDrawGetWindowSize(draw,&winw,&winh));
     137           1 :     winw = PetscMax(winw,1); winh = PetscMax(winh,1);
     138           1 :     PetscCall(PetscDrawClear(draw));
     139           1 :     PetscCall(PetscDrawSetTitle(draw,"Interval region"));
     140           1 :     PetscCall(PetscDrawAxisCreate(draw,&axis));
     141           1 :     a  = ctx->a*rg->sfactor;
     142           1 :     b  = ctx->b*rg->sfactor;
     143           1 :     c  = ctx->c*rg->sfactor;
     144           1 :     d  = ctx->d*rg->sfactor;
     145           1 :     lx = b-a;
     146           1 :     ly = d-c;
     147           1 :     ab = (a+b)/2;
     148           1 :     cd = (c+d)/2;
     149           1 :     w  = scale*PetscMax(lx/winw,ly/winh)/2;
     150           1 :     PetscCall(PetscDrawAxisSetLimits(axis,ab-w*winw,ab+w*winw,cd-w*winh,cd+w*winh));
     151           1 :     PetscCall(PetscDrawAxisDraw(axis));
     152           1 :     PetscCall(PetscDrawAxisDestroy(&axis));
     153           1 :     PetscCall(PetscDrawRectangle(draw,a,c,b,d,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE,PETSC_DRAW_BLUE));
     154           1 :     PetscCall(PetscDrawFlush(draw));
     155           1 :     PetscCall(PetscDrawSave(draw));
     156           1 :     PetscCall(PetscDrawPause(draw));
     157             :   }
     158           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     159             : }
     160             : 
     161        9717 : static PetscErrorCode RGIsTrivial_Interval(RG rg,PetscBool *trivial)
     162             : {
     163        9717 :   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
     164             : 
     165        9717 :   PetscFunctionBegin;
     166        9744 :   if (rg->complement) *trivial = (ctx->a==ctx->b && ctx->c==ctx->d)? PETSC_TRUE: PETSC_FALSE;
     167       10167 :   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        9717 :   PetscFunctionReturn(PETSC_SUCCESS);
     169             : }
     170             : 
     171          38 : static PetscErrorCode RGComputeContour_Interval(RG rg,PetscInt n,PetscScalar *cr,PetscScalar *ci)
     172             : {
     173          38 :   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
     174          38 :   PetscInt    i,N,Nv,Nh,k1,k0;
     175          38 :   PetscReal   hv,hh,t;
     176             : 
     177          38 :   PetscFunctionBegin;
     178          38 :   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          38 :   if (ctx->a==ctx->b || ctx->c==ctx->d) {
     180          34 :     if (ctx->a==ctx->b) {hv = (ctx->d-ctx->c)/(n-1); hh = 0.0;}
     181          34 :     else {hh = (ctx->b-ctx->a)/(n-1); hv = 0.0;}
     182      131494 :     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      131460 :       if (cr) cr[i] = ctx->a+hh*i;
     187      131460 :       if (ci) ci[i] = ctx->c+hv*i;
     188             : #endif
     189             :     }
     190             :   } else {
     191           4 :     PetscCheck(n>3,PetscObjectComm((PetscObject)rg),PETSC_ERR_SUP,"Minimum number of contour points: 4");
     192           4 :     N = n/2;
     193           4 :     t = ((ctx->d-ctx->c)/(ctx->d-ctx->c+ctx->b-ctx->a))*N;
     194           4 :     Nv = t-PetscFloorReal(t)>0.5?PetscCeilReal(t):PetscFloorReal(t);
     195           4 :     if (Nv==0) Nv++;
     196           2 :     else if (Nv==N) Nv--;
     197           4 :     Nh = N-Nv;
     198           4 :     hh = (ctx->b-ctx->a)/Nh;
     199           4 :     hv = (ctx->d-ctx->c)/Nv;
     200             :     /* positive imaginary part first */
     201           4 :     k1 = Nv/2+1;
     202           4 :     k0 = Nv-k1;
     203             : 
     204           4 :     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           0 :       if (cr) {cr[i-k1] = ctx->b;      cr[i-k1+N] = ctx->a;}
     210           0 :       if (ci) {ci[i-k1] = ctx->c+i*hv; ci[i-k1+N] = ctx->d-i*hv;}
     211             : #endif
     212             :     }
     213       10010 :     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       10006 :       if (cr) {cr[i+k0] = ctx->b-i*hh; cr[i+k0+N] = ctx->a+i*hh;}
     219       10006 :       if (ci) {ci[i+k0] = ctx->d;      ci[i+k0+N] = ctx->c;}
     220             : #endif
     221             :     }
     222           8 :     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           4 :       if (cr) {cr[i+k0+Nh] = ctx->a; cr[i+k0+Nh+N] = ctx->b;}
     228           4 :       if (ci) {ci[i+k0+Nh] = ctx->d+i*hv; ci[i+k0+Nh+N] = ctx->c-i*hv;}
     229             : #endif
     230             :     }
     231             :   }
     232          38 :   PetscFunctionReturn(PETSC_SUCCESS);
     233             : }
     234             : 
     235           4 : static PetscErrorCode RGComputeBoundingBox_Interval(RG rg,PetscReal *a,PetscReal *b,PetscReal *c,PetscReal *d)
     236             : {
     237           4 :   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
     238             : 
     239           4 :   PetscFunctionBegin;
     240           4 :   if (a) *a = ctx->a;
     241           4 :   if (b) *b = ctx->b;
     242           4 :   if (c) *c = ctx->c;
     243           4 :   if (d) *d = ctx->d;
     244           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     245             : }
     246             : 
     247           9 : static PetscErrorCode RGComputeQuadrature_Interval(RG rg,RGQuadRule quad,PetscInt n,PetscScalar *z,PetscScalar *zn,PetscScalar *w)
     248             : {
     249           9 :   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
     250           9 :   PetscReal   theta,max_w=0.0,radius=1.0;
     251           9 :   PetscScalar tmp,tmp2,center=0.0;
     252           9 :   PetscInt    i,j;
     253             : 
     254           9 :   PetscFunctionBegin;
     255           9 :   if (quad == RG_QUADRULE_CHEBYSHEV) {
     256           8 :     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         244 :     for (i=0;i<n;i++) {
     258         236 :       theta = PETSC_PI*(i+0.5)/n;
     259         236 :       zn[i] = PetscCosReal(theta);
     260         236 :       w[i]  = PetscCosReal((n-1)*theta)/n;
     261         236 :       if (ctx->c==ctx->d) z[i] = ((ctx->b-ctx->a)*(zn[i]+1.0)/2.0+ctx->a)*rg->sfactor;
     262           0 :       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         236 :         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           1 :     center = rg->sfactor*(ctx->b+ctx->a)/2.0;
     275             : #endif
     276           2 :     radius = PetscSqrtReal(PetscPowRealInt(rg->sfactor*(ctx->b-ctx->a)/2.0,2)+PetscPowRealInt(rg->sfactor*(ctx->d-ctx->c)/2.0,2));
     277          25 :     for (i=0;i<n;i++) {
     278          24 :       zn[i] = (z[i]-center)/radius;
     279          24 :       tmp = 1.0; tmp2 = 1.0;
     280         600 :       for (j=0;j<n;j++) {
     281         576 :         tmp *= z[j];
     282         576 :         if (i != j) tmp2 *= z[j]-z[i];
     283             :       }
     284          24 :       w[i] = tmp/tmp2;
     285          37 :       max_w = PetscMax(PetscAbsScalar(w[i]),max_w);
     286             :     }
     287          25 :     for (i=0;i<n;i++) w[i] /= (PetscScalar)max_w;
     288             :   }
     289           9 :   PetscFunctionReturn(PETSC_SUCCESS);
     290             : }
     291             : 
     292      130793 : static PetscErrorCode RGCheckInside_Interval(RG rg,PetscReal dx,PetscReal dy,PetscInt *inside)
     293             : {
     294      130793 :   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
     295             : 
     296      130793 :   PetscFunctionBegin;
     297      130793 :   if (dx>ctx->a && dx<ctx->b) *inside = 1;
     298       74415 :   else if (dx==ctx->a || dx==ctx->b) *inside = 0;
     299       74400 :   else *inside = -1;
     300      130793 :   if (*inside>=0) {
     301       56393 :     if (dy>ctx->c && dy<ctx->d) ;
     302       17058 :     else if (dy==ctx->c || dy==ctx->d) *inside = 0;
     303        9785 :     else *inside = -1;
     304             :   }
     305      130793 :   PetscFunctionReturn(PETSC_SUCCESS);
     306             : }
     307             : 
     308           5 : static PetscErrorCode RGIsAxisymmetric_Interval(RG rg,PetscBool vertical,PetscBool *symm)
     309             : {
     310           5 :   RG_INTERVAL *ctx = (RG_INTERVAL*)rg->data;
     311             : 
     312           5 :   PetscFunctionBegin;
     313           5 :   if (vertical) *symm = (ctx->a == -ctx->b)? PETSC_TRUE: PETSC_FALSE;
     314           2 :   else *symm = (ctx->c == -ctx->d)? PETSC_TRUE: PETSC_FALSE;
     315           5 :   PetscFunctionReturn(PETSC_SUCCESS);
     316             : }
     317             : 
     318         835 : static PetscErrorCode RGSetFromOptions_Interval(RG rg,PetscOptionItems *PetscOptionsObject)
     319             : {
     320         835 :   PetscBool      flg;
     321         835 :   PetscInt       k;
     322         835 :   PetscReal      array[4]={0,0,0,0};
     323             : 
     324         835 :   PetscFunctionBegin;
     325         835 :   PetscOptionsHeadBegin(PetscOptionsObject,"RG Interval Options");
     326             : 
     327         835 :     k = 4;
     328         835 :     PetscCall(PetscOptionsRealArray("-rg_interval_endpoints","Interval endpoints (two or four real values separated with a comma without spaces)","RGIntervalSetEndpoints",array,&k,&flg));
     329         835 :     if (flg) {
     330          49 :       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          49 :       PetscCall(RGIntervalSetEndpoints(rg,array[0],array[1],array[2],array[3]));
     332             :     }
     333             : 
     334         835 :   PetscOptionsHeadEnd();
     335         835 :   PetscFunctionReturn(PETSC_SUCCESS);
     336             : }
     337             : 
     338         876 : static PetscErrorCode RGDestroy_Interval(RG rg)
     339             : {
     340         876 :   PetscFunctionBegin;
     341         876 :   PetscCall(PetscFree(rg->data));
     342         876 :   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",NULL));
     343         876 :   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",NULL));
     344         876 :   PetscFunctionReturn(PETSC_SUCCESS);
     345             : }
     346             : 
     347         876 : SLEPC_EXTERN PetscErrorCode RGCreate_Interval(RG rg)
     348             : {
     349         876 :   RG_INTERVAL    *interval;
     350             : 
     351         876 :   PetscFunctionBegin;
     352         876 :   PetscCall(PetscNew(&interval));
     353         876 :   interval->a = -PETSC_MAX_REAL;
     354         876 :   interval->b = PETSC_MAX_REAL;
     355         876 :   interval->c = -PETSC_MAX_REAL;
     356         876 :   interval->d = PETSC_MAX_REAL;
     357         876 :   rg->data = (void*)interval;
     358             : 
     359         876 :   rg->ops->istrivial         = RGIsTrivial_Interval;
     360         876 :   rg->ops->computecontour    = RGComputeContour_Interval;
     361         876 :   rg->ops->computebbox       = RGComputeBoundingBox_Interval;
     362         876 :   rg->ops->computequadrature = RGComputeQuadrature_Interval;
     363         876 :   rg->ops->checkinside       = RGCheckInside_Interval;
     364         876 :   rg->ops->isaxisymmetric    = RGIsAxisymmetric_Interval;
     365         876 :   rg->ops->setfromoptions    = RGSetFromOptions_Interval;
     366         876 :   rg->ops->view              = RGView_Interval;
     367         876 :   rg->ops->destroy           = RGDestroy_Interval;
     368         876 :   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGIntervalSetEndpoints_C",RGIntervalSetEndpoints_Interval));
     369         876 :   PetscCall(PetscObjectComposeFunction((PetscObject)rg,"RGIntervalGetEndpoints_C",RGIntervalGetEndpoints_Interval));
     370         876 :   PetscFunctionReturn(PETSC_SUCCESS);
     371             : }