C     =================================================================
C     Module: NEWPACK subroutines
C     =================================================================

C     This version: November 10th, 2007.
C
C     This program uses ALGENCAN and the strategy described in:
C
C     E. G. Birgin and F. N. C. Sobral, "Minimizing the object dimen-
C     sions in circle and sphere packing problems", Computers & Oper-
C     ations Research, 2006, doi:10.1016/j.cor.2006.11.002.
C
C     (See also: J. M. Martínez and D. P. Ronconi, "Optimizing the
C     Packing of Cylinders into a Rectangular Container: A Nonlinear
C     Approach", European Journal of Operational Research 160, pp.
C     19-33, 2005.)
C
C     to solve the problem of packing a fixed number of (fixed-size)
C     identical circles into a rectangular region, minimizing the area
C     of the region.

      subroutine newpack(nite,itrial,iterad,refsol,tottime,bestsol,x)

      implicit none

C     PRE-DEFINED PARAMETERS
C
C     nitemax := max. num. of items       (user-defined)
C     itermax := max. num. of steps       (user-defined)
C     ncntmax := dim. of contact array    (>= nitemax * (nitemax-1)/2)
C     lwork   := dim. of working array    (>= ncntmax + nitemax)
C     mmax    := max. num. of constraints (>= ncntmax + nitemax + 4)
C     nmax    := max. num. of variables   (>=       2 * nitemax + 2)
C     rcond   := used to compute the effective rank of the Jacobian
C                matrix, defined as the dimension of the largest
C                triangular matrix given by its QR factorization with
C                pivoting whose estimated condition number is smaller
C                than 1 / rcond
C     rsdeps  := largest acceptable value of the infinity norm of the
C                residual associated with the solution found for the
C                system of equations
      integer mmax,nmax,nitemax,itermax,ncntmax,lwork
      double precision rcond,rsdeps

      parameter (nitemax   =       100)
      parameter (itermax   =      1000)
      parameter (ncntmax   =     10000)
      parameter (lwork     =     10000)
      parameter (mmax      =     10000)
      parameter (nmax      =     10000)
      parameter (rcond     =   1.0d-08)
      parameter (rsdeps    =   1.0d-15)

C     INPUT PARAMETERS
      integer nite,itrial
      double precision iterad,refsol

C     OUTPUT PARAMETERS
      real tottime
      double precision bestsol,x(*)

C     LOCAL VARIABLES
      logical defjac,stapnt
      integer i,j,k,m,n,ncnti,ncnto,iter,info,rank
      integer cnti(2,ncntmax),cnto(2,nitemax)
      integer ncnt(nitemax),jpvt(nmax)
      real time,dtime,dum(2)
      double precision mcheps,tau1,tau2,ratio,area
      double precision maxover,maxviol(4),prevx(nmax)
      double precision d(mmax),e(nmax),f(mmax),w(lwork)
      double precision a(nmax,nmax),jac(mmax,nmax)

C     COMMON VARIABLES
      integer mynite,nloose
      double precision myiterad,myx(nmax)

C     COMMON BLOCKS
      common /mypackdata/ myx,myiterad,mynite,nloose

C     EXTERNAL FUNCTIONS
      double precision supnorm

C     EXTERNAL SUBROUTINES
      external dtime

C     COMPUTING MACHINE-DEPENDENT CONSTANTS
C     REQUIRED FOR THE MODCHL SUBROUTINE
      mcheps = 1.0

 10   if (1.0 + mcheps .ne. 1.0) then
         mcheps = mcheps / 2.0
         go to 10
      end if

      mcheps = mcheps * 2.0
      tau1   = mcheps**(1.0/3.0)
      tau2   = mcheps**(1.0/3.0)

C     INITIALIZING PROBLEM DATA
      mynite   = nite
      myiterad = iterad

      n = 2 * mynite + 2

      do i = 1,n
         myx(i) = x(i)
      end do

C     PRINTING USEFUL INFORMATION
      open(10,file='newpack.out',access='append')

      write(*, 1000) mynite
      write(10,1000) mynite

C     COMPUTING CONTACTS AMONG ITENS AND
C     BETWEEN AN ITEM AND THE OBJECT
      call compcnt(n,mynite,myiterad,myx,ncnti,ncnto,cnti,cnto,ncnt)

C     COMPUTING LOOSE ITEMS AND DROPPING
C     THEM FROM THE PROBLEM TEMPORARILY
      call compite(n,mynite,ncnti,ncnto,nloose,cnti,cnto,ncnt,myx)

C     INITIALIZING ADDITIONAL VARIABLES
      defjac = .false.

      m = ncnti + ncnto

      do i = 1,n
         prevx(i) = myx(i)
      end do

C     MAIN LOOP: TRYING TO IMPROVE SOLUTION
      iter = 1
      time = dtime(dum)

 100  continue

C     COMPUTING VALUE OF SYSTEM'S EQUATIONS
      call compf(m,n,ncnti,ncnto,cnti,cnto,myiterad,myx,f)

C     TESTING PRIMARY STOPPING CRITERION
      if (supnorm(f,1,m) .lt. rsdeps) then
         write(*, 1010)
         write(10,1010)
         go to 200
      end if

