C ================================================================= C File: location.f C ================================================================= C ================================================================= C Module: Subroutines that define the problem C ================================================================= C Last update of any of the component of this module: C September 4, 2007. C Users are encouraged to download periodically updated versions of C this code at the COLLECTION home page: C C www.ime.usp.br/~egbirgin/collection/ C C and periodically updated versions of the TANGO Project solvers at C the TANGO home page: C C www.ime.usp.br/~egbirgin/tango/ C ================================================================= C Location problem C ---------------- C This is a variant of the family of location problems introduced C in [1]. In the original problem, given a set of np disjoint C polygons P1, P2, ..., Pnp in R^2 one wishes to find the point z1 C in P1 that minimizes the sum of the distances to the other C polygons. In this variant [2], we have, in addition to the np C polygons, nc circles. Moreover, there is an ellipse which has a C non empty intersection with P1 and such that z1 must be inside C the ellipse and zi, i = 2, ..., np+nc must be outside. The C detailed formulation can be found in [2]. C C [1] E. G. Birgin, J. M. Martínez and M. Raydan, "Algorithm 813: C SPG - software for convex-constrained optimization", ACM C Transactions on Mathematical Software 27, pp. 340-349, 2001. C C [2] R. Andreani, E. G. Birgin, J. M. Martínez and M. L. Schuverdt, C "On Augmented Lagrangian methods with general lower-level C constraints", submitted, 2005. C ****************************************************************** C ****************************************************************** subroutine inip(n,x,l,u,m,lambda,equatn,linear) implicit none C SCALAR ARGUMENTS integer m,n C ARRAY ARGUMENTS logical equatn(*),linear(*) double precision l(*),lambda(*),u(*),x(*) C This subroutine must set some problem data. For achieving this C objective YOU MUST MODIFY it according to your problem. See below C where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C This subroutine has no input parameters. C C On Return C C n integer, C number of variables, C C x double precision x(n), C initial point, C C l double precision l(n), C lower bounds on x, C C u double precision u(n), C upper bounds on x, C C m integer, C number of constraints (excluding the bounds), C C lambda double precision lambda(m), C initial estimation of the Lagrange multipliers, C C equatn logical equatn(m) C for each constraint j, set equatn(j) = .true. if it is an C equality constraint of the form c_j(x) = 0, and set C equatn(j) = .false. if it is an inequality constraint of C the form c_j(x) <= 0, C C linear logical linear(m) C for each constraint j, set linear(j) = .true. if it is a C linear constraint, and set linear(j) = .false. if it is a C nonlinear constraint. C PARAMETERS integer mmax,ncmax,nmax,npmax,nvsmax parameter ( ncmax = 1000 ) parameter ( npmax = 1000 ) parameter ( nvsmax = 13 * npmax ) parameter ( nmax = 2 * ( npmax + ncmax ) ) parameter ( mmax = npmax + ncmax ) C COMMON SCALARS integer nc,np,nx,ny,totnvs double precision rmax,xdisp,xscal,xstep,ydisp,yscal,ystep C COMMON ARRAYS integer nvs(npmax),pol(nvsmax) double precision ccent(2*ncmax),edges(3*nvsmax),pcent(2*npmax), + radii(ncmax),vert(2*nvsmax) C LOCAL SCALARS integer i,inform,nvmapp,nvmipp,pnum double precision procit,propol,rmin,secmar,step C EXTERNAL SUBROUTINES external genpro C COMMON BLOCKS common /const/ nx,ny,xstep,ystep,rmax 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 C DEFINE SOME PROBLEM PARAMETERS C GRID STEP (THIS IS JUST FOR SCALING THE PROBLEM) step = 1.0d0 C SECURITY MARGIN (BETWEEN 0 AND 1) TO AVOID CITIES OVERLAPPING secmar = 0.1d0 C DIFFERENT STEPS MAY BE USED IN THE GRID xstep = step ystep = step C MINIMUM AND MAXIMUM RADII FOR THE CIRCLE- AND POLYGON-CITIES C GENERATION rmin = 0.5d0 * ( step / 2.0d0 ) rmax = ( 1.0d0 - secmar ) * ( step / 2.0d0 ) C READ PROBLEM VARIABLE DATA C THIS PROBLEM DATA DEFINES DIFFERENT INSTANCIES OF THE PROBLEM AND C IS RELATED TO THE NUMBER OF VARIABLES AND CONSTRAINTS WHICH C APPROXIMATELY ARE: C C - THE NUMBER OF POLYGON-CITIES IS APPROXIMATELY: C NPC = NX * NY * PROCIT * PROPOL C C - THE NUMBER OF CIRCLE-CITIES IS APPROXIMATELY: C NCC = MX 8 NY * PROCIT * ( 1 - PROPOL ) C C - THE TOTAL NUMBER IF CITIES IS APPROXIMATELY: C NC = NPC + NCC C C - NUMBER OF VARIABLES: C 2 * ( NC + 1 ) C C - NUMBER OF LINEAR CONSTRAINTS: C NPC * (NVMAPP - NVMIPP ) / 2 + 4 C C - NUMBER OF NONLINEAR CONSTRAINTS: C NCC + NC + 1 C C WHERE: C C - PNUM IS A NUMBER THAT IDENTIFIES THE PROBLEM INSTANCE, C - NX IS THE NUMBER OF POINTS IN THE GRID ABSCISSA, C - NY IS THE NUMBER OF POINTS IN THE GRID ORDINATE, C - PROCIT IS THE PROBABILITY OF HAVING A CITY AT A GRID POINT, C - PROPOL IS THE PROBABILITY OF A CITY TO BE A POLYGON (THE CITIES C WHICH ARE NOT POLYGONS ARE CIRCLES), C - NVMIPP IS THE MINIMUM NUNBER OF VERTICES OF A POLYGON, C - NVMAPP IS THE MAXIMUM NUMBER OF VERTICES OF A POLYGON. read (*,fmt=*) pnum,nx,ny,procit,propol,nvmipp,nvmapp if ( procit .lt. 0.0d0 .or. procit .gt. 1.0d0 ) then write (*,fmt=*) 'PROCIT MUST BE BETWEEN ZERO AND ONE' stop end if if (propol .lt. 0.0d0 .or. propol .gt. 1.0d0 ) then write (*,fmt=*) 'PROPOL MUST BE BETWEEN ZERO AND ONE' stop end if if ( nvmipp .lt. 3 ) then write (*,fmt=*) 'NVMIPP MUST BE GREATER THAN OR EQUAL TO 3' stop end if if ( nvmapp .lt. nvmipp ) then write (*,fmt=*) 'NVMAPP MUST BE GREATER THAN OR EQUAL TO' write (*,fmt=*) 'NVMIPP. INCREASE NVMAPP OR REDUCE NVMIPP.' stop end if C GENERATE THE PROBLEM call genpro(nx,ny,xstep,ystep,procit,propol,nvmipp,nvmapp,rmin, +rmax,inform) if ( inform .ne. 0 ) then write (*,fmt=*) 'THERE WAS AN ERROR IN THE PROBLEM GENERATION' write (*,fmt=*) '(PROBABLY DUE TO MEMORY SPACE AVAILABILITY).' stop end if C Number of variables n = 2 * ( np + nc ) C Initial point do i = 1,np x(2*i-1) = pcent(2*i-1) x(2*i ) = pcent(2*i ) end do do i = 1,nc x(2*(np+i)-1) = ccent(2*i-1) x(2*(np+i) ) = ccent(2*i ) end do C Lower and upper bounds do i = 1,n l(i) = - 1.0d+20 u(i) = 1.0d+20 end do C Number of constraints (equalities plus inequalities) m = np + nc + totnvs + nc C Lagrange multipliers approximation. do i = 1,m lambda(i) = 0.0d0 end do C For each constraint i, set equatn(i) = .true. if it is an equality C constraint of the form c_i(x) = 0, and set equatn(i) = .false. if C it is an inequality constraint of the form c_i(x) <= 0. do i = 1,m equatn(i) = .false. end do C For each constraint i, set linear(i) = .true. if it is a linear C constraint, otherwise set linear(i) = .false. C Ellipse in and out do i = 1,np + nc linear(i) = .false. end do C Polygons do i = np + nc + 1,np + nc + totnvs linear(i) = .true. end do C Circles do i = np + nc + totnvs + 1,np + nc + totnvs + nc linear(i) = .false. end do end C ****************************************************************** C ****************************************************************** subroutine evalf(n,x,f,flag) implicit none C SCALAR ARGUMENTS integer flag,n double precision f C ARRAY ARGUMENTS double precision x(n) C This subroutine must compute the objective function. For achieving C this objective YOU MUST MODIFY it according to your problem. See C below where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C On Return C C f double precision, C objective function value at x, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the objective C function. (For example, trying to compute the square root C of a negative number, dividing by zero or a very small C number, etc.) If everything was o.k. you must set it C equal to zero. C LOCAL SCALARS integer i,ndist double precision diff1,diff2,dist C INTRINSIC FUNCTIONS intrinsic sqrt flag = 0 f = 0.0d0 ndist = n / 2 - 1 do i = 1,ndist diff1 = x(1) - x(2 * i + 1) diff2 = x(2) - x(2 * i + 2) dist = sqrt( diff1 ** 2 + diff2 ** 2) if ( dist .le. 1.0d-4 ) then write (*,fmt=*) + 'ERROR IN PROBLEM DEFINITION (DIST TOO SMALL)' flag = - 1 end if f = f + dist end do f = f / ndist end C ****************************************************************** C ****************************************************************** subroutine evalg(n,x,g,flag) implicit none C SCALAR ARGUMENTS integer flag,n C ARRAY ARGUMENTS double precision g(n),x(n) C This subroutine must compute the gradient vector of the objective C function. For achieving these objective YOU MUST MODIFY it in the C way specified below. However, if you decide to use numerical C derivatives (we dont encourage this option at all!) you dont need C to modify evalg. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C On Return C C g double precision g(n), C gradient vector of the objective function evaluated at x, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of any component C of the gradient vector. (For example, trying to compute C the square root of a negative number, dividing by zero or C a very small number, etc.) If everything was o.k. you C must set it equal to zero. C LOCAL SCALARS integer i,ndist double precision diff1,diff2,dist C INTRINSIC FUNCTIONS intrinsic sqrt flag = 0 g(1) = 0.0d0 g(2) = 0.0d0 ndist = n / 2 - 1 do i = 1,ndist diff1 = x(1) - x(2 * i + 1) diff2 = x(2) - x(2 * i + 2) dist = sqrt( diff1 ** 2 + diff2 ** 2) g(2*i+1) = - ( diff1 / dist ) / ndist g(2*i+2) = - ( diff2 / dist ) / ndist g(1) = g(1) - g(2 * i + 1) g(2) = g(2) - g(2 * i + 2) end do end C ****************************************************************** C ****************************************************************** subroutine evalh(n,x,hlin,hcol,hval,hnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,n,hnnz C ARRAY ARGUMENTS integer hcol(*),hlin(*) double precision hval(*),x(n) C This subroutine might compute the Hessian matrix of the objective C function. For achieving this objective YOU MAY MODIFY it according C to your problem. To modify this subroutine IS NOT MANDATORY. See C below where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C On Return C C hnnz integer, C number of perhaps-non-null elements of the computed C Hessian, C C hlin integer hlin(hnnz), C see below, C C hcol integer hcol(hnnz), C see below, C C hval double precision hval(hnnz), C the non-null value of the (hlin(k),hcol(k)) position C of the Hessian matrix of the objective function must C be saved at hval(k). Just the lower triangular part of C Hessian matrix must be computed, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the Hessian C matrix of the objective funtion. (For example, trying C to compute the square root of a negative number, C dividing by zero or a very small number, etc.) If C everything was o.k. you must set it equal to zero. C LOCAL SCALARS integer i,ndist double precision diff1,diff2,dist,g11,g22,g21 C INTRINSIC FUNCTIONS intrinsic sqrt flag = 0 ndist = n / 2 - 1 hnnz = 0 do i = 1,ndist diff1 = x(1) - x(2 * i + 1) diff2 = x(2) - x(2 * i + 2) dist = sqrt( diff1 ** 2 + diff2 ** 2 ) g11 = - ( - 1.0d0 + diff1 * diff1 / dist ** 2 ) / dist / ndist g22 = - ( - 1.0d0 + diff2 * diff2 / dist ** 2 ) / dist / ndist g21 = - ( + diff1 * diff2 / dist ** 2 ) / dist / ndist hlin(hnnz+1) = 1 hcol(hnnz+1) = 1 hval(hnnz+1) = g11 hlin(hnnz+2) = 2 * i + 1 hcol(hnnz+2) = 2 * i + 1 hval(hnnz+2) = g11 hlin(hnnz+3) = 2 * i + 1 hcol(hnnz+3) = 1 hval(hnnz+3) = - g11 hlin(hnnz+4) = 2 hcol(hnnz+4) = 2 hval(hnnz+4) = g22 hlin(hnnz+5) = 2 * i + 2 hcol(hnnz+5) = 2 * i + 2 hval(hnnz+5) = g22 hlin(hnnz+6) = 2 * i + 2 hcol(hnnz+6) = 2 hval(hnnz+6) = - g22 hlin(hnnz+7) = 2 hcol(hnnz+7) = 1 hval(hnnz+7) = g21 hlin(hnnz+8) = 2 * i + 2 hcol(hnnz+8) = 2 * i + 1 hval(hnnz+8) = g21 hlin(hnnz+9) = 2 * i + 1 hcol(hnnz+9) = 2 hval(hnnz+9) = - g21 hlin(hnnz+10) = 2 * i + 2 hcol(hnnz+10) = 1 hval(hnnz+10) = - g21 hnnz = hnnz + 10 end do end C ****************************************************************** C ****************************************************************** subroutine evalc(n,x,ind,c,flag) implicit none C SCALAR ARGUMENTS integer ind,flag,n double precision c C ARRAY ARGUMENTS double precision x(n) C This subroutine computes the ind-th constraint. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C ind integer, C index of the constraint to be computed, C C On Return C C c double precision, C ind-th constraint evaluated at x, C C flag integer C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C constraint. (For example, trying to compute the square C root of a negative number, dividing by zero or a very C small number, etc.) If everything was o.k. you must set C it equal to zero. C PARAMETERS integer mmax,ncmax,nmax,npmax,nvsmax parameter ( ncmax = 1000 ) parameter ( npmax = 1000 ) parameter ( nvsmax = 13 * npmax ) parameter ( nmax = 2 * ( npmax + ncmax ) ) parameter ( mmax = npmax + ncmax ) C COMMON SCALARS integer nc,np,totnvs double precision xdisp,xscal,ydisp,yscal C COMMON ARRAYS integer nvs(npmax),pol(nvsmax) double precision ccent(2*ncmax),edges(3*nvsmax),pcent(2*npmax), + radii(ncmax),vert(2*nvsmax) C LOCAL SCALARS integer i,j C LOCAL ARRAYS double precision p(2) 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 flag = 0 C Ellipse_in if ( ind .eq. 1 ) then p(1) = ( x(2 * ind - 1) - xdisp ) / xscal p(2) = ( x(2 * ind) - ydisp ) / yscal c = p(1) ** 2 + p(2) ** 2 - 1.0d0 C Ellipse_out else if ( ind .le. np + nc ) then p(1) = ( x(2 * ind - 1) - xdisp ) / xscal p(2) = ( x(2 * ind) - ydisp ) / yscal c = 1.0d0 - p(1) ** 2 - p(2) ** 2 C Polygons else if ( ind .le. np + nc + totnvs ) then j = ind - ( np + nc ) i = pol(j) c = edges(3*j-2) * x(2*i-1) + edges(3*j-1) * x(2*i) + + edges(3*j) C Circles else ! if ( ind .le. no + nc + totnvs + nc ) then j = ind - ( np + nc + totnvs ) i = j + np c = ( x(2*i-1) - ccent(2*j-1) ) ** 2 + + ( x(2*i ) - ccent(2*j ) ) ** 2 - radii(j) ** 2 end if end C ****************************************************************** C ****************************************************************** subroutine evaljac(n,x,ind,jcvar,jcval,jcnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,ind,n,jcnnz C ARRAY ARGUMENTS integer jcvar(n) double precision x(n),jcval(n) C This subroutine computes the gradient of the ind-th constraint. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C ind integer, C index of the constraint whose gradient will be computed, C C On Return C C jcnnz integer, C number of perhaps-non-null elements of the computed C gradient, C C jcvar integer jcvar(jcnnz), C see below, C C jcval double precision jcval(jcnnz), C the non-null value of the partial derivative of the C ind-th constraint with respect to the jcvar(k)-th C variable must be saved at jcval(k). C C flag integer C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C constraint. (For example, trying to compute the square C root of a negative number, dividing by zero or a very C small number, etc.) If everything was o.k. you must set C it equal to zero. C PARAMETERS integer mmax,ncmax,nmax,npmax,nvsmax parameter ( ncmax = 1000 ) parameter ( npmax = 1000 ) parameter ( nvsmax = 13 * npmax ) parameter ( nmax = 2 * ( npmax + ncmax ) ) parameter ( mmax = npmax + ncmax ) C COMMON SCALARS integer nc,np,totnvs double precision xdisp,xscal,ydisp,yscal C COMMON ARRAYS integer nvs(npmax),pol(nvsmax) double precision ccent(2*ncmax),edges(3*nvsmax),pcent(2*npmax), + radii(ncmax),vert(2*nvsmax) C LOCAL SCALARS integer i,j C LOCAL ARRAYS double precision p(2) 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 flag = 0 C Ellipse_in if ( ind .eq. 1 ) then p(1) = ( x(2 * ind - 1) - xdisp ) / xscal p(2) = ( x(2 * ind) - ydisp ) / yscal jcnnz = 2 jcvar(1) = 2 * ind - 1 jcval(1) = 2.0d0 * p(1) / xscal jcvar(2) = 2 * ind jcval(2) = 2.0d0 * p(2) / yscal C Ellipse_out else if ( ind .le. np + nc ) then p(1) = ( x(2 * ind - 1) - xdisp ) / xscal p(2) = ( x(2 * ind) - ydisp ) / yscal jcnnz = 2 jcvar(1) = 2 * ind - 1 jcval(1) = - 2.0d0 * p(1) / xscal jcvar(2) = 2 * ind jcval(2) = - 2.0d0 * p(2) / yscal C Polygons else if ( ind .le. np + nc + totnvs ) then j = ind - ( np + nc ) i = pol(j) jcnnz = 2 jcvar(1) = 2 * i - 1 jcval(1) = edges(3*j-2) jcvar(2) = 2 * i jcval(2) = edges(3*j-1) C Circles else ! if ( ind .le. no + nc + totnvs + nc ) then j = ind - ( np + nc + totnvs ) i = j + np jcnnz = 2 jcvar(1) = 2 * i - 1 jcval(1) = 2.0d0 * ( x(2*i-1) - ccent(2*j-1) ) jcvar(2) = 2 * i jcval(2) = 2.0d0 * ( x(2*i ) - ccent(2*j ) ) end if end C ****************************************************************** C ****************************************************************** subroutine evalhc(n,x,ind,hclin,hccol,hcval,hcnnz,flag) implicit none C SCALAR ARGUMENTS integer flag,ind,n,hcnnz C ARRAY ARGUMENTS integer hccol(*),hclin(*) double precision hcval(*),x(n) C This subroutine might compute the Hessian matrix of the ind-th C constraint. For achieving this objective YOU MAY MODIFY it C according to your problem. To modify this subroutine IS NOT C MANDATORY. See below where your modifications must be inserted. C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C ind integer, C index of the constraint whose Hessian will be computed, C C On Return C C hcnnz integer, C number of perhaps-non-null elements of the computed C Hessian, C C hclin integer hclin(hcnnz), C see below, C C hccol integer hccol(hcnnz), C see below, C C hcval double precision hcval(hcnnz), C the non-null value of the (hclin(k),hccol(k)) position C of the Hessian matrix of the ind-th constraint must C be saved at hcval(k). Just the lower triangular part of C Hessian matrix must be computed, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the Hessian C matrix of the ind-th constraint. (For example, trying C to compute the square root of a negative number, C dividing by zero or a very small number, etc.) If C everything was o.k. you must set it equal to zero. C PARAMETERS integer mmax,ncmax,nmax,npmax,nvsmax parameter ( ncmax = 1000 ) parameter ( npmax = 1000 ) parameter ( nvsmax = 13 * npmax ) parameter ( nmax = 2 * ( npmax + ncmax ) ) parameter ( mmax = npmax + ncmax ) C COMMON SCALARS integer nc,np,totnvs double precision xdisp,xscal,ydisp,yscal C COMMON ARRAYS integer nvs(npmax),pol(nvsmax) double precision ccent(2*ncmax),edges(3*nvsmax),pcent(2*npmax), + radii(ncmax),vert(2*nvsmax) C LOCAL SCALARS integer i,j 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 flag = 0 C Ellipse_in if ( ind .eq. 1 ) then hcnnz = 2 hclin(1) = 2 * ind - 1 hccol(1) = 2 * ind - 1 hcval(1) = 2.0d0 / xscal ** 2 hclin(2) = 2 * ind hccol(2) = 2 * ind hcval(2) = 2.0d0 / yscal ** 2 C Ellipse_out else if ( ind .le. np + nc ) then hcnnz = 2 hclin(1) = 2 * ind - 1 hccol(1) = 2 * ind - 1 hcval(1) = - 2.0d0 / xscal ** 2 hclin(2) = 2 * ind hccol(2) = 2 * ind hcval(2) = - 2.0d0 / yscal ** 2 C Polygons else if ( ind .le. np + nc + totnvs ) then hcnnz = 0 C Circles else ! if ( ind .le. no + nc + totnvs + nc ) then j = ind - ( np + nc + totnvs ) i = j + np hcnnz = 2 hclin(1) = 2 * i - 1 hccol(1) = 2 * i - 1 hcval(1) = 2.0d0 hclin(2) = 2 * i hccol(2) = 2 * i hcval(2) = 2.0d0 end if end C ****************************************************************** C ****************************************************************** subroutine evalhlp(n,x,m,lambda,p,hp,goth,flag) implicit none C SCALAR ARGUMENTS logical goth integer flag,m,n C ARRAY ARGUMENTS double precision hp(n),lambda(m),p(n),x(n) C This subroutine might compute the product of the Hessian of the C Lagrangian times vector p (just the Hessian of the objective C function in the unconstrained or bound-constrained case). C C Parameters of the subroutine: C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C m integer, C number of constraints, C C lambda double precision lambda(m), C vector of Lagrange multipliers, C C p double precision p(n), C vector of the matrix-vector product, C C goth logical, C can be used to indicate if the Hessian matrices were C computed at the current point. It is set to .false. C by the optimization method every time the current C point is modified. Sugestion: if its value is .false. C then compute the Hessians, save them in a common C structure and set goth to .true.. Otherwise, just use C the Hessians saved in the common block structure, C C On Return C C hp double precision hp(n), C Hessian-vector product, C C goth logical, C see above, C C flag integer, C You must set it to any number different of 0 (zero) if C some error ocurred during the evaluation of the C Hessian-vector product. (For example, trying to compute C the square root of a negative number, dividing by zero C or a very small number, etc.) If everything was o.k. you C must set it equal to zero. flag = - 1 end C ****************************************************************** C ****************************************************************** subroutine endp(n,x,l,u,m,lambda,equatn,linear) implicit none C SCALAR ARGUMENTS integer m,n C ARRAY ARGUMENTS logical equatn(m),linear(m) double precision l(n),lambda(m),u(n),x(n) C This subroutine can be used to do some extra job after the solver C has found the solution,like some extra statistics, or to save the C solution in some special format or to draw some graphical C representation of the solution. If the information given by the C solver is enough for you then leave the body of this subroutine C empty. C C Parameters of the subroutine: C C The paraemters of this subroutine are the same parameters of C subroutine inip. But in this subroutine there are not output C parameter. All the parameters are input parameters. call drawsol(n,x) end 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 = 1000 ) integer npmax parameter ( npmax = 1000 ) 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),pol(nvsmax) double precision ccent(2*ncmax),edges(3*nvsmax),pcent(2*npmax), + 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 = 1000 ) 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 = 1000 ) 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),pol(nvsmax) double precision edges(3*nvsmax),pcent(2*npmax),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 C ***************************************************************** C ***************************************************************** subroutine drawsol(n,x) C SCALAR ARGUMENTS integer n C ARRAY ARGUMENTS double precision x(n) C This subroutine generate a metapost file with the C graphical representation of the problem. C On Entry: C C n integer, C number of elements of array x, C C x double precision x(n) C current approximation to the solution. C PARAMETERS integer mmax,ncmax,nmax,npmax,nvsmax parameter ( ncmax = 1000 ) parameter ( npmax = 1000 ) parameter ( nvsmax = 13 * npmax ) parameter ( nmax = 2 * ( npmax + ncmax ) ) parameter ( mmax = npmax + ncmax ) C COMMON SCALARS integer nc,np,nx,ny,totnvs double precision rmax,xdisp,xscal,xstep,ydisp,yscal,ystep C COMMON ARRAYS integer nvs(npmax),pol(nvsmax) double precision ccent(2*ncmax),edges(3*nvsmax),pcent(2*npmax), + radii(ncmax),vert(2*nvsmax) C LOCAL SCALARS integer base,i,j double precision scale,vx,vy C COMMON BLOCKS common /const/ nx,ny,xstep,ystep,rmax 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 scale = min( 21.6d0 / ( ( nx - 1 ) * xstep + 2.0d0 * rmax ), + 27.9d0 / ( ( ny - 1 ) * ystep + 2.0d0 * rmax ) ) C DRAW THE PROBLEM open(unit=10,file='solution.mp') write(10,10) scale C ELLIPSE write(10,50) 2.0d0 * xscal, 2.0d0 * yscal,xdisp,ydisp C CIRCLES do i = 1,nc write(10,50) 2.0d0 * radii(i),2.0d0 * radii(i), + ccent(2 * i - 1),ccent(2 * i) end do C POLYGONS base = 0 do i = 1,np if ( nvs(i) .ne. 0 ) then write(10,20) do j = base + 1,base + nvs(i) vx = vert(2 * j - 1) vy = vert(2 * j) write(10,60) vx,vy end do write(10,40) base = base + nvs(i) end if end do C LINES do i = 1,np + nc write(10,70) x(1),x(2),x(2 * i - 1),x(2 * i) end do C DOTS WITH LABELS C write(10,80) C do i = 1,np + nc C write(10,90) i,x(2 * i - 1),x(2 * i) C end do C DOTS WITHOUT LABELS write(10,80) do i = 1,np + nc write(10,95) x(2 * i - 1),x(2 * i) end do C END write(10,100) close(10) C NON-EXECUTABLE STATEMENTS 10 format('beginfig(1);'/,'u = ',f20.10,' cm;') 20 format('draw') 30 format(5X,'(',f20.10,'u,',f20.10,'u)..') 40 format(5X,'cycle;') 50 format('draw fullcircle', + /,5X,'xscaled ',f20.10,'u', + /,5X,'yscaled ',f20.10,'u', + /,5X,'shifted (',f20.10,'u,',f20.10,'u);') 60 format(5X,'(',f20.10,'u,',f20.10,'u)--') 70 format('draw (',f20.10,'u,',f20.10,'u)--', + /,5X,'(',f20.10,'u,',f20.10,'u) dashed evenly;') 80 format('pickup pencircle scaled 4pt;') 90 format('dotlabel.top (btex $',i4,'$ etex, (',f20.10,'u,', + f20.10,'u) );') 95 format('draw (',f20.10,'u,',f20.10,'u);') 100 format('endfig;'/,'end;') end C ****************************************************************** C ****************************************************************** c double precision function drand(ix) c double precision ix c double precision a,p,b15,b16,xhi,xalo,leftlo,fhi,k c data a/16807.d0/,b15/32768.d0/,b16/65536.d0/,p/2147483647.d0/ c xhi= ix/b16 c xhi= xhi - dmod(xhi,1.d0) c xalo= (ix-xhi*b16)*a c leftlo= xalo/b16 c leftlo= leftlo - dmod(leftlo,1.d0) c fhi= xhi*a + leftlo c k= fhi/b15 c k= k - dmod(k,1.d0) c ix= (((xalo-leftlo*b16)-p)+(fhi-k*b15)*b16)+k c if (ix.lt.0) ix= ix + p c drand= ix*4.656612875d-10 c return c end