module ProblemType logical, private :: rectangle contains subroutine defRectanglePackingProblem() rectangle = .true. end subroutine defRectanglePackingProblem subroutine defCirclePackingProblem() rectangle = .false. end subroutine defCirclePackingProblem logical function isRectanglePackingProblem() isRectanglePackingProblem = rectangle end function isRectanglePackingProblem end module ProblemType ! ****************************************************************** ! ****************************************************************** subroutine inip(n,x,l,u,m,lambda,equatn,linear,coded,checkder,& nbv,niv,variableTypes,prob_id,num_rec) use ProblemType use Random implicit none integer mmax,nmax,ncirmax parameter ( mmax = 500000 ) parameter ( nmax = 500000 ) parameter ( ncirmax = 100 ) ! SCALAR ARGUMENTS integer m,n,nbv,niv,prob_id,num_rec logical checkder ! ARRAY ARGUMENTS logical :: coded(10) logical :: equatn(mmax),linear(mmax) double precision :: l(nmax),lambda(mmax),u(nmax),x(nmax) integer :: variableTypes(nmax) ! COMMON SCALARS double precision dx,dy,xmin,xmax,ymin,ymax integer probid,nrec,r,ncir double precision circle(ncirmax,3) double precision PI parameter ( PI = 3.141592653589793d0 ) integer CONTINUOUS,BINARY,INTEGER parameter ( CONTINUOUS = 0 ) parameter ( BINARY = 3 ) parameter ( INTEGER = 5 ) ! LOCAL SCALARS integer i logical, save :: problemInfoPrinted = .false. double precision lx,ux,ly,uy,objare,iteare,seed integer nreclb,nrecub ! COMMON BLOCKS common /packdata/ dx,dy,xmin,xmax,ymin,ymax,probid,nrec,r common /packcircledata/ circle,ncir ! Set problem data probid = prob_id nrec = num_rec call defpro(probid,dx,dy,lx,ux,ly,uy,xmin,xmax,ymin,ymax,r,& objare,iteare,nreclb,nrecub) if(.not. problemInfoPrinted) then write(*,100) probid,objare,iteare,nreclb,nrecub problemInfoPrinted = .true. end if call defdat(dx,dy,ncir,circle) call defRectanglePackingProblem() ! Number of variables n = 3 * nrec + 1 nbv = nrec niv = 0 ! Initial point seed = 7877.0d0 do i = 1,nrec x(2*i-1) = drand(seed) * (xmax - xmin) + xmin x(2*i ) = drand(seed) * (ymax - ymin) + ymin variableTypes(2*i-1) = CONTINUOUS variableTypes(2*i ) = CONTINUOUS end do do i = 1,nrec x(2*nrec + i) = drand(seed) variableTypes(2*nrec + i) = BINARY end do x(n) = drand(seed) * PI variableTypes(n) = CONTINUOUS ! Lower and upper bounds do i = 1,nrec l(2*i-1) = lx + min(dx,dy) u(2*i-1) = ux - min(dx,dy) l(2*i ) = ly + min(dx,dy) u(2*i ) = uy - min(dx,dy) end do do i = 1,nrec l(2*nrec + i) = 0.0d0 u(2*nrec + i) = 1.0d0 end do ! Theta l(n) = - 1.0d+20 u(n) = 1.0d+20 ! Number of constraints (equalities plus inequalities) m = 0 ! Lagrange multipliers approximation. do i = 1,m lambda(i) = 0.0d0 end do ! For each constraint i, set equatn(i) = .true. if it is an equality ! constraint of the form c_i(x) = 0, and set equatn(i) = .false. if ! it is an inequality constraint of the form c_i(x) <= 0. do i = 1,m equatn(i) = .true. end do ! For each constraint i, set linear(i) = .true. if it is a linear ! constraint, otherwise set linear(i) = .false. do i = 1,m linear(i) = .false. end do ! Indicate which subroutines did you code. coded( 1) = .true. ! evalf coded( 2) = .true. ! evalg coded( 3) = .false. ! evalh coded( 4) = .false. ! evalc coded( 5) = .false. ! evaljac coded( 6) = .false. ! evalhc coded( 7) = .false. ! evalfc coded( 8) = .false. ! evalgjac coded( 9) = .false. ! evalhl coded(10) = .false. ! evalhlp ! Set checkder = .true. if you code some derivatives and you would ! like them to be tested by finite differences. It is highly ! recommended. checkder = .false. 100 format( ' Problem = ',I8, & /,' Estimated object area = ',F8.2,& /,' Item area = ',F8.2,& /,' Known lower bound = ',I8, & /,' Estimated upper bound = ',I8,/) end subroutine inip ! ****************************************************************** ! ****************************************************************** subroutine evalf(n,x,f,flag) use ProblemType implicit none ! SCALAR ARGUMENTS integer flag,n double precision f ! ARRAY ARGUMENTS double precision x(n) ! COMMON SCALARS integer probid,nrec,r double precision dx,dy,xmin,xmax,ymin,ymax ! COMMON BLOCKS common /packdata/ dx,dy,xmin,xmax,ymin,ymax,probid,nrec,r if(isRectanglePackingProblem()) then call evalal(n,x,f,flag) else call evalal_circle(n,x,f,flag) end if end subroutine evalf ! ****************************************************************** ! ****************************************************************** subroutine evalw(probid,p,r,w) implicit none ! This subroutine computes, for each problem, the smooth functions ! that define the convex region. ! SCALAR ARGUMENTS integer probid,r ! ARRAY ARGUMENTS double precision p(2),w(r) ! LOCAL SCALARS double precision tmp if ( probid .eq. 1 ) then w(1) = - p(1) w(2) = - p(2) w(3) = - p(1) - p(2) + 3.0d0 w(4) = p(1) ** 2.0d0 + p(2) ** 2.0d0 - 100.0d0 else if ( probid .eq. 2 ) then w(1) = -7.0d0 * p(1) + 6.0d0 * p(2) - 24.0d0 w(2) = 7.0d0 * p(1) + 6.0d0 * p(2) - 108.0d0 w(3) = ( p(1) - 6.0d0 ) ** 2.0d0 + ( p(2) - 8.0d0 ) ** 2.0d0 - 9.0d0 else if ( probid .eq. 3 ) then w(1) = - p(1) w(2) = p(1) - 8.0d0 w(3) = ( p(1) - 6.0d0 ) ** 2.0d0 + p(2) ** 2.0d0 - 81.0d0 w(4) = ( p(1) - 1.7d0 ) ** 2.0d0 + ( p(2) - 10.0d0 ) ** 2.0d0 - 81.0d0 else if ( probid .eq. 4 ) then w(1) = p(1) ** 2.0d0 - p(2) w(2) = 0.25d0 * p(1) ** 2.0d0 + p(2) - 5.0d0 else if ( probid .eq. 5 ) then w(1) = p(1) ** 2.0d0 - p(2) w(2) = - p(1) + p(2) ** 2.0d0 - 6.0d0 * p(2) + 6.0d0 w(3) = p(1) + p(2) - 6.0d0 else if ( probid .eq. 6 ) then w(1) = - p(1) + p(2) ** 2.0d0 - 6.0d0 * p(2) + 6.0d0 w(2) = p(1) + p(2) ** 2.0d0 - 3.0d0 * p(2) - 0.75d0 else if ( probid .eq. 7 ) then w(1) = 0.25d0 * ( p(1) - 2.0d0 ) ** 2.0d0 + & 0.0625d0 * ( p(2) - 4.0d0 ) ** 2.0d0 - 1.0d0 else if ( probid .eq. 8 ) then w(1) = 0.25d0 * ( p(1) - 6.0d0 ) ** 2.0d0 + & ( p(2) - 6.0d0 ) ** 2.0d0 / 36.0d0 - 1.0d0 w(2) = ( p(1) - 6.0d0 ) ** 2.0d0 / 36.0d0 + & 0.25d0 * ( p(2) - 6.0d0 ) ** 2.0d0 - 1.0d0 w(3) = p(1) - p(2) - 3.0d0 w(4) = -p(1) + p(2) - 2.0d0 else if ( probid .eq. 9 ) then w(1) = 0.25d0 * ( p(1) - 3.0d0 ) ** 2.0d0 + & 0.0625d0 * ( p(2) - 4.0d0 ) ** 2.0d0 - 1.0d0 w(2) = 0.25d0 * ( p(1) - 2.65d0 ) ** 2.0d0 + & 0.0625d0 * ( p(2) - 4.0d0 ) ** 2.0d0 - 1.0d0 w(3) = - p(1) + 1.0d0 w(4) = p(1) - p(2) - 1.0d0 w(5) = p(1) + p(2) - 9.0d0 else if ( probid .eq. 10 ) then w(1) = ( p(1) - 6.0d0 ) ** 2.0d0 / 36.0d0 + & 0.25d0 * ( p(2) - 6.0d0 ) ** 2.0d0 - 1.0d0 w(2) = ( p(1) - 6.0d0 ) ** 2.0d0 / 9.0d0 + & ( p(2) - 8.0d0 ) ** 2.0d0 / 9.0d0 - 1.0d0 else if ( probid .eq. 11 ) then w(1) = ( p(1) / 6.0d0 ) ** 4 + ( p(2) / 2.0d0 ) ** 4 - 1.0d0 w(2) = 8.0d0 * p(1) - 11.0d0 * p(2) - 26.0d0 else if ( probid .ge. 12 .and. probid .le. 21 ) then if ( probid .eq. 12 .or. probid .eq. 17 ) then tmp = 8.619d0 else if ( probid .eq. 13 .or. probid .eq. 18 ) then tmp = 8.774d0 else if ( probid .eq. 14 .or. probid .eq. 19 ) then tmp = 9.155d0 else if ( probid .eq. 15 .or. probid .eq. 20 ) then tmp = 9.302d0 else if ( probid .eq. 16 .or. probid .eq. 21 ) then tmp = 9.310d0 else tmp = 9.310d0 write(*,*) 'Undefined problem!' stop end if w(1) = sqrt( 3.0d0 ) * p(1) + p(2) - sqrt( 3.0d0 ) * tmp w(2) = - sqrt( 3.0d0 ) * p(1) + p(2) w(3) = - p(2) end if end subroutine evalw ! ****************************************************************** ! ****************************************************************** subroutine evaldw(probid,p,r,dwdp1,dwdp2) implicit none ! This subroutine computes, for each problem, the derivatives of the ! smooth functions that define the convex region. ! SCALAR ARGUMENTS integer probid,r ! ARRAY ARGUMENTS double precision p(2) double precision dwdp1(r),dwdp2(r) if ( probid .eq. 1 ) then dwdp1(1) = - 1.0d0 dwdp1(2) = 0.0d0 dwdp1(3) = - 1.0d0 dwdp1(4) = 2.0d0 * p(1) dwdp2(1) = 0.0d0 dwdp2(2) = - 1.0d0 dwdp2(3) = - 1.0d0 dwdp2(4) = 2.0d0 * p(2) else if ( probid .eq. 2 ) then dwdp1(1) = - 7.0d0 dwdp1(2) = 7.0d0 dwdp1(3) = 2.0d0 * ( p(1) - 6.0d0 ) dwdp2(1) = 6.0d0 dwdp2(2) = 6.0d0 dwdp2(3) = 2.0d0 * ( p(2) - 8.0d0 ) else if ( probid .eq. 3 ) then dwdp1(1) = - 1.0d0 dwdp1(2) = 1.0d0 dwdp1(3) = 2.0d0 * ( p(1) - 6.0d0 ) dwdp1(4) = 2.0d0 * ( p(1) - 1.7d0 ) dwdp2(1) = 0.0d0 dwdp2(2) = 0.0d0 dwdp2(3) = 2.0d0 * ( p(2) - 0.0d0 ) dwdp2(4) = 2.0d0 * ( p(2) - 10.0d0 ) else if ( probid .eq. 4 ) then dwdp1(1) = 2.0d0 * p(1) dwdp1(2) = 0.5d0 * p(1) dwdp2(1) = - 1.0d0 dwdp2(2) = 1.0d0 else if ( probid .eq. 5 ) then dwdp1(1) = 2.0d0 * p(1) dwdp1(2) = - 1.0d0 dwdp1(3) = 1.0d0 dwdp2(1) = - 1.0d0 dwdp2(2) = 2.0d0 * p(2) - 6.0d0 dwdp2(3) = 1.0d0 else if ( probid .eq. 6 ) then dwdp1(1) = - 1.0d0 dwdp1(2) = 1.0d0 dwdp2(1) = 2.0d0 * p(2) - 6.0d0 dwdp2(2) = 2.0d0 * p(2) - 3.0d0 else if ( probid .eq. 7 ) then dwdp1(1) = 0.5d0 * ( p(1) - 2.0d0 ) dwdp2(1) = 0.125d0 * ( p(2) - 4.0d0 ) else if ( probid .eq. 8 ) then dwdp1(1) = 0.5d0 * ( p(1) - 6.0d0 ) dwdp1(2) = ( p(1) - 6.0d0 ) / 18.0d0 dwdp1(3) = 1.0d0 dwdp1(4) = -1.0d0 dwdp2(1) = ( p(2) - 6.0d0 ) / 18.0d0 dwdp2(2) = 0.5d0 * ( p(2) - 6.0d0 ) dwdp2(3) = -1.0d0 dwdp2(4) = 1.0d0 else if ( probid .eq. 9 ) then dwdp1(1) = 0.5d0 * ( p(1) - 3.0d0 ) dwdp1(2) = 0.5d0 * ( p(1) - 2.65d0 ) dwdp1(3) = -1.0d0 dwdp1(4) = 1.0d0 dwdp1(5) = 1.0d0 dwdp2(1) = 0.125d0 * ( p(2) - 4.0d0 ) dwdp2(2) = 0.125d0 * ( p(2) - 4.0d0 ) dwdp2(3) = 0.0d0 dwdp2(4) = -1.0d0 dwdp2(5) = 1.0d0 else if ( probid .eq. 10 ) then dwdp1(1) = ( p(1) - 6.0d0 ) / 18.0d0 dwdp1(2) = 2.0d0 * ( p(1) - 6.0d0 ) / 9.0d0 dwdp2(1) = 0.5d0 * ( p(2) - 6.0d0 ) dwdp2(2) = 2.0d0 * ( p(2) - 8.0d0 ) / 9.0d0 else if ( probid .eq. 11 ) then dwdp1(1) = 2.0d0 * ( p(1) / 6.0d0 ) ** 3.0d0 / 3.0d0 dwdp1(2) = 8.0d0 dwdp2(1) = 2.0d0 * ( p(2) / 2.0d0 ) ** 3.0d0 dwdp2(2) = - 11.0d0 else if ( probid .ge. 12 .and. probid .le. 21 ) then dwdp1(1) = sqrt( 3.0d0 ) dwdp1(2) = -sqrt( 3.0d0 ) dwdp1(3) = 0.0d0 dwdp2(1) = 1.0d0 dwdp2(2) = 1.0d0 dwdp2(3) = -1.0d0 end if end subroutine evaldw ! ****************************************************************** ! ****************************************************************** subroutine evalal(n,x,f,flag) implicit none ! On Entry: ! ! n integer ! number of variables ! ! x double precision x(n) ! current point ! ! NOTE: arguments m, lambda and rho are useful when GENCAN is ! being used for solving the box-constrained subproblems of an ! Augmented Lagrangian framework. When GENCAN is being used ! stand-alone for solving a bound-constrained problem, these ! arguments are dummy arguments and must be ignored. ! ! On Return ! ! f double precision ! objective function value at x ! ! flag integer ! 0 means ''no errors'' ! any other value means ''there was an error in the ! objective function calculation''. ! ! PARAMETERS integer rmax parameter ( rmax = 1000 ) ! COMMON SCALARS integer probid,nrec,r double precision dx,dy,xmin,xmax,ymin,ymax ! SCALAR ARGUMENTS integer flag,n double precision f ! ARRAY ARGUMENTS double precision x(n) ! LOCAL SCALARS integer i,j double precision dxi,dxj,dyi,dyj,tmpx,tmpy double precision theta,sin_t,cos_t double precision c1i,c2i,pi,c1j,c2j,pj ! LOCAL ARRAYS double precision q(2),w(rmax) ! COMMON BLOCKS common /packdata/ dx,dy,xmin,xmax,ymin,ymax,probid,nrec,r flag = 0 f = 0.0d0 theta = x(n) sin_t = sin(theta) cos_t = cos(theta) do i = 1,nrec c1i = x(2*i-1) c2i = x(2*i ) pi = x(2*nrec + i) dxi = (pi*dx + (1.0d0 - pi)*dy) dyi = (pi*dy + (1.0d0 - pi)*dx) ! Rectangles overlapping do j = i + 1,nrec c1j = x(2*j-1) c2j = x(2*j ) pj = x(2*nrec + j) dxj = (pj*dx + (1.0d0 - pj)*dy) dyj = (pj*dy + (1.0d0 - pj)*dx) tmpx = (dxi + dxj)**2.0d0 - & (c1i*cos_t + c2i*sin_t - c1j*cos_t - c2j*sin_t)**2.0d0 tmpy = (dyi + dyj)**2.0d0 - & (- c1i*sin_t + c2i*cos_t + c1j*sin_t - c2j*cos_t)**2.0d0 f = f + max(0.0d0, tmpx) ** 2.0d0 * max(0.0d0, tmpy) ** 2.0d0 end do ! Region constraints ! SW corner q(1) = - dxi * cos_t + dyi * sin_t + c1i q(2) = - dxi * sin_t - dyi * cos_t + c2i call evalw(probid,q,r,w) do j = 1,r f = f + max( 0.0d0, w(j) ) ** 2 end do ! SE corner q(1) = dxi * cos_t + dyi * sin_t + c1i q(2) = dxi * sin_t - dyi * cos_t + c2i call evalw(probid,q,r,w) do j = 1,r f = f + max( 0.0d0, w(j) ) ** 2 end do ! NE corner q(1) = dxi * cos_t - dyi * sin_t + c1i q(2) = dxi * sin_t + dyi * cos_t + c2i call evalw(probid,q,r,w) do j = 1,r f = f + max( 0.0d0, w(j) ) ** 2 end do ! NW corner q(1) = - dxi * cos_t - dyi * sin_t + c1i q(2) = - dxi * sin_t + dyi * cos_t + c2i call evalw(probid,q,r,w) do j = 1,r f = f + max( 0.0d0, w(j) ) ** 2 end do end do end subroutine evalal ! ****************************************************************** ! ****************************************************************** subroutine evalnal(n,x,g,flag) implicit none ! On Entry: ! n integer ! number of variables ! ! x double precision x(n) ! current point ! ! NOTE: arguments m, lambda and rho are useful when GENCAN is ! being used for solving the box-constrained subproblems of an ! Augmented Lagrangian framework. When GENCAN is being used ! stand-alone for solving a bound-constrained problem, these ! arguments are dummy arguments and must be ignored. ! ! On Return ! ! g double precision g(n) ! gradient of the objective function at x ! ! flag integer ! 0 means ''no errors'', ! any other value means ''there was an error in the ! gradient calculation''. ! ! PARAMETERS integer rmax parameter ( rmax = 1000 ) ! COMMON SCALARS integer probid,nrec,r double precision dx,dy,xmin,xmax,ymin,ymax ! SCALAR ARGUMENTS integer flag,n ! ARRAY ARGUMENTS double precision g(n),x(n) ! LOCAL SCALARS integer i,j double precision dxi,dxj,dyi,dyj,diffx,diffy,tmpx,tmpy,gradx,grady double precision basex,basey,gradp,tmp_t,dq1dv,dq2dv,dwdq double precision c1i,c2i,pi,c1j,c2j,pj double precision cos_t,sin_t,theta ! LOCAL ARRAYS double precision q(2),w(rmax) double precision partial_dwdq1(rmax),partial_dwdq2(rmax) ! COMMON BLOCKS common /packdata/ dx,dy,xmin,xmax,ymin,ymax,probid,nrec,r flag = 0 do i = 1,n g(i) = 0.0d0 end do theta = x(n) sin_t = sin(theta) cos_t = cos(theta) do i = 1,nrec c1i = x(2*i-1) c2i = x(2*i ) pi = x(2*nrec + i) dxi = (pi*dx + (1.0d0 - pi)*dy) dyi = (pi*dy + (1.0d0 - pi)*dx) ! Rectangles overlapping do j = i + 1,nrec c1j = x(2*j-1) c2j = x(2*j ) pj = x(2*nrec + j) dxj = (pj*dx + (1.0d0 - pj)*dy) dyj = (pj*dy + (1.0d0 - pj)*dx) diffx = c1i*cos_t + c2i*sin_t - c1j*cos_t - c2j*sin_t diffy = - c1i*sin_t + c2i*cos_t + c1j*sin_t - c2j*cos_t tmpx = ( dxi + dxj ) ** 2.0d0 - diffx ** 2.0d0 tmpy = ( dyi + dyj ) ** 2.0d0 - diffy ** 2.0d0 if( tmpx .gt. 0.0d0 .and. tmpy .gt. 0.0d0 ) then basex = - 4.0d0 * tmpx * diffx basey = - 4.0d0 * tmpy * diffy gradx = (basex * cos_t * (tmpy**2.0d0)) - & ((tmpx**2.0d0) * basey * sin_t) grady = (basex * sin_t * (tmpy**2.0d0)) + & ((tmpx**2.0d0) * basey * cos_t) ! C_i g(2*i-1) = g(2*i-1) + gradx g(2*i ) = g(2*i ) + grady ! C_j g(2*j-1) = g(2*j-1) - gradx g(2*j ) = g(2*j ) - grady ! p gradp = 4.0d0 *(tmpx * (dxi+dxj) * (dx-dy) * (tmpy**2.0d0) + & (tmpx**2.0d0) * tmpy * (dyi+dyj) * (dy-dx)) ! p_i g(2*nrec + i) = g(2*nrec + i) + gradp ! p_j g(2*nrec + j) = g(2*nrec + j) + gradp ! theta tmp_t = - c1i*sin_t + c2i*cos_t + c1j*sin_t - c2j*cos_t g(n) = g(n) - 4.0d0 * tmp_t * tmpx * diffx * (tmpy**2.0d0) tmp_t = - c1i*cos_t - c2i*sin_t + c1j*cos_t + c2j*sin_t g(n) = g(n) - 4.0d0 * tmp_t * tmpy * diffy * (tmpx**2.0d0) end if end do ! Region constraints ! SW corner q(1) = - dxi * cos_t + dyi * sin_t + c1i q(2) = - dxi * sin_t - dyi * cos_t + c2i call evalw(probid,q,r,w) call evaldw(probid,q,r,partial_dwdq1,partial_dwdq2) ! C do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) g(2*i-1) = g(2*i-1) + 2.0d0 * w(j) * dwdq dwdq = partial_dwdq2(j) g(2*i ) = g(2*i ) + 2.0d0 * w(j) * dwdq end if end do ! p_i dq1dv = - cos_t * (dx - dy) + sin_t * (dy - dx) dq2dv = - sin_t * (dx - dy) - cos_t * (dy - dx) do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(2*nrec + i) = g(2*nrec + i) + 2.0d0 * w(j) * dwdq end if end do ! theta dq1dv = dxi * sin_t + dyi * cos_t dq2dv = - dxi * cos_t + dyi * sin_t do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(n) = g(n) + 2.0d0 * w(j) * dwdq end if end do ! SE corner q(1) = dxi * cos_t + dyi * sin_t + c1i q(2) = dxi * sin_t - dyi * cos_t + c2i call evalw(probid,q,r,w) call evaldw(probid,q,r,partial_dwdq1,partial_dwdq2) ! C do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) g(2*i-1) = g(2*i-1) + 2.0d0 * w(j) * dwdq dwdq = partial_dwdq2(j) g(2*i ) = g(2*i ) + 2.0d0 * w(j) * dwdq end if end do ! p_i dq1dv = cos_t * (dx - dy) + sin_t * (dy - dx) dq2dv = sin_t * (dx - dy) - cos_t * (dy - dx) do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(2*nrec + i) = g(2*nrec + i) + 2.0d0 * w(j) * dwdq end if end do ! theta dq1dv = - dxi * sin_t + dyi * cos_t dq2dv = dxi * cos_t + dyi * sin_t do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(n) = g(n) + 2.0d0 * w(j) * dwdq end if end do ! NE corner q(1) = dxi * cos_t - dyi * sin_t + c1i q(2) = dxi * sin_t + dyi * cos_t + c2i call evalw(probid,q,r,w) call evaldw(probid,q,r,partial_dwdq1,partial_dwdq2) ! C do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) g(2*i-1) = g(2*i-1) + 2.0d0 * w(j) * dwdq dwdq = partial_dwdq2(j) g(2*i ) = g(2*i ) + 2.0d0 * w(j) * dwdq end if end do ! p_i dq1dv = cos_t * (dx - dy) - sin_t * (dy - dx) dq2dv = sin_t * (dx - dy) + cos_t * (dy - dx) do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(2*nrec + i) = g(2*nrec + i) + 2.0d0 * w(j) * dwdq end if end do ! theta dq1dv = - dxi * sin_t - dyi * cos_t dq2dv = dxi * cos_t - dyi * sin_t do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(n) = g(n) + 2.0d0 * w(j) * dwdq end if end do ! NW corner q(1) = - dxi * cos_t - dyi * sin_t + c1i q(2) = - dxi * sin_t + dyi * cos_t + c2i call evalw(probid,q,r,w) call evaldw(probid,q,r,partial_dwdq1,partial_dwdq2) ! C do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) g(2*i-1) = g(2*i-1) + 2.0d0 * w(j) * dwdq dwdq = partial_dwdq2(j) g(2*i ) = g(2*i ) + 2.0d0 * w(j) * dwdq end if end do ! p_i dq1dv = - cos_t * (dx - dy) - sin_t * (dy - dx) dq2dv = - sin_t * (dx - dy) + cos_t * (dy - dx) do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(2*nrec + i) = g(2*nrec + i) + 2.0d0 * w(j) * dwdq end if end do ! theta dq1dv = dxi * sin_t - dyi * cos_t dq2dv = - dxi * cos_t - dyi * sin_t do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(n) = g(n) + 2.0d0 * w(j) * dwdq end if end do end do end subroutine evalnal ! ****************************************************************** ! ****************************************************************** subroutine defpro(probid,dx,dy,lx,ux,ly,uy,xmin,xmax,ymin,ymax,r,& objare,iteare,nreclb,nrecub) implicit none ! This subroutine defines, for each problem, the items dimensions. ! It also defines the number of smooth functions and the box that ! define the convex region. ! PARAMETERS integer rmax parameter ( rmax = 1000 ) ! SCALARS ARGUMENTS integer nreclb,nrecub,probid,r double precision dx,dy,lx,ux,ly,uy,xmin,xmax,ymin,ymax,objare,iteare ! LOCAL SCALARS logical feas integer i,j,k,nintx,ninty double precision tmp,xstep,ystep,unita ! LOCAL ARRAYS double precision q(2),w(rmax) if ( probid .eq. 1 ) then dx = 0.5d0 * 2.0d0 dy = 0.5d0 * 1.0d0 lx = 0.0d0 ux = 1.0d+99 ly = 0.0d0 uy = 1.0d+99 xmin = 0.0d0 xmax = 10.0d0 ymin = 0.0d0 ymax = 10.0d0 r = 4 else if ( probid .eq. 2 ) then dx = 0.5d0 * 1.10d0 dy = 0.5d0 * 0.55d0 lx = - 1.0d+99 ux = 1.0d+99 ly = - 1.0d+99 uy = 1.0d+99 xmin = 3.0d0 xmax = 9.0d0 ymin = 5.0d0 ymax = 11.0d0 r = 3 else if ( probid .eq. 3 ) then dx = 0.5d0 * 2.0d0 dy = 0.5d0 * 0.6d0 lx = 0.0d0 ux = 8.0d0 ly = - 1.0d+99 uy = 1.0d+99 xmin = 0.0d0 xmax = 8.0d0 ymin = 1.0d0 ymax = 9.0d0 r = 4 else if ( probid .eq. 4 ) then dx = 0.5d0 * 1.0d0 dy = 0.5d0 * 0.4d0 lx = - 1.0d+99 ux = 1.0d+99 ly = - 1.0d+99 uy = 1.0d+99 xmin = -3.0d0 xmax = 3.0d0 ymin = 0.0d0 ymax = 5.0d0 r = 2 else if ( probid .eq. 5 ) then dx = 0.5d0 * 0.90d0 dy = 0.5d0 * 0.30d0 lx = - 1.0d+99 ux = 1.0d+99 ly = - 1.0d+99 uy = 1.0d+99 xmin = -3.0d0 xmax = 3.0d0 ymin = 0.0d0 ymax = 6.0d0 r = 3 else if ( probid .eq. 6 ) then dx = 0.5d0 * 0.90d0 dy = 0.5d0 * 0.30d0 lx = - 1.0d+99 ux = 1.0d+99 ly = - 1.0d+99 uy = 1.0d+99 xmin = -3.0d0 xmax = 3.0d0 ymin = 0.0d0 ymax = 6.0d0 r = 2 else if ( probid .eq. 7 ) then dx = 0.5d0 * 2.0d0 dy = 0.5d0 * 0.5d0 lx = - 1.0d+99 ux = 1.0d+99 ly = - 1.0d+99 uy = 1.0d+99 xmin = 0.0d0 xmax = 4.0d0 ymin = 0.0d0 ymax = 8.0d0 r = 1 else if ( probid .eq. 8 ) then dx = 0.5d0 * 0.7d0 dy = 0.5d0 * 0.5d0 lx = - 1.0d+99 ux = 1.0d+99 ly = - 1.0d+99 uy = 1.0d+99 xmin = 4.0d0 xmax = 8.0d0 ymin = 4.0d0 ymax = 8.0d0 r = 4 else if ( probid .eq. 9 ) then dx = 0.5d0 * 0.8d0 dy = 0.5d0 * 0.6d0 lx = 1.0d0 ux = 1.0d+99 ly = - 1.0d+99 uy = 1.0d+99 xmin = 1.00d0 xmax = 4.65d0 ymin = 0.00d0 ymax = 8.00d0 r = 5 else if ( probid .eq. 10 ) then dx = 0.5d0 * 0.95d0 dy = 0.5d0 * 0.35d0 lx = 0.0d0 ux = 1.0d+99 ly = - 1.0d+99 uy = 1.0d+99 xmin = 3.0d0 xmax = 9.0d0 ymin = 5.0d0 ymax = 8.0d0 r = 2 else if ( probid .eq. 11 ) then dx = 0.5d0 * 1.9d0 dy = 0.5d0 * 0.5d0 lx = - 7.0d0 ux = 7.0d0 ly = - 3.0d0 uy = 3.0d0 xmin = -6.0d0 xmax = 6.0d0 ymin = -2.0d0 ymax = 2.0d0 r = 2 else if ( probid .ge. 12 .and. probid .le. 21 ) then if ( probid .le. 16 ) then dx = 0.5d0 * 1.0d0 dy = 0.5d0 * 1.0d0 else dx = 0.5d0 * 1.0d0 dy = 0.5d0 * 0.5d0 end if lx = - 1.0d+99 ux = 1.0d+99 ly = 0.0d0 uy = 1.0d+99 if ( probid .eq. 12 .or. probid .eq. 17 ) then tmp = 8.619d0 else if ( probid .eq. 13 .or. probid .eq. 18 ) then tmp = 8.774d0 else if ( probid .eq. 14 .or. probid .eq. 19 ) then tmp = 9.155d0 else if ( probid .eq. 15 .or. probid .eq. 20 ) then tmp = 9.302d0 else if ( probid .eq. 16 .or. probid .eq. 21 ) then tmp = 9.310d0 else tmp = 8.619d0 write(*,*) 'Undefined problem!' stop end if xmin = 0.0d0 xmax = tmp ymin = 0.0d0 ymax = 0.5d0 * tmp * sqrt( 3.0d0 ) r = 3 else write(*,*) 'Undefined problem!' stop end if ! ESTIMATE OBJECT AREA nintx = 1000 ninty = 1000 xstep = ( xmax - xmin ) / real(nintx,kind=8) ystep = ( ymax - ymin ) / real(ninty,kind=8) unita = xstep * ystep objare = 0.0d0 do i = 1,nintx do j = 1,ninty q(1) = xmin + ( real(i,kind=8) - 0.5d0 ) * xstep q(2) = ymin + ( real(j,kind=8) - 0.5d0 ) * ystep call evalw(probid,q,r,w) feas = .true. do k = 1,r if ( w(k) .gt. 0.0d0 ) then feas = .false. end if end do if ( feas ) then objare = objare + unita end if end do end do ! COMPUTE ITEM AREA iteare = 4.0d0 * dx * dy ! SET TRIVIAL BOUNDS nreclb = 0 nrecub = int( objare / iteare ) end subroutine defpro ! ****************************************************************** ! ****************************************************************** subroutine evalg(n,x,g,flag) use ProblemType implicit none ! SCALAR ARGUMENTS integer flag,n ! ARRAY ARGUMENTS double precision g(n),x(n) flag = 0 if(isRectanglePackingProblem()) then call evalnal(n,x,g,flag) else call evalnal_circle(n,x,g,flag) end if end subroutine evalg ! ****************************************************************** ! ****************************************************************** subroutine evalh(n,x,hlin,hcol,hval,hnnz,flag) implicit none ! SCALAR ARGUMENTS integer flag,hnnz,n ! ARRAY ARGUMENTS integer hcol(*),hlin(*) double precision hval(*),x(n) flag = - 1 end subroutine evalh ! ****************************************************************** ! ****************************************************************** subroutine evalc(n,x,ind,c,flag) implicit none ! SCALAR ARGUMENTS integer ind,flag,n double precision c ! ARRAY ARGUMENTS double precision x(n) ! Constraints flag = - 1 end subroutine evalc ! ****************************************************************** ! ****************************************************************** subroutine evaljac(n,x,ind,jcvar,jcval,jcnnz,flag) implicit none ! SCALAR ARGUMENTS integer flag,ind,jcnnz,n ! ARRAY ARGUMENTS integer jcvar(n) double precision x(n),jcval(n) ! Sparse gradient vector of the ind-th constraint flag = - 1 end subroutine evaljac ! ****************************************************************** ! ****************************************************************** subroutine evalhc(n,x,ind,hclin,hccol,hcval,hcnnz,flag) implicit none ! SCALAR ARGUMENTS integer flag,hcnnz,ind,n ! ARRAY ARGUMENTS integer hccol(*),hclin(*) double precision hcval(*),x(n) flag = - 1 end subroutine evalhc ! ****************************************************************** ! ****************************************************************** subroutine evalfc(n,x,f,m,c,flag) implicit none ! SCALAR ARGUMENTS integer flag,m,n double precision f ! ARRAY ARGUMENTS double precision c(m),x(n) flag = - 1 end subroutine evalfc ! ****************************************************************** ! ****************************************************************** subroutine evalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,flag) implicit none ! SCALAR ARGUMENTS integer flag,jcnnz,m,n ! ARRAY ARGUMENTS integer jcfun(*),jcvar(*) double precision g(n),jcval(*),x(n) flag = - 1 end subroutine evalgjac ! ****************************************************************** ! ****************************************************************** subroutine evalhl(n,x,m,lambda,scalef,scalec,hllin,hlcol,hlval,hlnnz,flag) implicit none ! SCALAR ARGUMENTS integer flag,hlnnz,m,n double precision scalef ! ARRAY ARGUMENTS integer hlcol(*),hllin(*) double precision hlval(*),lambda(m),scalec(m),x(n) flag = - 1 end subroutine evalhl ! ****************************************************************** ! ****************************************************************** subroutine evalhlp(n,x,m,lambda,scalef,scalec,p,hp,goth,flag) implicit none ! SCALAR ARGUMENTS logical goth integer flag,m,n double precision scalef ! ARRAY ARGUMENTS double precision hp(n),lambda(m),p(n),scalec(m),x(n) flag = - 1 end subroutine evalhlp ! ****************************************************************** ! ****************************************************************** subroutine endp(n,x,l,u,m,lambda,equatn,linear) implicit none ! SCALAR ARGUMENTS integer m,n ! ARRAY ARGUMENTS logical :: equatn(m),linear(m) double precision :: l(n),lambda(m),u(n),x(n) end subroutine endp ! ****************************************************************** ! ****************************************************************** subroutine randomInitialPoint(n,x,seed) use Random implicit none ! SCALAR ARGUMENTS integer , intent(in) :: n real(kind=8), intent(in out) :: seed ! ARRAY ARGUMENTS real(kind=8), intent(out) :: x(n) double precision PI parameter ( PI = 3.141592653589793d0 ) ! COMMON SCALARS double precision dx,dy,xmin,xmax,ymin,ymax integer probid,nrec,r ! COMMON BLOCKS common /packdata/ dx,dy,xmin,xmax,ymin,ymax,probid,nrec,r ! LOCAL SCALARS integer :: i ! C do i = 1,nrec x(2*i - 1) = drand(seed) * (xmax - xmin) + xmin x(2*i ) = drand(seed) * (ymax - ymin) + ymin end do ! P do i = 1,nrec x(2*nrec + i) = drand(seed) end do ! THETA x(n) = drand(seed) * PI end subroutine randomInitialPoint ! ****************************************************************** ! ****************************************************************** subroutine evalal_circle(n,x,f,flag) implicit none ! On Entry: ! ! n integer ! number of variables ! ! x double precision x(n) ! current point ! ! NOTE: arguments m, lambda and rho are useful when GENCAN is ! being used for solving the box-constrained subproblems of an ! Augmented Lagrangian framework. When GENCAN is being used ! stand-alone for solving a bound-constrained problem, these ! arguments are dummy arguments and must be ignored. ! ! On Return ! ! f double precision ! objective function value at x ! ! flag integer ! 0 means ''no errors'' ! any other value means ''there was an error in the ! objective function calculation''. ! ! PARAMETERS integer rmax,ncirmax,nobjmax parameter ( rmax = 1000 ) parameter ( ncirmax = 100 ) parameter ( nobjmax = 1000 ) ! COMMON SCALARS integer probid,nrec,r,ncir double precision dx,dy,xmin,xmax,ymin,ymax double precision circle(ncirmax,3) ! SCALAR ARGUMENTS integer flag,n double precision f ! ARRAY ARGUMENTS double precision x(n) ! LOCAL SCALARS integer i,j,k,l double precision dxi,dxj,dyi,dyj,tmpx,tmpy,cx,cy double precision theta,sin_t,cos_t double precision c1i,c2i,pi,c1j,c2j,pj double precision dist,dismin ! LOCAL ARRAYS double precision q(2),w(rmax) double precision cc(nobjmax,ncirmax,2) ! COMMON BLOCKS common /packdata/ dx,dy,xmin,xmax,ymin,ymax,probid,nrec,r common /packcircledata/ circle,ncir flag = 0 f = 0.0d0 theta = x(n) sin_t = sin(theta) cos_t = cos(theta) ! Compute the center of the circles. do i = 1,nrec c1i = x(2*i-1) c2i = x(2*i ) pi = x(2*nrec + i) do k = 1,ncir cx = pi * circle(k,1) + (1.0d0 - pi) * circle(k,2) cy = (1.0d0 - pi) * circle(k,1) + pi * circle(k,2) cc(i,k,1) = cx * cos_t - cy * sin_t + c1i cc(i,k,2) = cx * sin_t + cy * cos_t + c2i end do end do ! Circles overlapping. do i = 1,nrec c1i = x(2*i-1) c2i = x(2*i ) pi = x(2*nrec + i) dxi = (pi*dx + (1.0d0 - pi)*dy) dyi = (pi*dy + (1.0d0 - pi)*dx) do j = i + 1,nrec c1j = x(2*j-1) c2j = x(2*j ) pj = x(2*nrec + j) dxj = (pj*dx + (1.0d0 - pj)*dy) dyj = (pj*dy + (1.0d0 - pj)*dx) tmpx = (dxi + dxj)**2.0d0 - & (c1i*cos_t + c2i*sin_t - c1j*cos_t - c2j*sin_t)**2.0d0 tmpy = (dyi + dyj)**2.0d0 - & (- c1i*sin_t + c2i*cos_t + c1j*sin_t - c2j*cos_t)**2.0d0 ! If the rectangles i and j overlap one another, check the ! circles overlapping. if(tmpx .gt. 0.0d0 .and. tmpy .gt. 0.0d0) then do k = 1,ncir do l = 5,ncir dist = (cc(i,k,1) - cc(j,l,1)) ** 2.0d0 + & (cc(i,k,2) - cc(j,l,2)) ** 2.0d0 dismin = (circle(k,3) + circle(l,3))**2.0d0 f = f + 0.5d0 * max( 0.0d0, dismin - dist ) ** 2 end do end do end if end do ! Region constraints ! SW corner q(1) = - dxi * cos_t + dyi * sin_t + c1i q(2) = - dxi * sin_t - dyi * cos_t + c2i call evalw(probid,q,r,w) do j = 1,r f = f + max( 0.0d0, w(j) ) ** 2 end do ! SE corner q(1) = dxi * cos_t + dyi * sin_t + c1i q(2) = dxi * sin_t - dyi * cos_t + c2i call evalw(probid,q,r,w) do j = 1,r f = f + max( 0.0d0, w(j) ) ** 2 end do ! NE corner q(1) = dxi * cos_t - dyi * sin_t + c1i q(2) = dxi * sin_t + dyi * cos_t + c2i call evalw(probid,q,r,w) do j = 1,r f = f + max( 0.0d0, w(j) ) ** 2 end do ! NW corner q(1) = - dxi * cos_t - dyi * sin_t + c1i q(2) = - dxi * sin_t + dyi * cos_t + c2i call evalw(probid,q,r,w) do j = 1,r f = f + max( 0.0d0, w(j) ) ** 2 end do end do end subroutine evalal_circle ! ****************************************************************** ! ****************************************************************** subroutine evalnal_circle(n,x,g,flag) implicit none ! On Entry: ! n integer ! number of variables ! ! x double precision x(n) ! current point ! ! NOTE: arguments m, lambda and rho are useful when GENCAN is ! being used for solving the box-constrained subproblems of an ! Augmented Lagrangian framework. When GENCAN is being used ! stand-alone for solving a bound-constrained problem, these ! arguments are dummy arguments and must be ignored. ! ! On Return ! ! g double precision g(n) ! gradient of the objective function at x ! ! flag integer ! 0 means ''no errors'', ! any other value means ''there was an error in the ! gradient calculation''. ! ! PARAMETERS integer rmax,ncirmax,nobjmax parameter ( rmax = 1000 ) parameter ( ncirmax = 100 ) parameter ( nobjmax = 1000 ) ! COMMON SCALARS integer probid,nrec,r,ncir double precision dx,dy,xmin,xmax,ymin,ymax,circle(ncirmax,3) ! SCALAR ARGUMENTS integer flag,n ! ARRAY ARGUMENTS double precision g(n),x(n) ! LOCAL SCALARS integer i,j,k,l double precision dxi,dxj,dyi,dyj,tmpx,tmpy double precision gradp,dq1dv,dq2dv,dwdq double precision c1i,c2i,pi,c1j,c2j,pj double precision cos_t,sin_t,theta double precision ck1,ck2,cl1,cl2,dismin,dist,d1dv,d2dv,cy,cx double precision cxi,cyi,cxj,cyj,tmp1,tmp2 ! LOCAL ARRAYS double precision q(2),w(rmax) double precision partial_dwdq1(rmax),partial_dwdq2(rmax) double precision cc(nobjmax,ncirmax,2) ! COMMON BLOCKS common /packdata/ dx,dy,xmin,xmax,ymin,ymax,probid,nrec,r common /packcircledata/ circle,ncir flag = 0 do i = 1,n g(i) = 0.0d0 end do theta = x(n) sin_t = sin(theta) cos_t = cos(theta) ! Compute the center of the circles. do i = 1,nrec c1i = x(2*i-1) c2i = x(2*i ) pi = x(2*nrec + i) do k = 1,ncir cx = pi * circle(k,1) + (1.0d0 - pi) * circle(k,2) cy = (1.0d0 - pi) * circle(k,1) + pi * circle(k,2) cc(i,k,1) = cx * cos_t - cy * sin_t + c1i cc(i,k,2) = cx * sin_t + cy * cos_t + c2i end do end do ! Circles overlapping. do i = 1,nrec c1i = x(2*i-1) c2i = x(2*i ) pi = x(2*nrec + i) dxi = (pi*dx + (1.0d0 - pi)*dy) dyi = (pi*dy + (1.0d0 - pi)*dx) do j = i + 1,nrec c1j = x(2*j-1) c2j = x(2*j ) pj = x(2*nrec + j) dxj = (pj*dx + (1.0d0 - pj)*dy) dyj = (pj*dy + (1.0d0 - pj)*dx) tmpx = (dxi + dxj)**2.0d0 - & (c1i*cos_t + c2i*sin_t - c1j*cos_t - c2j*sin_t)**2.0d0 tmpy = (dyi + dyj)**2.0d0 - & (- c1i*sin_t + c2i*cos_t + c1j*sin_t - c2j*cos_t)**2.0d0 ! If the rectangles i and j overlap one another, check the ! circles overlapping. if(tmpx .gt. 0.0d0 .and. tmpy .gt. 0.0d0) then do k = 1,ncir ck1 = circle(k,1) ck2 = circle(k,2) cxi = cc(i,k,1) cyi = cc(i,k,2) do l = 5,ncir cl1 = circle(l,1) cl2 = circle(l,2) cxj = cc(j,l,1) cyj = cc(j,l,2) dist = (cxi - cxj) ** 2.0d0 + (cyi - cyj) ** 2.0d0 dismin = (circle(k,3) + circle(l,3))**2.0d0 if(dismin .gt. dist) then tmp1 = - 2.0d0 * (dismin - dist) * (cxi - cxj) tmp2 = - 2.0d0 * (dismin - dist) * (cyi - cyj) ! C_i g(2*i-1) = g(2*i-1) + tmp1 g(2*i ) = g(2*i ) + tmp2 ! C_j g(2*j-1) = g(2*j-1) - tmp1 g(2*j ) = g(2*j ) - tmp2 ! p ! p_i d1dv = (ck1 - ck2) * (cos_t + sin_t) d2dv = (ck1 - ck2) * (sin_t - cos_t) gradp = - 2.0d0 * (dismin - dist) * & ((cxi - cxj)*d1dv + (cyi - cyj)*d2dv) g(2*nrec + i) = g(2*nrec + i) + gradp ! p_j d1dv = - (cl1 - cl2) * (cos_t + sin_t) d2dv = - (cl1 - cl2) * (sin_t - cos_t) gradp = - 2.0d0 * (dismin - dist) * & ((cxi - cxj)*d1dv + (cyi - cyj)*d2dv) g(2*nrec + j) = g(2*nrec + j) + gradp ! theta d1dv = - (pi*ck1 + (1.0d0 - pi)*ck2) * sin_t - & ((1.0d0 - pi)*ck1 + pi*ck2) * cos_t - & (- (pj*cl1 + (1.0d0 - pj)*cl2) * sin_t - & ((1.0d0 - pj)*cl1 + pj*cl2) * cos_t) d2dv = (pi*ck1 + (1.0d0 - pi)*ck2) * cos_t - & ((1.0d0 - pi)*ck1 + pi*ck2) * sin_t - & ((pj*cl1 + (1.0d0 - pj)*cl2) * cos_t - & ((1.0d0 - pj)*cl1 + pj*cl2) * sin_t) tmp1 = - 2.0d0 * (dismin - dist) * & ((cxi - cxj)*d1dv + (cyi - cyj)*d2dv) g(n) = g(n) + tmp1 end if end do end do end if end do ! Region constraints ! SW corner q(1) = - dxi * cos_t + dyi * sin_t + c1i q(2) = - dxi * sin_t - dyi * cos_t + c2i call evalw(probid,q,r,w) call evaldw(probid,q,r,partial_dwdq1,partial_dwdq2) ! C do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) g(2*i-1) = g(2*i-1) + 2.0d0 * w(j) * dwdq dwdq = partial_dwdq2(j) g(2*i ) = g(2*i ) + 2.0d0 * w(j) * dwdq end if end do ! p_i dq1dv = - cos_t * (dx - dy) + sin_t * (dy - dx) dq2dv = - sin_t * (dx - dy) - cos_t * (dy - dx) do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(2*nrec + i) = g(2*nrec + i) + 2.0d0 * w(j) * dwdq end if end do ! theta dq1dv = dxi * sin_t + dyi * cos_t dq2dv = - dxi * cos_t + dyi * sin_t do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(n) = g(n) + 2.0d0 * w(j) * dwdq end if end do ! SE corner q(1) = dxi * cos_t + dyi * sin_t + c1i q(2) = dxi * sin_t - dyi * cos_t + c2i call evalw(probid,q,r,w) call evaldw(probid,q,r,partial_dwdq1,partial_dwdq2) ! C do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) g(2*i-1) = g(2*i-1) + 2.0d0 * w(j) * dwdq dwdq = partial_dwdq2(j) g(2*i ) = g(2*i ) + 2.0d0 * w(j) * dwdq end if end do ! p_i dq1dv = cos_t * (dx - dy) + sin_t * (dy - dx) dq2dv = sin_t * (dx - dy) - cos_t * (dy - dx) do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(2*nrec + i) = g(2*nrec + i) + 2.0d0 * w(j) * dwdq end if end do ! theta dq1dv = - dxi * sin_t + dyi * cos_t dq2dv = dxi * cos_t + dyi * sin_t do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(n) = g(n) + 2.0d0 * w(j) * dwdq end if end do ! NE corner q(1) = dxi * cos_t - dyi * sin_t + c1i q(2) = dxi * sin_t + dyi * cos_t + c2i call evalw(probid,q,r,w) call evaldw(probid,q,r,partial_dwdq1,partial_dwdq2) ! C do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) g(2*i-1) = g(2*i-1) + 2.0d0 * w(j) * dwdq dwdq = partial_dwdq2(j) g(2*i ) = g(2*i ) + 2.0d0 * w(j) * dwdq end if end do ! p_i dq1dv = cos_t * (dx - dy) - sin_t * (dy - dx) dq2dv = sin_t * (dx - dy) + cos_t * (dy - dx) do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(2*nrec + i) = g(2*nrec + i) + 2.0d0 * w(j) * dwdq end if end do ! theta dq1dv = - dxi * sin_t - dyi * cos_t dq2dv = dxi * cos_t - dyi * sin_t do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(n) = g(n) + 2.0d0 * w(j) * dwdq end if end do ! NW corner q(1) = - dxi * cos_t - dyi * sin_t + c1i q(2) = - dxi * sin_t + dyi * cos_t + c2i call evalw(probid,q,r,w) call evaldw(probid,q,r,partial_dwdq1,partial_dwdq2) ! C do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) g(2*i-1) = g(2*i-1) + 2.0d0 * w(j) * dwdq dwdq = partial_dwdq2(j) g(2*i ) = g(2*i ) + 2.0d0 * w(j) * dwdq end if end do ! p_i dq1dv = - cos_t * (dx - dy) - sin_t * (dy - dx) dq2dv = - sin_t * (dx - dy) + cos_t * (dy - dx) do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(2*nrec + i) = g(2*nrec + i) + 2.0d0 * w(j) * dwdq end if end do ! theta dq1dv = dxi * sin_t - dyi * cos_t dq2dv = - dxi * cos_t - dyi * sin_t do j = 1,r if ( w(j) .gt. 0.0d0 ) then dwdq = partial_dwdq1(j) * dq1dv + partial_dwdq2(j) * dq2dv g(n) = g(n) + 2.0d0 * w(j) * dwdq end if end do end do end subroutine evalnal_circle ! ****************************************************************** ! ****************************************************************** subroutine defdat(dx,dy,ncir,circle) ! Given the rectangular-item dimensions, this subroutines computes ! some additional data to simplificate further calculations. In ! particular, it defines the sentinel sets, a partial covering of ! a rectangular item by circles and some objects constructed as a ! combination of rectangular items. ! The sentries or sentinels are used in the nonlinear formulation of ! the packing problem as explained in the article. The partial ! covering by circles is used to define an auxiliar problem, easier ! than the original problem, that will be used to find an initial ! approximation of ther solution. Finally, the objects are used ! in the problem definition. ! PARAMETERS integer ncirmax parameter ( ncirmax = 100 ) ! SCALAR ARGUMENTS integer ncir double precision dx,dy ! ARRAY ARGUMENTS double precision circle(ncirmax,3) ! LOCAL SCALARS integer addcir,i double precision tmp ! For each circle, we define its center and its diameter. circle( 1,1) = dx circle( 1,2) = dy circle( 1,3) = 0.0d0 circle( 2,1) = dx circle( 2,2) = - dy circle( 2,3) = 0.0d0 circle( 3,1) = - dx circle( 3,2) = dy circle( 3,3) = 0.0d0 circle( 4,1) = - dx circle( 4,2) = - dy circle( 4,3) = 0.0d0 circle( 5,1) = dx - 0.25d0 * dy circle( 5,2) = dy - 0.25d0 * dy circle( 5,3) = 0.25d0 * dy circle( 6,1) = dx - 0.25d0 * dy circle( 6,2) = - dy + 0.25d0 * dy circle( 6,3) = 0.25d0 * dy circle( 7,1) = - dx + 0.25d0 * dy circle( 7,2) = dy - 0.25d0 * dy circle( 7,3) = 0.25d0 * dy circle( 8,1) = - dx + 0.25d0 * dy circle( 8,2) = - dy + 0.25d0 * dy circle( 8,3) = 0.25d0 * dy circle( 9,1) = dx - 0.50d0 * dy circle( 9,2) = dy - 0.50d0 * dy circle( 9,3) = 0.50d0 * dy circle(10,1) = dx - 0.50d0 * dy circle(10,2) = - dy + 0.50d0 * dy circle(10,3) = 0.50d0 * dy circle(11,1) = - dx + 0.50d0 * dy circle(11,2) = dy - 0.50d0 * dy circle(11,3) = 0.50d0 * dy circle(12,1) = - dx + 0.50d0 * dy circle(12,2) = - dy + 0.50d0 * dy circle(12,3) = 0.50d0 * dy if ( dx .eq. dy ) then ncir = 13 circle(13,1) = 0.0d0 circle(13,2) = 0.0d0 circle(13,3) = dy else addcir = 2 * int( dx / dy ) ncir = 12 + addcir tmp = 2.0d0 * ( dx - dy ) / ( addcir - 1 ) do i = 1,addcir circle(12+i,1) = - dx + dy + ( i - 1 ) * tmp circle(12+i,2) = 0.0d0 circle(12+i,3) = dy end do end if end subroutine defdat ! ****************************************************************** ! ****************************************************************** subroutine drawsol(n,x,solfile) implicit none ! This subroutine draws, for each problem, the convex region and the ! packed items. ! SCALARS ARGUMENTS integer n ! ARRAY ARGUMENTS double precision x(n) character * 11 solfile ! COMMON SCALARS integer probid,nrec,r double precision dx,dy,xmin,xmax,ymin,ymax ! LOCAL SCALARS integer i double precision scale,x1,y1,x2,y2,x3,y3,x4,y4,r_,a,b,c,ai,bi,side,& x5,y5,x6,y6,c1,c2,dxi,dyi,p,theta ! LOCAL ARRAYS double precision q1(2),q2(2),q3(2),q4(2),xr(n) ! COMMON BLOCKS common /packdata/ dx,dy,xmin,xmax,ymin,ymax,probid,nrec,r nrec = ( n - 1 ) / 3 open(unit=10,file=solfile) ! SCALING if ( xmax - xmin .ge. ymax - ymin ) then scale = 13.0d0 / ( xmax - xmin ) else scale = 20.0d0 / ( ymax - ymin ) end if ! CONVEX REGIONS if ( probid .eq. 1 ) then scale = 0.5d0 write(10,10) scale write(10,21) x1 = 0.0d0 y1 = 10.0d0 write(10,22) x1,y1 x2 = 0.0d0 y2 = 3.0d0 write(10,22) x2,y2 x3 = 3.0d0 y3 = 0.0d0 write(10,22) x3,y3 x4 = 10.0d0 y4 = 0.0d0 write(10,25) x4,y4 x5 = 10.0d0 / sqrt(2.0d0) y5 = 10.0d0 / sqrt(2.0d0) write(10,25) x5,y5 write(10,23) else if ( probid .eq. 2 ) then scale = 0.9d0 write(10,10) scale write(10,21) x1 = 6.0d0 y1 = 11.0d0 write(10,22) x1,y1 x2 = 3.0353d0 y2 = 7.5412d0 write(10,25) x2,y2 x3 = 6.0d0 y3 = 5.0d0 write(10,25) x3,y3 x4 = 8.9647d0 y4 = 7.5412d0 write(10,22) x4,y4 write(10,23) else if ( probid .eq. 3 ) then scale = 0.65d0 write(10,10) scale write(10,21) c1 = 6.0d0 c2 = 0.0d0 r_ = 9.0d0 x1 = 0.0d0 y1 = c2 + sqrt( r_ ** 2 - ( c1 - x1 ) ** 2 ) write(10,25) x1,y1 x2 = 4.0d0 y2 = c2 + sqrt( r_ ** 2 - ( c1 - x2 ) ** 2 ) write(10,25) x2,y2 x3 = 8.0d0 y3 = c2 + sqrt( r_ ** 2 - ( c1 - x3 ) ** 2 ) write(10,22) x3,y3 c1 = 1.7d0 c2 = 10.0d0 r_ = 9.0d0 x4 = 8.0d0 y4 = c2 - sqrt( r_ ** 2 - ( c1 - x4 ) ** 2 ) write(10,25) x4,y4 x5 = 4.0d0 y5 = c2 - sqrt( r_ ** 2 - ( c1 - x5 ) ** 2 ) write(10,25) x5,y5 x6 = 0.0d0 y6 = c2 - sqrt( r_ ** 2 - ( c1 - x6 ) ** 2 ) write(10,22) x6,y6 write(10,23) else if ( probid .eq. 4 ) then scale = 1.0d0 write(10,10) scale write(10,20) xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin a = 1.0d0 b = 0.0d0 c = 0.0d0 x1 = -sqrt(5.0d0) do while ( x1 .lt. sqrt(5.0d0) ) y1 = a * x1 ** 2 + b * x1 + c x2 = x1 + 1.0d-2 y2 = a * x2 ** 2 + b * x2 + c write(10,110) x1,y1,x2,y2 x1 = x1 + 1.0d-2 end do a = - 0.25d0 b = 0.00d0 c = 5.00d0 x1 = - 3.0d0 do while ( x1 .lt. 3.0d0 ) y1 = a * x1 ** 2 + b * x1 + c x2 = x1 + 1.0d-2 y2 = a * x2 ** 2 + b * x2 + c write(10,110) x1,y1,x2,y2 x1 = x1 + 1.0d-2 end do else if ( probid .eq. 5 ) then scale = 1.0d0 write(10,10) scale write(10,20) xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin x1 = -sqrt(6.0d0) do while ( x1 .lt. sqrt(6.0d0) ) y1 = x1 ** 2 x2 = x1 + 1.0d-2 y2 = x2 ** 2 write(10,110) x1,y1,x2,y2 x1 = x1 + 1.0d-2 end do y1 = 0.0d0 do while ( y1 .lt. 6.0d0 ) x1 = y1 ** 2 - 6.0d0 * y1 + 6.0d0 y2 = y1 + 1.0d-2 x2 = y2 ** 2 - 6.0d0 * y2 + 6.0d0 if ( x1 .le. 3.0d0 .and. x2 .le. 3.0d0 ) then write(10,110) x1,y1,x2,y2 end if y1 = y1 + 1.0d-2 end do x1 = 0.0d0 y1 = 6.0d0 x2 = 3.0d0 y2 = 3.0d0 write(10,70) x1,y1,x2,y2 else if ( probid .eq. 6 ) then scale = 1.0d0 write(10,10) scale write(10,20) xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin y1 = 0.0d0 do while ( y1 .lt. 6.0d0 ) x1 = y1 ** 2 - 6.0d0 * y1 + 6.0d0 y2 = y1 + 1.0d-2 x2 = y2 ** 2 - 6.0d0 * y2 + 6.0d0 if ( x1 .le. 3.0d0 .and. x2 .le. 3.0d0 ) then write(10,110) x1,y1,x2,y2 end if y1 = y1 + 1.0d-2 end do y1 = -1.5d0 do while ( y1 .lt. 4.5d0 ) x1 = - y1 ** 2 + 3.0d0 * y1 + 0.75d0 y2 = y1 + 1.0d-2 x2 = - y2 ** 2 + 3.0d0 * y2 + 0.75d0 if ( x1 .ge. -3.0d0 .and. x2 .ge. -3.0d0 ) then write(10,110) x1,y1,x2,y2 end if y1 = y1 + 1.0d-2 end do else if ( probid .eq. 7 ) then scale = 0.8d0 write(10,10) scale a = 2.0d0 b = 4.0d0 x1 = 2.0d0 y1 = 4.0d0 write(10,50) 2.0d0*a,2.0d0*b,x1,y1 else if ( probid .eq. 8 ) then scale = 0.6d0 write(10,10) scale a = 2.0d0 b = 6.0d0 x1 = 6.0d0 y1 = 6.0d0 write(10,50) 2.0d0*a,2.0d0*b,x1,y1 a = 6.0d0 b = 2.0d0 x1 = 6.0d0 y1 = 6.0d0 write(10,50) 2.0d0*a,2.0d0*b,x1,y1 x1 = 4.0d0 y1 = 6.0d0 x2 = 6.0d0 y2 = 8.0d0 write(10,70) x1,y1,x2,y2 x1 = 7.0d0 y1 = 4.0d0 x2 = 8.0d0 y2 = 5.0d0 write(10,70) x1,y1,x2,y2 else if ( probid .eq. 9 ) then scale = 0.8d0 write(10,10) scale a = 2.0d0 b = 4.0d0 x1 = 3.0d0 y1 = 4.0d0 write(10,50) 2.0d0*a,2.0d0*b,x1,y1 a = 2.0d0 b = 4.0d0 x1 = 2.65d0 y1 = 4.0d0 write(10,50) 2.0d0*a,2.0d0*b,x1,y1 write(10,21) x1 = 1.0d0 y1 = 0.0d0 write(10,22) x1,y1 x2 = 5.0d0 y2 = 4.0d0 write(10,22) x2,y2 x3 = 1.0d0 y3 = 8.0d0 write(10,27) x3,y3 else if ( probid .eq. 10 ) then scale = 1.0d0 write(10,10) scale a = 6.0d0 b = 2.0d0 x1 = 6.0d0 y1 = 6.0d0 write(10,50) 2.0d0*a,2.0d0*b,x1,y1 a = 3.0d0 b = 3.0d0 x1 = 6.0d0 y1 = 8.0d0 write(10,50) 2.0d0*a,2.0d0*b,x1,y1 else if ( probid .eq. 11 ) then scale = 0.9d0 write(10,10) scale write(10,20) xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin c = 4.0d0 a = 6.0d0 b = 2.0d0 x1 = 0.0d0 y1 = 0.0d0 ai = - a do while ( ai .lt. a ) bi = - b do while ( bi .lt. b ) if ( abs( ( ai / a )**c + ( bi / b )**c - 1.0d0 ) .le. & 1.0d-3 ) then write(10,100) x1+ai,y1+bi end if bi = bi + 3.0d-3 end do ai = ai + 3.0d-3 end do x1 = 0.5d0 y1 = - 2.0d0 x2 = 6.0d0 y2 = 2.0d0 write(10,70) x1,y1,x2,y2 else if ( probid .ge. 12 .and. probid .le. 21 ) then scale = 0.75d0 write(10,10) scale if ( probid .eq. 12 .or. probid .eq. 17 ) then side = 8.619d0 else if ( probid .eq. 13 .or. probid .eq. 18 ) then side = 8.774d0 else if ( probid .eq. 14 .or. probid .eq. 19 ) then side = 9.155d0 else if ( probid .eq. 15 .or. probid .eq. 20 ) then side = 9.302d0 else if ( probid .eq. 16 .or. probid .eq. 21 ) then side = 9.310d0 end if x1 = 0.0d0 y1 = 0.0d0 x2 = 0.5d0 * side y2 = 0.5d0 * side * sqrt( 3.0d0 ) x3 = side y3 = 0.0d0 write(10,80) x1,y1,x2,y2,x3,y3 end if ! OBJECTS theta = x(n) do i = 1,nrec p = x(2*nrec + i) dxi = (p*dx + (1.0d0 - p)*dy) dyi = (p*dy + (1.0d0 - p)*dx) call rotate2(x(2*i-1),x(2*i),xr(2*i-1),xr(2*i),-theta) q1(1) = xr(2*i-1) - dxi q1(2) = xr(2*i ) - dyi q2(1) = xr(2*i-1) + dxi q2(2) = xr(2*i ) - dyi q3(1) = xr(2*i-1) + dxi q3(2) = xr(2*i ) + dyi q4(1) = xr(2*i-1) - dxi q4(2) = xr(2*i ) + dyi call rotate(q1(1),q1(2),theta) call rotate(q2(1),q2(2),theta) call rotate(q3(1),q3(2),theta) call rotate(q4(1),q4(2),theta) write(10,141) q1(1),q1(2),q2(1),q2(2),q3(1),q3(2),q4(1),q4(2) write(10,130) q1(1),q1(2),q2(1),q2(2),q3(1),q3(2),q4(1),q4(2) end do write(10,90) close(10) return 10 format('input parameters;'/,'beginfig(1);' /,'u = ',f20.10,' cm;') 20 format('draw (',f20.10,'u,',f20.10,'u)--',& '(',f20.10,'u,',f20.10,'u)--', & '(',f20.10,'u,',f20.10,'u)--', & '(',f20.10,'u,',f20.10,'u)--cycle;') 21 format('draw ') 22 format(' (',f20.10,'u,',f20.10,'u)--') 25 format(' (',f20.10,'u,',f20.10,'u)..') 23 format(' cycle;') 27 format(' (',f20.10,'u,',f20.10,'u);') 50 format('draw fullcircle ', & 'xscaled ',f20.10,'u yscaled ',f20.10,'u ',& 'shifted (',f20.10,'u,',f20.10,'u);') 70 format('draw (',f20.10,'u,',f20.10,'u)--',& '(',f20.10,'u,',f20.10,'u);') 80 format('draw (',f20.10,'u,',f20.10,'u)--',& '(',f20.10,'u,',f20.10,'u)--', & '(',f20.10,'u,',f20.10,'u)--cycle;') 90 format('endfig;',/,'end;') 100 format('draw (',f20.10,'u,',f20.10,'u);') 110 format('draw (',f20.10,'u,',f20.10,'u)..(',f20.10,'u,',& f20.10,'u);') 130 format('draw (',f20.10,'u,',f20.10,'u)--',& '(',f20.10,'u,',f20.10,'u)--', & '(',f20.10,'u,',f20.10,'u)--', & '(',f20.10,'u,',f20.10,'u)--cycle;') 141 format('fill (',f20.10,'u,',f20.10,'u)--',& '(',f20.10,'u,',f20.10,'u)--', & '(',f20.10,'u,',f20.10,'u)--', & '(',f20.10,'u,',f20.10,'u)--cycle ', & 'withcolor boxColor;') end subroutine drawsol ! ****************************************************************** ! ****************************************************************** subroutine rotate(x,y,theta) double precision x,y,theta,xr,yr xr = x * cos(theta) - y * sin(theta) yr = x * sin(theta) + y * cos(theta) x = xr y = yr end subroutine rotate2(x,y,xr,yr,theta) double precision x,y,theta,xr,yr xr = x * cos(theta) - y * sin(theta) yr = x * sin(theta) + y * cos(theta) end ! ****************************************************************** ! ******************************************************************