C     COMPUTING VALUE OF SYSTEM'S JACOBIAN
      call compjac(m,n,mmax,ncnti,ncnto,cnti,cnto,myiterad,myx,jac)

C     COMPUTING RHS OF NEWTON LINEAR SYSTEM
      do i = 1,m
         d(i) = -f(i)
      end do

C     COMPUTING STEP ALONG NEWTON DIRECTION
      call dgelsy(m,n,1,jac,mmax,d,mmax,jpvt,rcond,rank,w,lwork,info)

C     CHECKING RANK DEFICIENCY OF JACOBIAN
      if (rank .lt. min(m,n)) then
         defjac = .true.

C        CHANGING TO MODIFIED CHOLESKY
C        FACTORIZATION ALGORITHM
         call compjac(m,n,mmax,ncnti,ncnto,cnti,cnto,myiterad,myx,jac)

         do j = 1,n
            do i = j,n
               a(i,j) = 0.0
               do k = 1,m
                  a(i,j) = a(i,j) + jac(k,i) * jac(k,j)
               end do
            end do
         end do

         do i = 1,n
            d(i) = 0.0
            do j = 1,m
               d(i) = d(i) - jac(j,i) * f(j)
            end do
        end do

        call modchl(nmax,n,a,w,mcheps,tau1,tau2,jpvt,e)
        call solsys(nmax,n,a,jpvt,d,d,w,e)
      end if

C     COMPUTING NEXT SOLUTION
      stapnt = .true.

      do i = 1,n
         myx(i) = myx(i) + d(i)

         if (myx(i) .ne. prevx(i)) then
            stapnt = .false.
         end if

         prevx(i) = myx(i)
      end do

C     TESTING SECONDARY STOPPING CRITERIA
      if (stapnt) then
         write(*, 1011)
         write(10,1011)
         go to 200
      end if

      if (iter .eq. itermax) then
         write(*, 1012)
         write(10,1012)
         go to 200
      end if

C     ITERATING
      iter = iter + 1
      go to 100

C     LEAVING MAIN LOOP
 200  continue

C     SCALING SOLUTION IN ORDER TO
C     ELIMINATE ITEM OVERLAPPINGS
      call compcon(n,mynite,myiterad,myx,maxover,maxviol)

      ratio = (2.0 * myiterad) / (2.0*myiterad - maxover)

      do i = 1,n
         myx(i) = myx(i) * ratio
      end do

C     SCALING SOLUTION IN ORDER TO
C     ELIMINATE OBJECT VIOLATIONS
      call compcon(n,mynite,myiterad,myx,maxover,maxviol)

      do i = 1,mynite
         myx(2*i-1) = myx(2*i-1) + maxviol(1)
         myx(2*i  ) = myx(2*i  ) + maxviol(3)
      end do

      myx(n  ) = myx(n  ) + (maxviol(1)+maxviol(2))
      myx(n-1) = myx(n-1) + (maxviol(3)+maxviol(4))

C     UPDATING EXECUTION TIME
      time    = dtime(dum)
      tottime = tottime + dum(1)

C     PRINTING SOLUTION STATISTICS
      area = myx(n) * myx(n-1)

      write(*, 1020) iter,tottime,area,bestsol,area-bestsol
      write(10,1020) iter,tottime,area,bestsol,area-bestsol

C     CLOSING OUTPUT FILE
      close(10)

C     REPLACING DROPPED ITEMS
      call compnlp(time)

      n       = 2 * mynite + 2
      tottime = tottime + time

C     SCALING SOLUTION IN ORDER TO
C     ELIMINATE ITEM OVERLAPPINGS
      call compcon(n,mynite,myiterad,myx,maxover,maxviol)

      ratio = (2.0 * myiterad) / (2.0*myiterad - maxover)

      do i = 1,n
         myx(i) = myx(i) * ratio
      end do

C     SCALING SOLUTION IN ORDER TO
C     ELIMINATE OBJECT VIOLATIONS
      call compcon(n,mynite,myiterad,myx,maxover,maxviol)

      do i = 1,mynite
         myx(2*i-1) = myx(2*i-1) + maxviol(1)
         myx(2*i  ) = myx(2*i  ) + maxviol(3)
      end do

      myx(n  ) = myx(n  ) + (maxviol(1)+maxviol(2))
      myx(n-1) = myx(n-1) + (maxviol(3)+maxviol(4))

C     SAVING PROBLEM SOLUTION
      area = myx(n) * myx(n-1)

      if (area .lt. bestsol) then
         bestsol = area

         call compcon(n,mynite,myiterad,myx,maxover,maxviol)
         call drawsol(n,mynite,ncnti,ncnto,cnti,cnto,myiterad,myx)
         call savesol(n,mynite,ncnti,ncnto,nloose,defjac,stapnt,itrial,
     $                iter,myiterad,refsol,bestsol,maxover,
     $                max(maxviol(1),maxviol(2),maxviol(3),
     $                maxviol(4)),tottime,myx)
      end if

