C ***************************************************************** C ***************************************************************** subroutine genpro(nx,ny,xstep,ystep,procit,propol,nvmipp,nvmapp, +rmin,rmax,inform) C This subroutine generates a location problem (see the Figure). C First, a regular grid with nx horizontal points and ny vertical C points in the positive orthant is considered. The points of the C grid start at the origin with an horizontal distance of xstep and C a vertical distance of ystept. This grid will be the space at C which the cities (represented by polygons) will be distributed. C Before building the cities, an area (rectangle) of preservation C where almost nothing can be done, is defined. This area of C preservation will receive, after the construction of the cities, C an hydraulic plant of energy generation (to supply the energy to C the cities). Then, in the rest of the space, the cities are built. C At each point of the grid (out of the central region) a city C (represented by a polygon) will be built with probability procit. C The definition of the polygon uses variables nvmipp, nvmapp, rmin C and rmax in a way described in the genpol (generate polygon) C subroutine. To transmit the energy from the plant to the cities, a C tower inside each city and a tower inside the central region must C be built. The objective of the problem is to determine the C location of this towers in order to minimize the sum of the C distances from each city tower to the central one. C C On Entry: C C nx integer, C number of horizontal points in the grid, C C ny integer, C number of vertical points in the grid, C C xstep double precision, C horizontal distance between points of the grid, C C ystep double precision, C vertical distance between points of the grid, C C procit double precision, C probability of defining a city at point of the grid C (0 <= procit <= 1), C C propol double precision, C probability of the city been defined by a polygon (0 <= C propol <= 1), if the city is not defined by a polygon then C it is defined by a circle, C C nvmipp integer, C parameter for the polygon generation (described in genpol C subroutine), C C nvmapp integer, C parameter for the polygon generation (described in genpol C subroutine), C C rmin double precision, C parameter for the polygon generation (described in genpol C subroutine), C C rmax double precision, C parameter for the polygon generation (described in genpol C subroutine). C C On output: C C As described in the genpol subroutine, the output is saved in the C polyg common block. C PARAMETERS integer ncmax parameter ( ncmax = 630000 ) integer npmax parameter ( npmax = 1000000 ) integer nvsmax parameter ( nvsmax = 13 * npmax ) C SCALAR ARGUMENTS integer inform,nvmapp,nvmipp,nx,ny double precision procit,propol,rmax,rmin,xstep,ystep C COMMON SCALARS integer nc,np,totnvs double precision xdisp,xscal,ydisp,yscal C COMMON ARRAYS integer nvs(npmax) double precision ccent(2*ncmax),edges(3*nvsmax),pcent(2*npmax), + pol(nvsmax),radii(ncmax),vert(2*nvsmax) C LOCAL SCALARS integer i,j,k logical inscre,outcre,insrec double precision c,cx,cy,seed,vx,vy,lx,ly,ux,uy C EXTERNAL FUNCTIONS double precision drand external drand C EXTERNAL SUBROUTINES external genpol C COMMON BLOCKS common /circl/ ccent,radii,nc common /poly1/ nvs,np,totnvs common /poly2/ vert common /poly3/ edges common /poly4/ pcent common /poly5/ pol common /ellip/ xdisp,xscal,ydisp,yscal inform = 0 C SEED FOR THE RANDOM GENERATION seed = 760013.0d0 C DEFINE CENTRAL REGION xdisp = 0.20d0 * ( nx - 1 ) * xstep ydisp = 0.50d0 * ( ny - 1 ) * ystep xscal = 0.15d0 * ( nx - 1 ) * xstep yscal = 0.50d0 * ( ny - 1 ) * ystep C DEFINE CENTRAL POINT POLYGON nvs(1) = 4 lx = xdisp + 0.50d0 * xscal ux = xdisp + 1.50d0 * xscal ly = ydisp - 0.75d0 * yscal uy = ydisp + 0.75d0 * yscal vert(1) = lx vert(2) = ly vert(3) = lx vert(4) = uy vert(5) = ux vert(6) = uy vert(7) = ux vert(8) = ly edges(1) = - 1.0d0 edges(2) = 0.0d0 edges(3) = lx edges(4) = 0.0d0 edges(5) = 1.0d0 edges(6) = - uy edges(7) = 1.0d0 edges(8) = 0.0d0 edges(9) = - ux edges(10) = 0.0d0 edges(11) = - 1.0d0 edges(12) = ly np = 1 totnvs = 4 pcent(1) = ( lx + ux ) / 2.0d0 pcent(2) = ( ly + uy ) / 2.0d0 do i = 1,4 pol(i) = 1 end do C THESE THREE LINES BELOW SHOULD BE USED INSTEAD OF THE 27 LINES C ABOVE IF THE CENTRAL RECTANGLE WOULD NOT TO BE CONSIDERED. C nvs(1) = 0 C np = 1 C totnvs = 0 nc = 0 C DEFINE CITY-POLYGONS AND CITY-CIRCLES CENTERED AT THE GRID C POINTS. THE CITIES MUST BE OUTSIDE THE CENTRAL REGION AND NOT C COMPLETELY INSIDE THE CASSINIAN CURVE do i = 0,nx - 1 do j = 0,ny - 1 cx = i * xstep cy = j * ystep if ( drand(seed) .le. propol ) then ! GENERATE A CITY-POLYGON call genpol(cx,cy,nvmipp,nvmapp,rmin,rmax,seed,inform) if ( inform .ne. 0 ) then return end if ! TEST WHETHER THE POLYGON IS IN THE INTERIOR, IT IS ! OUTSIDE OR IT CUTS CENTRAL REGION ! 1) CITIES IN THE INTERIOR ARE FORBIDDEN TO AVOID ! INFEASIBLE PROBLEMS ! 2) CITIES THAT TOUCH THE CENTRAL REGION ARE DESIRED ! WITH THE AIM OF MAKING THE CENTRAL-REGION CONSTRAINT ! ACTIVE AT THE SOLUTION AS MANY TIMES AS POSSIBLE ! 3) CITIES THAT DO NOT TOUCH THE CENTRAL REGION ARE ! ARE INTRODUCED INTO THE PROBLEM WITH PROBABILITY PROCIT ! IT IS NOT EASY TO DETERMINE THE PREVIOUS RELATIONS (INSIDE, ! OUTSIDE, HALF AND A HALF) BETWEEN THE POLYGON AND THE ! NON-C0NVEX CENTRAL REGION. SO WE WILL CONSIDER THAT: ! A POLYGON IS INSIDE THE CENTRAL REGION IF ALL ITS VERTICES ! ARE INSIDE OR IN THE BORDER OF THE CENTRAL REGION. ! A POLYGON IS OUTSIDE THE CENTRAL REGION IF ALL ITS VERTICES ! ARE OUTSIDE OR IN THE BORDER OF THE CENTRAL REGION. ! A POLYGON CUTS THE CENTRAL REGION IF IT HAS AT LEAST A ! VERTEX INSIDE OR IN THE BORDER OF THE CENTRAL REGION AND ! AT LEAST A VERTEX OUTSIDE OR IN THE BORDER OF THE CENTRAL ! REGION. SO, TO BE HALF AND HALF IS EQUIVALENT TO NOT TO ! INSIDE, NOR OUTISDE. inscre = .true. outcre = .true. do k = 1,nvs(np+1) vx = ( vert(2*(totnvs+k)-1) - xdisp ) / xscal vy = ( vert(2*(totnvs+k)) - ydisp ) / yscal c = vx ** 2 + vy ** 2 - 1.0d0 if ( c .gt. 0.0d0 ) then inscre = .false. end if if ( c .lt. 0.0d0 ) then outcre = .false. end if end do if ( cx .ge. lx - rmax .and. cx .le. ux + rmax .and. + cy .ge. ly - rmax .and. cy .le. uy + rmax ) then insrec = .true. else insrec = .false. end if ! ADD THE CITY-POLYGON TO THE PROBLEM if ( .not. insrec .and. + ( ( .not. inscre .and. .not. outcre ) .or. + ( outcre .and. drand(seed) .le. procit ) ) ) then np = np + 1 totnvs = totnvs + nvs(np) end if else ! GENERATE A CITY-CIRCLE call gencir(cx,cy,rmin,rmax,seed,inform) if ( inform .ne. 0 ) then return end if ! VERIFY WHERE IT IS inscre = .true. outcre = .true. do k = 1,4 if ( k .eq. 1 ) then vx = ( cx - radii(nc + 1) - xdisp ) / xscal vy = ( cy - ydisp ) / yscal else if ( k .eq. 2 ) then vx = ( cx + radii(nc + 1) - xdisp ) / xscal vy = ( cy - ydisp ) / yscal else if ( k .eq. 3 ) then vx = ( cx - xdisp ) / xscal vy = ( cy - radii(nc + 1) - ydisp ) / yscal else if ( k .eq. 4 ) then vx = ( cx - xdisp ) / xscal vy = ( cy + radii(nc + 1) - ydisp ) / yscal end if c = vx ** 2 + vy ** 2 - 1.0d0 if ( c .gt. 0.0d0 ) then inscre = .false. end if if ( c .lt. 0.0d0 ) then outcre = .false. end if end do if ( cx .ge. lx - rmax .and. cx .le. ux + rmax .and. + cy .ge. ly - rmax .and. cy .le. uy + rmax ) then insrec = .true. else insrec = .false. end if ! ADD THE CITY-CIRCLE TO THE PROBLEM if ( .not. insrec .and. + ( ( .not. inscre .and. .not. outcre ) .or. + ( outcre .and. drand(seed) .le. procit ) ) ) then nc = nc + 1 end if end if end do end do end C ***************************************************************** C ***************************************************************** subroutine gencir(cx,cy,rmin,rmax,seed,inform) C PARAMETERS integer ncmax parameter ( ncmax = 630000 ) C SCALAR ARGUMENTS integer inform double precision cx,cy,rmax,rmin,seed C COMMON SCALARS integer nc C COMMON ARRAYS double precision ccent(2*ncmax),radii(ncmax) C EXTERNAL FUNCTIONS double precision drand external drand C COMMON BLOCKS common /circl/ ccent,radii,nc inform = 0 C VERIFY SPACE AVAILABILITY FOR CIRCLES if ( nc .eq. ncmax ) then write (*,fmt=*) 'THE MAXIMUM NUMBER OF CIRCLES WAS ACHIEVED.' write (*,fmt=*) 'INCREASE NCMAX OR REDUCE THE NUMBER OF GRID ' write (*,fmt=*) 'POINTS (NX*NY) OR THE PROBABILITY OF HAVING ' write (*,fmt=*) 'A CITY-CIRCLE AT A GRID POINT' write (*,fmt=*) '(PROCIT*(1-PROPOL)).' inform = - 1 return end if C SAVE THE CENTER ccent(2 * nc + 1) = cx ccent(2 * nc + 2) = cy C GENERATE THE RADIUS radii(nc + 1) = rmin + ( rmax - rmin ) * drand(seed) end C ***************************************************************** C ***************************************************************** subroutine genpol(cx,cy,nvmipp,nvmapp,rmin,rmax,seed,inform) C This subroutine generates a polygon in R^2 with its C vertices in a sphere centered at point (cx,cy). The C number of vertices is randomly generated C satisfying nvmipp <= number of vertices <= C nvmapp. The ratio of the sphere is also randomly C generated satisfying rmin <= ratio < rmax. The C generated polygon is stored in the common block "polyg". C C On Entry: C C cx double precision, C first coordinate of the center of the sphere, C C cy double precision, C second coordinate of the center of the sphere, C C nvmipp integer, C minimum number of vertices, C C nvmapp integer, C maximum number of vertices, C C rmin double precision, C minimum ratio of the sphere, C C rmax double precision, C maximum ratio of the sphere, C C seed double precision, C seed for the random generation. C C On Output: C C inform integer C termination parameter C C 0 = everything was ok C C < 0 = there was not enough memory space. C C The generated polygon is stored in the polyg common block C described below. C C Common block polyg: C C common /polyg/nvs,vert,edges,np,totnvs C C This structure represents, at any time, np polygons. C Position i of array nvs indicates the number of vertices C of polygon i. Arrays vert and edges store the C vertices and edges of the polygons. C C For example, if nvs(1) = 3 it indicates that the first C polygon has 3 vertices (edges). Then, if the vertices C are (x1,y1), (x2,y2) and (x3,y3), we have that vert(1) C = x1, vert(2) = y1, vert(3) = x2, vert(4) = y2, vert(5) C = x3, and vert(6) = y3. And, if the edges (written as C ax + by + c = 0) are a1 x + b1 y + c1 = 0, a2 x + b2 y C + c2 = 0, and a3 x + b3 y + c3 = 0 then edges(1) = a1, C edges(2) = b1, edges(3) = c1, edges(4) = a2, edges(5) = C b2, edges(6) = c2, edges(7) = a3, edges(8) = b3 and C edges(9) = c3. C C totnvs indicates the total number of vertices C of the set of polygons. This information is used when C a new polygon is created to know the first free position C of arrays vert and edges at which the vertices and edges C of the new polygon will be saved. C C Two additional details: C C 1) For each polygon, the vertices are ordered clockwise C and edge i corresponds to the edge between vertices C i and i+1 (0 if i=n). C C 2) For each edge of the form ax + bx + c = 0, constants C a, b and c are chosen in such a way that C (|a| = 1 or |b| = 1) and (a cx + b cy + c <= 0). C PARAMETERS integer npmax parameter ( npmax = 1000000 ) integer nvsmax parameter ( nvsmax = 13 * npmax ) double precision pi parameter ( pi = 3.1415926535898d0 ) C SCALAR ARGUMENTS integer inform,nvmapp,nvmipp double precision cx,cy,rmax,rmin,seed C COMMON SCALARS integer np,totnvs C COMMON ARRAYS integer nvs(npmax) double precision edges(3*nvsmax),pcent(2*npmax),pol(nvsmax), + vert(2*nvsmax) C LOCAL SCALARS integer i double precision dist,lseed,mindist,r C LOCAL ARRAYS double precision angl(nvsmax) C EXTERNAL FUNCTIONS double precision drand external drand C EXTERNAL SUBROUTINES external class,constr C INTRINSIC FUNCTIONS intrinsic cos,int,sin C COMMON BLOCKS common /poly1/ nvs,np,totnvs common /poly2/ vert common /poly3/ edges common /poly4/ pcent common /poly5/ pol inform = 0 C VERIFY SPACE AVAILABILITY FOR POLYGONS if ( np .eq. npmax ) then write (*,fmt=*) 'THE MAXIMUM NUMBER OF POLYGONS WAS ACHIEVED.' write (*,fmt=*) 'INCREASE NPMAX OR REDUCE THE NUMBER OF GRID ' write (*,fmt=*) 'POINTS (NX*NY) OR THE PROBABILITY OF HAVING ' write (*,fmt=*) 'A CITY-POLYGON AT A GRID POINT' write (*,fmt=*) '(PROCIT*PROPOL).' inform = - 1 return end if C GENERATE THE NUMBER OF VERTICES lseed = 157318.0d0 + seed nvs(np+1) = nvmipp + int((nvmapp-nvmipp+1)*drand(lseed)) C VERIFY SPACE AVAILABILITY FOR VERTICES AND SIDES if ( totnvs + nvs(np+1) .gt. nvsmax ) then write (*,fmt=*) 'THE MAXIMUM NUMBER OF POLYGONS VERTICES WAS' write (*,fmt=*) 'ACHIEVED. INCREASE NVSMAX OR REDUCE THE' write (*,fmt=*) 'AVERAGE NUMBER OF VERTICES (NVMAPP-NVMIPP)/2' write (*,fmt=*) 'OR THE NUMBER OF POLYGONS. TO REDUCE THE' write (*,fmt=*) 'NUMBER OF POLYGONS, REDUCE THE NUMBER OF' write (*,fmt=*) 'GRID POINTS (NX*NY) OR THE PROBABILITY OF' write (*,fmt=*) 'HAVING A CITY-POLYGONS AT A GRID POINT' write (*,fmt=*) '(PROCIT*PROPOL).' inform = - 1 return end if C SAVE THE "CENTER" TO BE USED AS INITIAL POINT pcent(2*np+1) = cx pcent(2*np+2) = cy C IDENTIFY THE CONSTRAINTS WITH THE POLYGON do i = totnvs + 1,totnvs + nvs(np+1) pol(i) = np + 1 end do C GENERATE THE RADIUS OF THE SPHERE r = rmin + ( rmax - rmin ) * drand( seed ) C GENERATE ALL ANGLES SATISFYING 0 <= ANGLE_I < 2*PI 10 continue do i = 1,nvs(np + 1) angl(i) = 2 * pi * drand( lseed ) end do C CLASSIFY THE ANGLES IN DECREASING ORDER call class(nvs(np + 1),angl) C CONSTRUCT THE VERTICES do i = 1,nvs(np + 1) vert(2 * ( totnvs + i ) - 1) = cx + r * cos( angl(i) ) vert(2 * ( totnvs + i ) ) = cy + r * sin( angl(i) ) end do C FOR NUMERICAL STABILITY IN THE CONSTRAINT CALCULATION, C AVOID TOO SIMILAR ANGLES (EQUIVALENT TO TOO NEAR POINTS) mindist = 1.0D+99 do i = totnvs + 1,totnvs + nvs(np+1) - 1 dist = ( vert(2*i-1) - vert(2*i+1) ) ** 2 + + ( vert(2*i) - vert(2*i+2) ) ** 2 if ( dist .lt. mindist ) then mindist = dist end if end do i = totnvs + nvs(np+1) dist = ( vert(2*i-1) - vert(2*totnvs+1) ) ** 2 + + ( vert(2*i) - vert(2*totnvs+2) ) ** 2 if ( dist .lt. mindist ) then mindist = dist end if if ( mindist .lt. 1.0D-8 ) then C write (*,fmt=*) 'DISCARDING ANGLES.' go to 10 end if C CONSTRUCT THE EDGES do i = totnvs + 1,totnvs + nvs(np+1) - 2 call constr(vert(2*i-1),vert(2*i),vert(2*i+1),vert(2*i+2), + vert(2*i+3),vert(2*i+4),edges(3*i-2),edges(3*i-1), + edges(3*i),inform) if ( inform .ne. 0 ) then return end if end do i = totnvs + nvs(np+1) - 1 call constr(vert(2*i-1),vert(2*i),vert(2*i+1),vert(2*i+2), + vert(2*totnvs+1),vert(2*totnvs+2),edges(3*i-2), + edges(3*i-1),edges(3*i),inform) if ( inform .ne. 0 ) then return end if i = totnvs + nvs(np+1) call constr(vert(2*i-1),vert(2*i),vert(2*totnvs+1), + vert(2*totnvs+2),vert(2*totnvs+3),vert(2*totnvs+4), + edges(3*i-2),edges(3*i-1),edges(3*i),inform) if ( inform .ne. 0 ) then return end if end C ***************************************************************** C ***************************************************************** subroutine constr(x1,y1,x2,y2,x3,y3,a,b,c,inform) C This subroutine computes the real constants a, b and c of C the straight line ax + by + c = 0 in R^2 defined by the C points (x1,y1) and (x2,y2); such that the point (x3,y3) C satisfies the constraint ax + by + c <= 0. C C On Entry: C C x1 double precision, C first coordinate of point (x1,y1), C C y1 double precision, C second coordinate of point (x1,y1), C C x2 double precision, C first coordinate of point (x2,y2), C C y2 double precision, C second coordinate of point (x2,y2), C C x3 double precision, C first coordinate of point (x3,y3), C C y3 double precision, C second coordinate of point (x3,y3). C C On Return C C a,b,c double precision C the desired constants. C SCALAR ARGUMENTS integer inform double precision a,b,c,x1,x2,x3,y1,y2,y3 if (x1.eq.x2 .and. y1.eq.y2) then write (*,fmt=*) + 'ERROR IN FUNCTION CONSTRAINT: X1=X2 AND Y1=Y2' inform = - 1 return end if if (y1.ne.y2) then a = 1.0d0 b = - (x2-x1)/ (y2-y1) c = - (x1+b*y1) else a = 0.0d0 b = 1.0d0 c = -y1 end if if (a*x3+b*y3+c.gt.0.0d0) then a = -a b = -b c = -c end if end C ***************************************************************** C ***************************************************************** subroutine class(n,x) C This subroutine classifies the elements of a vector in C decreasing order, i.e., on output: x(1) >= x(2) >= C ... >= x(n). C C On Entry: C C n integer, C number of elements of the vector to be classified, C C x double precision x(n), C vector to be classified. C C On Return C C x double precision x(n), C classified vector. C SCALAR ARGUMENTS integer n C ARRAY ARGUMENTS double precision x(n) C LOCAL SCALARS double precision aux,xmax integer i,j,pos do i = 1,n xmax = x(i) pos = i do j = i + 1,n if (x(j).gt.xmax) then xmax = x(j) pos = j end if end do if (pos.ne.i) then aux = x(i) x(i) = x(pos) x(pos) = aux end if end do end