      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

!     ******************************************************************
!     ******************************************************************