C     NON-EXECUTABLE STATEMENTS
 1000 format(/,' Starting NEWPACK',/,
     $       18X,'Number of items               = ',4X,I8)

 1010 format(/,18X,'New satisfactory solution found')
 1011 format(/,18X,'Best attainable solution achieved')
 1012 format(/,18X,'Maximum number of iterations reached')

 1020 format(/,18X,'Number of iterations          = ',4X,I8,
     $       /,18X,'Total time elapsed (s)        = ',F12.2,
     $       /,18X,'Rectangular object area       = ',F12.8,
     $       /,18X,'Previous best solution        = ',F12.8,
     $       /,18X,'Difference (1st - 2nd)        = ',1P,D12.1)

      end

C     ******************************************************************
C     ******************************************************************

      double precision function supnorm(x,first,last)

      implicit none

C     INPUT PARAMETERS
      integer first,last
      double precision x(*)

C     LOCAL VARIABLES
      integer i

      supnorm = 0.0

      do i = first,last
         supnorm = max(abs(x(i)),supnorm)
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine compcnt(n,mynite,myiterad,myx,ncnti,ncnto,cnti,cnto,
     $                   ncnt)

      implicit none

C     PRE-DEFINED PARAMETERS
C
C     cntieps := minimum distance between two items borders in order for
C                them not to be considered in contact with each other
C     cntoeps := minimum distance between an item and the object borders
C                in order for them not to be considered in contact
      double precision cntieps,cntoeps

      parameter (cntieps   =   1.0d-04)
      parameter (cntoeps   =   1.0d-04)

C     INPUT PARAMETERS
      integer n,mynite
      double precision myiterad,myx(n)

C     OUTPUT PARAMETERS
      integer ncnti,ncnto
      integer cnti(2,*),cnto(2,*),ncnt(*)

C     LOCAL VARIABLES
      integer a,b
      double precision cnt

      ncnti = 0
      ncnto = 0
      do a = 1,mynite
         ncnt(a) = 0
      end do

      do a = 1,mynite
         do b = a+1,mynite
            cnt = (2.0 * myiterad) - sqrt((myx(2*a-1) - myx(2*b-1))**2 +
     $                                    (myx(2*a  ) - myx(2*b  ))**2)

            if (cnt + cntieps .gt. 0.0) then
               ncnti   = ncnti   + 1
               ncnt(a) = ncnt(a) + 1
               ncnt(b) = ncnt(b) + 1
               cnti(1,ncnti) = a
               cnti(2,ncnti) = b
            end if
         end do

         cnt = myiterad - myx(2*a-1)

         if (cnt + cntoeps .gt. 0.0) then
            ncnto   = ncnto   + 1
            ncnt(a) = ncnt(a) + 1
            cnto(1,ncnto) = a
            cnto(2,ncnto) = 1
         end if

         cnt = myx(2*a-1) + myiterad - myx(n)

         if (cnt + cntoeps .gt. 0.0) then
            ncnto   = ncnto   + 1
            ncnt(a) = ncnt(a) + 1
            cnto(1,ncnto) = a
            cnto(2,ncnto) = 2
         end if

         cnt = myiterad - myx(2*a)

         if (cnt + cntoeps .gt. 0.0) then
            ncnto   = ncnto   + 1
            ncnt(a) = ncnt(a) + 1
            cnto(1,ncnto) = a
            cnto(2,ncnto) = 3
         end if

         cnt = myx(2*a) + myiterad - myx(n-1)

         if (cnt + cntoeps .gt. 0.0) then
            ncnto   = ncnto   + 1
            ncnt(a) = ncnt(a) + 1
            cnto(1,ncnto) = a
            cnto(2,ncnto) = 4
         end if
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine compite(n,mynite,ncnti,ncnto,nloose,cnti,cnto,ncnt,myx)

      implicit none

C     OUTPUT PARAMETERS
      integer n,mynite,ncnti,ncnto,nloose
      integer cnti(2,ncnti),cnto(2,ncnto),ncnt(mynite)
      double precision myx(n)

C     LOCAL VARIABLES
      integer a,b,i,j
      double precision objsiz(2),tmp(2)

      nloose    = 0
      objsiz(1) = myx(n-1)
      objsiz(2) = myx(n  )

      a = 1
 10   if (a .le. mynite) then
 100     if (ncnt(a) .eq. 1) then
C           LOOSENING ITEMS THAT
C           MAKE JUST ONE CONTACT
            do i = 1,ncnto
               if (a .eq. cnto(1,i)) then
                  ncnto   = ncnto   - 1
                  ncnt(a) = ncnt(a) - 1
                  do j = i,ncnto
                     cnto(1,j) = cnto(1,j+1)
                     cnto(2,j) = cnto(2,j+1)
                  end do
                  go to 200
               end if
            end do

            do i = 1,ncnti
               if (a .eq. cnti(1,i) .or. a .eq. cnti(2,i)) then
                  if (a .eq. cnti(1,i)) then
                     b = cnti(2,i)
                  else
                     b = cnti(1,i)
                  end if

                  ncnti   = ncnti   - 1
                  ncnt(a) = ncnt(a) - 1
                  ncnt(b) = ncnt(b) - 1
                  do j = i,ncnti
                     cnti(1,j) = cnti(1,j+1)
                     cnti(2,j) = cnti(2,j+1)
                  end do

                  if (ncnt(b) .le. 1) then
                     a = min(a,b)
                     go to 100
                  else
                     go to 200
                  end if
               end if
            end do
         end if

 200     if (ncnt(a) .eq. 0) then
C           DROPPING LOOSE ITEMS FROM
C           THE PROBLEM TEMPORARILY
            nloose = nloose + 1
            mynite = mynite - 1

            tmp(1) = myx(2*a-1)
            tmp(2) = myx(2*a  )

            do i = a,mynite
               ncnt(i)    = ncnt(i+1)
               myx(2*i-1) = myx(2*i+1)
               myx(2*i  ) = myx(2*i+2)
            end do

            ncnt( mynite+1) = 0
            myx(2*mynite+1) = objsiz(1)
            myx(2*mynite+2) = objsiz(2)
            myx(2*mynite+3) = tmp(1)
            myx(2*mynite+4) = tmp(2)

            do i = 1,ncnti
               if (cnti(1,i) .gt. a) then
                  cnti(1,i) = cnti(1,i) - 1
               end if
               if (cnti(2,i) .gt. a) then
                  cnti(2,i) = cnti(2,i) - 1
               end if
            end do

            do i = 1,ncnto
               if (cnto(1,i) .gt. a) then
                  cnto(1,i) = cnto(1,i) - 1
               end if
            end do
         end if

         a = a + 1
         go to 10
      end if

      n = 2 * mynite + 2

      end

C     ******************************************************************
C     ******************************************************************

      subroutine compf(m,n,ncnti,ncnto,cnti,cnto,myiterad,myx,f)

      implicit none

C     INPUT PARAMETERS
      integer m,n,ncnti,ncnto
      integer cnti(2,ncnti),cnto(2,ncnto)
      double precision myiterad,myx(n)

C     OUTPUT PARAMETERS
      double precision f(m)

C     LOCAL VARIABLES
      integer a,b,i

      do i = 1,ncnti
         a = cnti(1,i)
         b = cnti(2,i)

         f(i) = (2.0 * myiterad)**2 - ((myx(2*a-1) - myx(2*b-1))**2 +
     $                                 (myx(2*a  ) - myx(2*b  ))**2)
      end do

      do i = ncnti+1,m
         a = cnto(1,i-ncnti)

         if (cnto(2,i-ncnti) .eq. 1) then

            f(i) = myiterad - myx(2*a-1)

         else if (cnto(2,i-ncnti) .eq. 2) then

            f(i) = myx(2*a-1) + myiterad - myx(n)

         else if (cnto(2,i-ncnti) .eq. 3) then

            f(i) = myiterad - myx(2*a)

         else

            f(i) = myx(2*a) + myiterad - myx(n-1)
         end if
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine compjac(m,n,mmax,ncnti,ncnto,cnti,cnto,myiterad,myx,
     $                   jac)

      implicit none

C     INPUT PARAMETERS
      integer m,n,mmax,ncnti,ncnto
      integer cnti(2,ncnti),cnto(2,ncnto)
      double precision myiterad,myx(n)

C     OUTPUT PARAMETERS
      double precision jac(mmax,n)

C     LOCAL VARIABLES
      integer a,b,c,i,j,k

      do j = 1,n
         do i = 1,m
            jac(i,j) = 0.0
         end do
      end do

      do i = 1,ncnti
         a = cnti(1,i)
         b = cnti(2,i)

         jac(i,2*a-1) = -2.0 * myx(2*a-1) + 2.0 * myx(2*b-1)
         jac(i,2*a  ) = -2.0 * myx(2*a  ) + 2.0 * myx(2*b  )
         jac(i,2*b-1) = -jac(i,2*a-1)
         jac(i,2*b  ) = -jac(i,2*a  )
      end do

      do i = ncnti+1,m
         a = cnto(1,i-ncnti)

         if (cnto(2,i-ncnti) .eq. 1) then

            jac(i,2*a-1) = -1.0

         else if (cnto(2,i-ncnti) .eq. 2) then

            jac(i,2*a-1) =  1.0
            jac(i,n    ) = -1.0

         else if (cnto(2,i-ncnti) .eq. 3) then

            jac(i,2*a) = -1.0

         else

            jac(i,2*a) =  1.0
            jac(i,n-1) = -1.0
         end if
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine compcon(n,mynite,myiterad,myx,maxover,maxviol)

      implicit none

C     INPUT PARAMETERS
      integer n,mynite
      double precision myiterad,myx(n)

C     OUTPUT PARAMETERS
      double precision maxover,maxviol(4)

C     LOCAL VARIABLES
      integer a,b,i
      double precision con

      maxover    = -1.0d+99
      maxviol(1) = -1.0d+99
      maxviol(2) = -1.0d+99
      maxviol(3) = -1.0d+99
      maxviol(4) = -1.0d+99

      do a = 1,mynite
         do b = a+1,mynite
            con = (2.0 * myiterad) - sqrt((myx(2*a-1) - myx(2*b-1))**2 +
     $                                    (myx(2*a  ) - myx(2*b  ))**2)
            maxover = max(con,maxover)
         end do

         con = myiterad - myx(2*a-1)
         maxviol(1) = max(con,maxviol(1))

         con = myx(2*a-1) + myiterad - myx(n)
         maxviol(2) = max(con,maxviol(2))

         con = myiterad - myx(2*a)
         maxviol(3) = max(con,maxviol(3))

         con = myx(2*a) + myiterad - myx(n-1)
         maxviol(4) = max(con,maxviol(4))
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine drawsol(n,mynite,ncnti,ncnto,cnti,cnto,myiterad,myx)

      implicit none

C     INPUT PARAMETERS
      integer n,mynite,ncnti,ncnto
      integer cnti(2,ncnti),cnto(2,ncnto)
      double precision myiterad,myx(n)

C     LOCAL VARIABLES
      integer i

      open(10,file='solution.mp')

      write(10,1000) 20.0 / max(myx(n-1),myx(n))

      write(10,1010) 0.0,0.0,
     $               0.0,myx(n-1),
     $               myx(n),myx(n-1),
     $               myx(n),0.0

      do i = 1,mynite
         write(10,1020) 2.0 * myiterad,
     $                  2.0 * myiterad,
     $                  myx(2*i-1),myx(2*i)
         write(10,1030) myx(2*i-1),myx(2*i)
      end do

      do i = 1,ncnti
         write(10,1040) myx(2*cnti(1,i)-1),myx(2*cnti(1,i)),
     $                  myx(2*cnti(2,i)-1),myx(2*cnti(2,i))
      end do

      do i = 1,ncnto
         if (cnto(2,i) .eq. 1) then
            write(10,1040) myx(2*cnto(1,i)-1),myx(2*cnto(1,i)),
     $                     0.0,myx(2*cnto(1,i))
         else if (cnto(2,i) .eq. 2) then
            write(10,1040) myx(2*cnto(1,i)-1),myx(2*cnto(1,i)),
     $                     myx(n),myx(2*cnto(1,i))
         else if (cnto(2,i) .eq. 3) then
            write(10,1040) myx(2*cnto(1,i)-1),myx(2*cnto(1,i)),
     $                     myx(2*cnto(1,i)-1),0.0
         else
            write(10,1040) myx(2*cnto(1,i)-1),myx(2*cnto(1,i)),
     $                     myx(2*cnto(1,i)-1),myx(n-1)
         end if
      end do

      write(10,1050)

      close(10)

C     NON-EXECUTABLE STATEMENTS
 1000 format('beginfig(1);',/,
     $       'u = ',F20.16,' cm;',/)

 1010 format('draw         (',F20.16,'u,',F20.16,'u)--',/,
     $       '             (',F20.16,'u,',F20.16,'u)--',/,
     $       '             (',F20.16,'u,',F20.16,'u)--',/,
     $       '             (',F20.16,'u,',F20.16,'u)--',/,
     $       '     cycle;',/)

 1020 format('draw fullcircle',/,
     $       '     xscaled  ',F20.16,'u',/,
     $       '     yscaled  ',F20.16,'u',/,
     $       '     shifted (',F20.16,'u,',F20.16,'u);')

 1030 format('draw         (',F20.16,'u,',F20.16,'u)',/,
     $       '     withpen pencircle scaled 4;',/)

 1040 format('draw         (',F20.16,'u,',F20.16,'u)--',/,
     $       '             (',F20.16,'u,',F20.16,'u)',/,
     $       '     dashed evenly;',/)

 1050 format('endfig;',/,'end;')

      end

C     ******************************************************************
C     ******************************************************************

      subroutine savesol(n,mynite,ncnti,ncnto,nloose,defjac,stapnt,
     $                   itrial,iter,myiterad,refsol,bestsol,maxover,
     $                   maxviol,tottime,myx)

      implicit none

C     INPUT PARAMETERS
      logical defjac,stapnt
      integer n,mynite,ncnti,ncnto
      integer nloose,itrial,iter
      real tottime
      double precision myiterad,refsol,bestsol
      double precision maxover,maxviol,myx(n)

C     LOCAL VARIABLES
      integer i

      open(10,file='solution.dat')

      write(10,1000) mynite,myiterad
      write(10,1010) myx(n),myx(n-1),bestsol,refsol
      write(10,1020) bestsol-refsol,maxover,maxviol
      write(10,1030) ncnti+ncnto,ncnti,ncnto
      write(10,1040) nloose,defjac,stapnt
      write(10,1050) itrial,iter,tottime

      write(10,1060)
      do i = 1,mynite
         write(10,1070) i,myx(2*i-1),myx(2*i)
      end do

      close(10)

      open(10,file='solution.csv')
      write(10,1080) mynite,myiterad,myx(n),myx(n-1),bestsol,refsol,
     $               bestsol-refsol,maxover,maxviol,ncnti+ncnto,ncnti,
     $               ncnto,nloose,defjac,stapnt,itrial,iter,tottime
      close(10)

C     NON-EXECUTABLE STATEMENTS
 1000 format('Number of items                    = ',12X,I8,/,
     $       'Radius of packed items             = ',F20.16,/)

 1010 format('Length of outside object           = ',F20.16,/,
     $       'Width of outside object            = ',F20.16,/,
     $       'Area of outside object             = ',F20.16,/,
     $       'Best reference solution            = ',F20.16,/)

 1020 format('Difference (1st - 2nd)             = ',1P,D20.1,/,
     $       'Maximum item overlapping           = ',1P,D20.1,/,
     $       'Maximum object violation           = ',1P,D20.1,/)

 1030 format('Number of detected contacts        = ',12X,I8,/,
     $       'Number of contacts among items     = ',12X,I8,/,
     $       'Number of contacts with object     = ',12X,I8,/)

 1040 format('Number of loose items              = ',12X,I8,/,
     $       'Jacobian matrix is rank deficient  = ',19X,L1,/,
     $       'Newton method stopped convergence  = ',19X,L1,/)

 1050 format('Number of GENPACK trials           = ',12X,I8,/,
     $       'Number of NEWPACK iterations       = ',12X,I8,/,
     $       'Total time elapsed since start (s) = ',F20.2,/)

 1060 format('Items location:',/)

 1070 format(I8,2(1X,F20.16))

 1080 format(I8,',',5(F20.16,','),1P,3(D20.1,','),0P,4(I8,','),
     $       2(L1,','),2(I8,','),F20.2)

      end

C     ******************************************************************
C     ******************************************************************

C     =================================================================
C     Module: ALGENCAN subroutines
C     =================================================================

      subroutine compnlp(time)

      implicit none

C     PRE-DEFINED PARAMETERS
      integer mmax,nmax

      parameter (mmax   =   100000)
      parameter (nmax   =   100000)

C     OUTPUT PARAMETERS
      real time

C     LOCAL VARIABLES
      logical checkder
      logical equatn(mmax),linear(mmax)
      integer gtype,hptype,iprint,inform,m,maxoutit,maxtotfc,maxtotit,n,
     $        ncomp,outiter,totcgcnt,totfcnt,totgcnt,totiter,wi1(nmax)
      double precision epsfeas,epsopt,f,nalpsupn,snorm
      double precision l(nmax),lambda(mmax),rho(mmax),u(nmax),wd1(mmax),
     $                 wd2(mmax),wd3(nmax),wd4(mmax),wd5(mmax),
     $                 wd6(nmax),wd7(nmax),wd8(nmax),wd9(nmax),
     $                 wd10(nmax),wd11(nmax),wd12(nmax),wd13(nmax),
     $                 wd14(nmax),wd15(nmax),wd16(nmax),wd17(mmax),
     $                 wd18(nmax),wd19(nmax),x(nmax)
      character * 6 precond

C     EXTERNAL SUBROUTINES
      external easyalgencan

      call myinip(n,x,l,u,m,lambda,rho,equatn,linear)

      call param(gtype,hptype,precond,checkder,epsfeas,epsopt,maxoutit,
     $           maxtotit,maxtotfc,iprint,ncomp)

      call solver(n,x,l,u,m,lambda,rho,equatn,linear,gtype,hptype,
     $            precond,checkder,epsfeas,epsopt,maxoutit,maxtotit,
     $            maxtotfc,iprint,ncomp,f,snorm,nalpsupn,outiter,
     $            totiter,totfcnt,totgcnt,totcgcnt,time,inform,
     $            wi1,wd1,wd2,wd3,wd4,wd5,wd6,wd7,wd8,wd9,wd10,
     $            wd11,wd12,wd13,wd14,wd15,wd16,wd17,wd18,wd19)

      call myendp(n,x,l,u,m,lambda,rho,equatn,linear)

      end

C     ******************************************************************
C     ******************************************************************

      subroutine myinip(n,x,l,u,m,lambda,rho,equatn,linear)

      implicit none

C     PRE-DEFINED PARAMETERS
      integer nmax

      parameter (nmax   =   10000)

C     OUTPUT PARAMETERS
      logical equatn(*),linear(*)
      integer m,n
      double precision l(*),lambda(*),rho(*),u(*),x(*)

C     LOCAL VARIABLES
      integer i,j,k

C     COMMON VARIABLES
      integer mynite,nloose,icon(3,nmax)
      double precision myiterad,myx(nmax)

C     COMMON BLOCKS
      common /mypackdata/ myx,myiterad,mynite,nloose
      common /myinipdata/ icon

C     Number of variables
      n = 3 * nloose + 1

C     Initial point
      do i = 1,2*nloose
         x(i) = myx(2*mynite+2+i)
      end do

      do i = 2*nloose+1,n
         x(i) = myiterad
      end do

C     Lower and upper bounds
      do i = 1,n
         l(i) = -1.0d+20
         u(i) =  1.0d+20
      end do

C     Number of constraints (equalities plus inequalities)
      m = nloose*mynite + nloose*(nloose-1)/2 + 5*nloose

C     Lagrange multipliers approximation
      do i = 1,m
         lambda(i) = 0.0
      end do

C     Initial penalty parameters
      do i = 1,m
         rho(i) = 10.0
      end do

C     For each constraint i, set equatn(i) = .true. if it is an equality
C     constraint of the form c_i(x) = 0, and set equatn(i) = .false. if
C     it is an inequality constraint of the form c_i(x) <= 0.
      do i = 1,m
         equatn(i) = .false.
      end do

C     For each constraint i, set linear(i) = .true. if it is a linear
C     constraint, otherwise set linear(i) = .false.
      do i = 1,m-5*nloose
         linear(i) = .false.
      end do

      do i = m-5*nloose+1,m
         linear(i) = .true.
      end do

C     Constraints especification
      k = 0

      do i = 1,nloose
         do j = 1,mynite
            k = k + 1
            icon(1,k) = 1
            icon(2,k) = i
            icon(3,k) = j
         end do
      end do

      do i = 1,nloose
         do j = i+1,nloose
            k = k + 1
            icon(1,k) = 2
            icon(2,k) = i
            icon(3,k) = j
         end do
      end do

      do i = 1,nloose
         do j = 1,4
            k = k + 1
            icon(1,k) = 3
            icon(2,k) = i
            icon(3,k) = j
         end do
      end do

      do i = 1,nloose
         k = k + 1
         icon(1,k) = 4
         icon(2,k) = i
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine myevalf(n,x,f,flag)

      implicit none

C     INPUT PARAMETERS
      integer n
      double precision x(n)

C     OUTPUT PARAMETERS
      integer flag
      double precision f

      f = -x(n)

      flag = 0

      end

C     ******************************************************************
C     ******************************************************************

      subroutine myevalg(n,x,g,flag)

      implicit none

C     INPUT PARAMETERS
      integer n
      double precision x(n)

C     OUTPUT PARAMETERS
      integer flag
      double precision g(n)

C     LOCAL VARIABLES
      integer i

      do i = 1,n-1
         g(i) = 0.0
      end do

      g(n) = -1.0

      flag = 0

      end

C     ******************************************************************
C     ******************************************************************

      subroutine myevalh(n,x,hlin,hcol,hval,nnzh,flag)

      implicit none

C     INPUT PARAMETERS
      integer n
      double precision x(n)

C     OUTPUT PARAMETERS
      integer flag,nnzh,hcol(*),hlin(*)
      double precision hval(*)

      nnzh = 0

      flag = 0

      end

C     ******************************************************************
C     ******************************************************************

      subroutine myevalc(n,x,ind,c,flag)

      implicit none

C     PRE-DEFINED PARAMETERS
      integer nmax

      parameter (nmax   =   10000)

C     INPUT PARAMETERS
      integer n,ind
      double precision x(n)

C     OUTPUT PARAMETERS
      integer flag
      double precision c

C     LOCAL VARIABLES
      integer a,b
      double precision objlen,objwid

C     COMMON VARIABLES
      integer mynite,nloose,icon(3,nmax)
      double precision myiterad,myx(nmax)

C     COMMON BLOCKS
      common /mypackdata/ myx,myiterad,mynite,nloose
      common /myinipdata/ icon

      objlen = myx(2*mynite+2)
      objwid = myx(2*mynite+1)

      if (icon(1,ind) .eq. 1) then
         a = icon(2,ind)
         b = icon(3,ind)

         c = (x(2*nloose+a) + myiterad)**2 -
     $       ((x(2*a-1) - myx(2*b-1))**2 + (x(2*a) - myx(2*b))**2)

      else if (icon(1,ind) .eq. 2) then
         a = icon(2,ind)
         b = icon(3,ind)

         c = (x(2*nloose+a) + x(2*nloose+b))**2 -
     $       ((x(2*a-1) - x(2*b-1))**2 + (x(2*a) - x(2*b))**2)

      else if (icon(1,ind) .eq. 3) then
         a = icon(2,ind)

         if (icon(3,ind) .eq. 1) then

            c = x(2*nloose+a) - x(2*a-1)

         else if (icon(3,ind) .eq. 2) then

            c = x(2*a-1) + x(2*nloose+a) - objlen

         else if (icon(3,ind) .eq. 3) then

            c = x(2*nloose+a) - x(2*a)

         else

            c = x(2*a) + x(2*nloose+a) - objwid
         end if

      else
         a = icon(2,ind)

         c = x(n) - x(2*nloose+a)
      end if

      flag = 0

      end

C     ******************************************************************
C     ******************************************************************

      subroutine myevaljac(n,x,ind,indjac,valjac,nnzjac,flag)

      implicit none

C     PRE-DEFINED PARAMETERS
      integer nmax

      parameter (nmax   =   10000)

C     INPUT PARAMETERS
      integer n,ind
      double precision x(n)

C     OUTPUT PARAMETERS
      integer flag,nnzjac,indjac(n)
      double precision valjac(n)

C     LOCAL VARIABLES
      integer a,b

C     COMMON VARIABLES
      integer mynite,nloose,icon(3,nmax)
      double precision myiterad,myx(nmax)

C     COMMON BLOCKS
      common /mypackdata/ myx,myiterad,mynite,nloose
      common /myinipdata/ icon

      if (icon(1,ind) .eq. 1) then
         a = icon(2,ind)
         b = icon(3,ind)

         nnzjac = 3

         indjac(1) = 2 * a - 1
         indjac(2) = 2 * a
         indjac(3) = 2 * nloose + a

         valjac(1) = -2.0 * x(2*a-1)      + 2.0 * myx(2*b-1)
         valjac(2) = -2.0 * x(2*a  )      + 2.0 * myx(2*b  )
         valjac(3) =  2.0 * x(2*nloose+a) + 2.0 * myiterad

      else if (icon(1,ind) .eq. 2) then
         a = icon(2,ind)
         b = icon(3,ind)

         nnzjac = 6

         indjac(1) = 2 * a - 1
         indjac(2) = 2 * a
         indjac(3) = 2 * nloose + a
         indjac(4) = 2 * b - 1
         indjac(5) = 2 * b
         indjac(6) = 2 * nloose + b

         valjac(1) = -2.0 * x(2*a-1)      + 2.0 * x(2*b-1)
         valjac(2) = -2.0 * x(2*a  )      + 2.0 * x(2*b  )
         valjac(3) =  2.0 * x(2*nloose+a) + 2.0 * x(2*nloose+b)
         valjac(4) = -valjac(1)
         valjac(5) = -valjac(2)
         valjac(6) =  valjac(3)

      else if (icon(1,ind) .eq. 3) then
         a = icon(2,ind)

         nnzjac = 2

         if (icon(3,ind) .eq. 1) then

            indjac(1) = 2 * a - 1
            indjac(2) = 2 * nloose + a

            valjac(1) = -1.0
            valjac(2) =  1.0

         else if (icon(3,ind) .eq. 2) then

            indjac(1) = 2 * a - 1
            indjac(2) = 2 * nloose + a

            valjac(1) =  1.0
            valjac(2) =  1.0

         else if (icon(3,ind) .eq. 3) then

            indjac(1) = 2 * a
            indjac(2) = 2 * nloose + a

            valjac(1) = -1.0
            valjac(2) =  1.0

         else

            indjac(1) = 2 * a
            indjac(2) = 2 * nloose + a

            valjac(1) =  1.0
            valjac(2) =  1.0
         end if

      else
         a = icon(2,ind)

         nnzjac = 2

         indjac(1) = 2 * nloose + a
         indjac(2) = n

         valjac(1) = -1.0
         valjac(2) =  1.0
      end if

      flag = 0

      end

C     ******************************************************************
C     ******************************************************************

      subroutine myevalhc(n,x,ind,hclin,hccol,hcval,nnzhc,flag)

      implicit none

C     PRE-DEFINED PARAMETERS
      integer nmax

      parameter (nmax   =   10000)

C     INPUT PARAMETERS
      integer n,ind
      double precision x(n)

C     OUTPUT PARAMETERS
      integer flag,nnzhc,hccol(*),hclin(*)
      double precision hcval(*)

C     LOCAL VARIABLES
      integer a,b

C     COMMON VARIABLES
      integer mynite,nloose,icon(3,nmax)
      double precision myiterad,myx(nmax)

C     COMMON BLOCKS
      common /mypackdata/ myx,myiterad,mynite,nloose
      common /myinipdata/ icon

      if (icon(1,ind) .eq. 1) then
         a = icon(2,ind)
         b = icon(3,ind)

         nnzhc = 3

         hclin(1) = 2 * a - 1
         hclin(2) = 2 * a
         hclin(3) = 2 * nloose + a

         hccol(1) = 2 * a - 1
         hccol(2) = 2 * a
         hccol(3) = 2 * nloose + a

         hcval(1) = -2.0
         hcval(2) = -2.0
         hcval(3) =  2.0

      else if (icon(1,ind) .eq. 2) then
         a = icon(2,ind)
         b = icon(3,ind)

         nnzhc = 9

         hclin(1) = 2 * a - 1
         hclin(2) = 2 * a
         hclin(3) = 2 * nloose + a
         hclin(4) = 2 * b - 1
         hclin(5) = 2 * b
         hclin(6) = 2 * nloose + b
         hclin(7) = 2 * b - 1
         hclin(8) = 2 * b
         hclin(9) = 2 * nloose + b

         hccol(1) = 2 * a - 1
         hccol(2) = 2 * a
         hccol(3) = 2 * nloose + a
         hccol(4) = 2 * b - 1
         hccol(5) = 2 * b
         hccol(6) = 2 * nloose + b
         hccol(7) = 2 * a - 1
         hccol(8) = 2 * a
         hccol(9) = 2 * nloose + a

         hcval(1) = -2.0
         hcval(2) = -2.0
         hcval(3) =  2.0
         hcval(4) = -2.0
         hcval(5) = -2.0
         hcval(6) =  2.0
         hcval(7) =  2.0
         hcval(8) =  2.0
         hcval(9) =  2.0

      else
         nnzhc = 0
      end if

      flag = 0

      end

C     ******************************************************************
C     ******************************************************************

      subroutine myendp(n,x,l,u,m,lambda,rho,equatn,linear)

      implicit none

C     PRE-DEFINED PARAMETERS
      integer nmax

      parameter (nmax   =   10000)

C     INPUT PARAMETERS
      logical equatn(*),linear(*)
      integer m,n
      double precision l(*),lambda(*),rho(*),u(*),x(*)

C     LOCAL VARIABLES
      integer i
      double precision objsiz(2)

C     COMMON VARIABLES
      integer mynite,nloose
      double precision myiterad,myx(nmax)

C     COMMON BLOCKS
      common /mypackdata/ myx,myiterad,mynite,nloose

      objsiz(1) = myx(2*mynite+1)
      objsiz(2) = myx(2*mynite+2)

      do i = 1,2*nloose
         myx(2*mynite+i) = x(i)
      end do

      mynite = mynite + nloose
      myx(2*mynite+1) = objsiz(1)
      myx(2*mynite+2) = objsiz(2)

      end
