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

      subroutine inip(n,x,l,u,m,lambda,rho,equatn,linear,
     +ml,meq,alin,acol,aval,annz,rhs,fknown)

      implicit none

C     SCALAR ARGUMENTS
      integer annz,m,meq,ml,n
      double precision fknown

C     ARRAY ARGUMENTS
      logical equatn(*),linear(*)
      integer alin(*),acol(*)
      double precision aval(*),l(*),lambda(*),rho(*),rhs(*),u(*),x(*)

C     LOCAL SCALARS
      integer i,j,k
      double precision seed

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer Q
      double precision drand
      external drand,Q

C     Problem data

      write(*,*) 'Enter a, b, and nballs: '
      read(*,*) a,b,nballs

c     a = 2.0d0
c     b = 1.0d0
      r = 1.0d0
c     nballs = 2

C     Variables

      n = 3 * nballs

      do i = 1,nballs
          l(Q('u',i)) = - a
          u(Q('u',i)) =   a
          l(Q('v',i)) = - b
          u(Q('v',i)) =   b
          l(Q('s',i)) = 0.0d0
          u(Q('s',i)) = 1.0d0
      end do

c     write(*,*) 'Enter seed: '
c     read(*,*) seed
      seed = 123456.0d0

      do i = 1,n
          x(i) = l(i) + ( u(i) - l(i) ) * drand(seed)
      end do

C     Upper-level Constraints

      m = 2 * nballs + nballs * ( nballs - 1 ) / 2

      k = 0
      do i = 1,nballs
          do j = i + 1,nballs
              k = k + 1
              pairs(1,k) = i
              pairs(2,k) = j
          end do
      end do

      do i = 1,m
          linear(i) = .false.
          lambda(i) =  0.0d0
          rho(i)    = 10.0d0
      end do

      do i = 1,nballs
          equatn(i) = .true.
      end do

      do i = nballs + 1,2 * nballs
          equatn(i) = .false.
      end do

      do i = 2 * nballs + 1,2 * nballs + nballs * ( nballs - 1 ) / 2
          equatn(i) = .false.
      end do

C     Lower-level linear constraints

      ml   = 0
      meq  = 0
      annz = 0

      fknown = 1.0d+20

c     go to 123
c     if ( a .eq. 4.0d0 .and. b .eq. 2.0d0 ) then
c         if ( nballs .eq. 2 ) then
c             fknown = 3.6000632560807304D+01
c         else if ( nballs .eq. 3 ) then
c             fknown = 5.8521317671731637D+01
c         else if ( nballs .eq. 4 ) then
c             fknown = 9.0072124538555983D+01
c         else if ( nballs .eq. 5 ) then
c             fknown = 1.1780270817412281D+02
c         end if
c     else if ( a .eq. 3.0d0 .and. b .eq. 2.0d0 ) then
c         if ( nballs .eq. 2 ) then
c             fknown = 1.6000808541953081D+01
c         else if ( nballs .eq. 3 ) then
c             fknown = 2.6239622360900952D+01
c         else if ( nballs .eq. 4 ) then
c             fknown = 4.0410065156787176D+01
c         end if
c     else if ( a .eq. 2.0d0 .and. b .eq. 2.0d0 ) then
c         if ( nballs .eq. 2 ) then
c             fknown = 3.9999950333831871D+00
c         end if
c     end if
c123  continue

      end

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

      subroutine evalf(n,x,f,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n
      double precision f

C     ARRAY ARGUMENTS
      double precision x(n)

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     LOCAL SCALARS
      integer i,j
      double precision cte,dist,ui,uj,vi,vj,si,sj,xxi,xxj,yyi,yyj

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer Q
      external Q

      flag = 0

      f = 0.0d0

      cte = b ** 2 / a ** 2

      do i = 1,nballs
          do j = i + 1,nballs
              ui = x(Q('u',i))
              vi = x(Q('v',i))
              si = x(Q('s',i))

              xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
              yyi = si * vi

              uj = x(Q('u',j))
              vj = x(Q('v',j))
              sj = x(Q('s',j))

              xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
              yyj = sj * vj

              dist = ( xxi - xxj ) ** 2 + ( yyi - yyj ) ** 2
              f = f - dist
          end do
      end do              

      end

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

      subroutine evalg(n,x,g,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n

C     ARRAY ARGUMENTS
      double precision g(n),x(n)

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     LOCAL SCALARS
      integer i,j
      double precision cte,diffx,diffy,ui,uj,vi,vj,si,sj,xxi,xxj,yyi,yyj

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer Q
      external Q

      flag = 0

      do i = 1,nballs
          g(Q('u',i)) = 0.0d0
          g(Q('v',i)) = 0.0d0
          g(Q('s',i)) = 0.0d0
      end do              

      cte = b ** 2 / a ** 2

      do i = 1,nballs
          do j = i + 1,nballs
              ui = x(Q('u',i))
              vi = x(Q('v',i))
              si = x(Q('s',i))

              xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
              yyi = si * vi

              uj = x(Q('u',j))
              vj = x(Q('v',j))
              sj = x(Q('s',j))

              xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
              yyj = sj * vj

              diffx = 2.0d0 * ( xxi - xxj )
              diffy = 2.0d0 * ( yyi - yyj )

              g(Q('u',i)) = 
     +        g(Q('u',i)) - diffx * ( 1.0d0 + ( si - 1.0d0 ) * cte )
              g(Q('u',j)) = 
     +        g(Q('u',j)) + diffx * ( 1.0d0 + ( sj - 1.0d0 ) * cte )

              g(Q('v',i)) = g(Q('v',i)) - diffy * si
              g(Q('v',j)) = g(Q('v',j)) + diffy * sj

              g(Q('s',i)) = g(Q('s',i)) - diffx * cte * ui - diffy * vi
              g(Q('s',j)) = g(Q('s',j)) + diffx * cte * uj + diffy * vj
          end do
      end do              

      end

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

      subroutine evalh(n,x,hlin,hcol,hval,hnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n,hnnz

C     ARRAY ARGUMENTS
      integer hcol(*),hlin(*)
      double precision hval(*),x(n)

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     LOCAL SCALARS
      integer i,j
      double precision cte,diffx,diffy,dxxidui,dxxjduj,ui,uj,vi,vj,si,
     +        sj,xxi,xxj,yyi,yyj

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer Q
      external Q

      flag = 0

      cte = b ** 2 / a ** 2

      hnnz = 0

      do i = 1,nballs
          do j = i + 1,nballs
              ui = x(Q('u',i))
              vi = x(Q('v',i))
              si = x(Q('s',i))

              xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
              yyi = si * vi

              uj = x(Q('u',j))
              vj = x(Q('v',j))
              sj = x(Q('s',j))

              xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
              yyj = sj * vj

              diffx = 2.0d0 * ( xxi - xxj )
              diffy = 2.0d0 * ( yyi - yyj )

              dxxidui = ( 1.0d0 + ( si - 1.0d0 ) * cte )
              dxxjduj = ( 1.0d0 + ( sj - 1.0d0 ) * cte )

              hlin(hnnz+1) = Q('u',i)
              hcol(hnnz+1) = Q('u',i)
              hval(hnnz+1) = - 2.0d0 * dxxidui ** 2

              hlin(hnnz+2) = Q('s',i)
              hcol(hnnz+2) = Q('u',i)
              hval(hnnz+2) = - 2.0d0 * cte * ui * dxxidui - diffx * cte

              hlin(hnnz+3) = Q('u',j)
              hcol(hnnz+3) = Q('u',i)
              hval(hnnz+3) = 2.0d0 * dxxjduj * dxxidui

              hlin(hnnz+4) = Q('s',j)
              hcol(hnnz+4) = Q('u',i)
              hval(hnnz+4) = 2.0d0 * cte * uj * dxxidui 

              hlin(hnnz+5) = Q('v',i)
              hcol(hnnz+5) = Q('v',i)
              hval(hnnz+5) = - 2.0d0 * si ** 2

              hlin(hnnz+6) = Q('s',i)
              hcol(hnnz+6) = Q('v',i)
              hval(hnnz+6) = - 4.0d0 * yyi +  2.0d0 * yyj

              hlin(hnnz+7) = Q('v',j)
              hcol(hnnz+7) = Q('v',i)
              hval(hnnz+7) = 2.0d0 * sj * si

              hlin(hnnz+8) = Q('s',j)
              hcol(hnnz+8) = Q('v',i)
              hval(hnnz+8) = 2.0d0 * vj * si

              hlin(hnnz+9) = Q('s',i)
              hcol(hnnz+9) = Q('s',i)
              hval(hnnz+9) = - 2.0d0 * ( cte ** 2 * ui ** 2 + vi ** 2 )

              hlin(hnnz+10) = Q('s',j)
              hcol(hnnz+10) = Q('s',i)
              hval(hnnz+10) = 2.0d0 * ( cte ** 2 * ui * uj + vj * vi )

              hlin(hnnz+11) = Q('s',i)
              hcol(hnnz+11) = Q('u',j)
              hval(hnnz+11) = 2.0d0 * cte * ui * dxxjduj

              hlin(hnnz+12) = Q('u',j)
              hcol(hnnz+12) = Q('u',j)
              hval(hnnz+12) = - 2.0d0 * dxxjduj ** 2 

              hlin(hnnz+13) = Q('s',j)
              hcol(hnnz+13) = Q('u',j)
              hval(hnnz+13) = - 2.0d0 * cte * uj * dxxjduj + diffx * cte

              hlin(hnnz+14) = Q('s',i)
              hcol(hnnz+14) = Q('v',j)
              hval(hnnz+14) = 2.0d0 * vi * sj

              hlin(hnnz+15) = Q('v',j)
              hcol(hnnz+15) = Q('v',j)
              hval(hnnz+15) = - 2.0d0 * sj ** 2

              hlin(hnnz+16) = Q('s',j)
              hcol(hnnz+16) = Q('v',j)
              hval(hnnz+16) = 2.0d0 * yyi - 4.0d0 * yyj

              hlin(hnnz+17) = Q('s',j)
              hcol(hnnz+17) = Q('s',j)
              hval(hnnz+17) = - 2.0d0 * ( cte ** 2 * uj ** 2 + vj ** 2 )

              hnnz = hnnz + 17
          end do
      end do              

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalc(n,x,ind,c,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer ind,flag,n
      double precision c

C     ARRAY ARGUMENTS
      double precision x(n)

C     LOCAL SCALARS
      integer i,j,k
      double precision cte,ui,uj,vi,vj,si,sj,xxi,xxj,yyi,yyj

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer Q
      external Q

      flag = 0

      if ( ind .le. nballs ) then
          i  = ind

          ui = x(Q('u',i))
          vi = x(Q('v',i))

          c = ( ui / a ) ** 2 + ( vi / b ) ** 2 - 1.0d0

      else if ( ind .le. 2 * nballs ) then
          i  = ind - nballs

          cte = ( b ** 2 / a ** 2 ) ** 2

          ui = x(Q('u',i))
          vi = x(Q('v',i))
          si = x(Q('s',i))

          c = r ** 2 - ( si - 1.0d0 ) ** 2 * ( cte * ui ** 2 + vi ** 2 )

      else
          k = ind - 2 * nballs
          i = pairs(1,k)
          j = pairs(2,k)

          cte = b ** 2 / a ** 2

          ui = x(Q('u',i))
          vi = x(Q('v',i))
          si = x(Q('s',i))

          xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
          yyi = si * vi

          uj = x(Q('u',j))
          vj = x(Q('v',j))
          sj = x(Q('s',j))

          xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
          yyj = sj * vj

          c = 4.0d0 * r ** 2 - ( xxi - xxj ) ** 2 - ( yyi - yyj ) ** 2
      end if

      end

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

      subroutine evaljac(n,x,ind,jcvar,jcval,jcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,jcnnz,n

C     ARRAY ARGUMENTS
      integer jcvar(n)
      double precision x(n),jcval(n)

C     LOCAL SCALARS
      integer i,j,k
      double precision cte,diffx,diffy,ui,uj,vi,vj,si,sj,xxi,xxj,yyi,yyj

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer Q
      external Q

      flag = 0

      if ( ind .le. nballs ) then
          i  = ind

          ui = x(Q('u',i))
          vi = x(Q('v',i))

          jcnnz = 2

          jcvar(1) = Q('u',i)
          jcval(1) = 2.0d0 * ui / a ** 2

          jcvar(2) = Q('v',i)
          jcval(2) = 2.0d0 * vi / b ** 2

      else if ( ind .le. 2 * nballs ) then
          i  = ind - nballs

          cte = ( b ** 2 / a ** 2 ) ** 2

          ui = x(Q('u',i))
          vi = x(Q('v',i))
          si = x(Q('s',i))

          jcnnz = 3

          jcvar(1) = Q('u',i)
          jcval(1) = - ( si - 1.0d0 ) ** 2 * 2.0d0 * cte * ui

          jcvar(2) = Q('v',i)
          jcval(2) = - ( si - 1.0d0 ) ** 2 * 2.0d0 * vi

          jcvar(3) = Q('s',i)
          jcval(3) = - 2.0d0 * ( si - 1.0d0 ) * 
     +                 ( cte * ui ** 2 + vi ** 2 )

      else
          k = ind - 2 * nballs
          i = pairs(1,k)
          j = pairs(2,k)

          cte = b ** 2 / a ** 2

          ui = x(Q('u',i))
          vi = x(Q('v',i))
          si = x(Q('s',i))

          xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
          yyi = si * vi

          uj = x(Q('u',j))
          vj = x(Q('v',j))
          sj = x(Q('s',j))

          xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
          yyj = sj * vj

          diffx = 2.0d0 * ( xxi - xxj )
          diffy = 2.0d0 * ( yyi - yyj )

          jcnnz = 6

          jcvar(1) = Q('u',i)
          jcval(1) = - diffx * ( 1.0d0 + ( si - 1.0d0 ) * cte )

          jcvar(2) = Q('v',i)
          jcval(2) = - diffy * si

          jcvar(3) = Q('s',i)
          jcval(3) = - diffx * cte * ui - diffy * vi

          jcvar(4) = Q('u',j)
          jcval(4) = diffx * ( 1.0d0 + ( sj - 1.0d0 ) * cte )

          jcvar(5) = Q('v',j)
          jcval(5) = diffy * sj

          jcvar(6) = Q('s',j)
          jcval(6) = diffx * cte * uj + diffy * vj
      end if

      end

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

      subroutine evalhc(n,x,ind,hclin,hccol,hcval,hcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,hcnnz,ind,n

C     ARRAY ARGUMENTS
      integer hccol(*),hclin(*)
      double precision hcval(*),x(n)

C     LOCAL SCALARS
      integer i,j,k
      double precision cte,diffx,diffy,dxxidui,dxxjduj,ui,uj,vi,vj,si,
     +        sj,xxi,xxj,yyi,yyj

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer Q
      external Q

      flag = 0

      if ( ind .le. nballs ) then
          i  = ind

          ui = x(Q('u',i))
          vi = x(Q('v',i))

          hcnnz = 2

          hclin(1) = Q('u',i)
          hccol(1) = Q('u',i)
          hcval(1) = 2.0d0 / a ** 2

          hclin(2) = Q('v',i)
          hccol(2) = Q('v',i)
          hcval(2) = 2.0d0 / b ** 2

      else if ( ind .le. 2 * nballs ) then
          i  = ind - nballs

          cte = ( b ** 2 / a ** 2 ) ** 2

          ui = x(Q('u',i))
          vi = x(Q('v',i))
          si = x(Q('s',i))

          hcnnz = 5

          hclin(1) = Q('u',i)
          hccol(1) = Q('u',i)
          hcval(1) = - ( si - 1.0d0 ) ** 2 * 2.0d0 * cte

          hclin(2) = Q('s',i)
          hccol(2) = Q('u',i)
          hcval(2) = - 4.0d0 * ( si - 1.0d0 ) * cte * ui

          hclin(3) = Q('v',i)
          hccol(3) = Q('v',i)
          hcval(3) = - ( si - 1.0d0 ) ** 2 * 2.0d0

          hclin(4) = Q('s',i)
          hccol(4) = Q('v',i)
          hcval(4) = - 4.0d0 * ( si - 1.0d0 ) * vi

          hclin(5) = Q('s',i)
          hccol(5) = Q('s',i)
          hcval(5) = - 2.0d0 * ( cte * ui ** 2 + vi ** 2 )

      else
          k = ind - 2 * nballs
          i = pairs(1,k)
          j = pairs(2,k)

          cte = b ** 2 / a ** 2

          ui = x(Q('u',i))
          vi = x(Q('v',i))
          si = x(Q('s',i))

          xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
          yyi = si * vi

          uj = x(Q('u',j))
          vj = x(Q('v',j))
          sj = x(Q('s',j))

          xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
          yyj = sj * vj

          diffx = 2.0d0 * ( xxi - xxj )
          diffy = 2.0d0 * ( yyi - yyj )

          dxxidui = ( 1.0d0 + ( si - 1.0d0 ) * cte )
          dxxjduj = ( 1.0d0 + ( sj - 1.0d0 ) * cte )

          hcnnz = 17

          hclin(1) = Q('u',i)
          hccol(1) = Q('u',i)
          hcval(1) = - 2.0d0 * dxxidui ** 2

          hclin(2) = Q('s',i)
          hccol(2) = Q('u',i)
          hcval(2) = - 2.0d0 * cte * ui * dxxidui - diffx * cte

          hclin(3) = Q('u',j)
          hccol(3) = Q('u',i)
          hcval(3) = 2.0d0 * dxxjduj * dxxidui

          hclin(4) = Q('s',j)
          hccol(4) = Q('u',i)
          hcval(4) = 2.0d0 * cte * uj * dxxidui 

          hclin(5) = Q('v',i)
          hccol(5) = Q('v',i)
          hcval(5) = - 2.0d0 * si ** 2

          hclin(6) = Q('s',i)
          hccol(6) = Q('v',i)
          hcval(6) = - 4.0d0 * yyi +  2.0d0 * yyj

          hclin(7) = Q('v',j)
          hccol(7) = Q('v',i)
          hcval(7) = 2.0d0 * sj * si

          hclin(8) = Q('s',j)
          hccol(8) = Q('v',i)
          hcval(8) = 2.0d0 * vj * si

          hclin(9) = Q('s',i)
          hccol(9) = Q('s',i)
          hcval(9) = - 2.0d0 * ( cte ** 2 * ui ** 2 + vi ** 2 )

          hclin(10) = Q('s',j)
          hccol(10) = Q('s',i)
          hcval(10) = 2.0d0 * ( cte ** 2 * ui * uj + vj * vi )

          hclin(11) = Q('s',i)
          hccol(11) = Q('u',j)
          hcval(11) = 2.0d0 * cte * ui * dxxjduj

          hclin(12) = Q('u',j)
          hccol(12) = Q('u',j)
          hcval(12) = - 2.0d0 * dxxjduj ** 2 

          hclin(13) = Q('s',j)
          hccol(13) = Q('u',j)
          hcval(13) = - 2.0d0 * cte * uj * dxxjduj + diffx * cte

          hclin(14) = Q('s',i)
          hccol(14) = Q('v',j)
          hcval(14) = 2.0d0 * vi * sj

          hclin(15) = Q('v',j)
          hccol(15) = Q('v',j)
          hcval(15) = - 2.0d0 * sj ** 2

          hclin(16) = Q('s',j)
          hccol(16) = Q('v',j)
          hcval(16) = 2.0d0 * yyi - 4.0d0 * yyj

          hclin(17) = Q('s',j)
          hccol(17) = Q('s',j)
          hcval(17) = - 2.0d0 * ( cte ** 2 * uj ** 2 + vj ** 2 )
      end if

      end

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

      subroutine evalhlp(n,x,m,lambda,p,hp,goth,flag)

      implicit none

C     SCALAR ARGUMENTS
      logical goth
      integer flag,m,n

C     ARRAY ARGUMENTS
      double precision hp(n),lambda(m),p(n),x(n)

      flag = - 1

      end

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

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

      implicit none

C     SCALAR ARGUMENTS
      integer m,n

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      double precision l(n),lambda(m),rho(m),u(n),x(n)

      end

C     ******************************************************************
C     ******************************************************************
 
      integer function Q(vtype,vind)

      implicit none

C     SCALAR ARGUMENTS
      integer vind
      character vtype

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

      if ( vtype .eq. 'u' ) then
          Q = vind
      else if ( vtype .eq. 'v' ) then
          Q = vind + nballs
      else if ( vtype .eq. 's' ) then
          Q = vind + 2 * nballs
      else
          write(*,*) 'Error in Q function!'
          stop
      end if

      return

      end

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

      subroutine addc(n,l,u,f,ml,meq,alin,acol,aval,annz,b,fknown,mla,
     +annza,ind)

C     SCALAR ARGUMENTS
      integer annz,annza,meq,ml,mla,n
      double precision f,fknown

C     ARRAY ARGUMENTS
      integer alin(annz),acol(annz),ind(ml)
      double precision aval(annz),b(ml),l(n),u(n)

      mla   = ml
      annza = annz

      end

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

      subroutine setlv(n,linvar)

      implicit none

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      logical linvar(n)

C     LOCAL SCALARS
      integer i

      do i = 1,n
          linvar(i) = .false.
      end do

      end

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

      subroutine evalfi(n,y,f,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n
      double precision f(2)

C     ARRAY ARGUMENTS
      double precision y(2*n)

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     LOCAL SCALARS
      integer i,j
      double precision cte

C     LOCAL ARRAYS
      double precision diffx(2),diffy(2),dist(2),ui(2),uj(2),vi(2),
     +        vj(2),si(2),sj(2),xxi(2),xxj(2),yyi(2),yyj(2),tmp(2)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer P,Q
      external P,Q

      flag = 0

c     f = 0.0d0
      call ivl1(0.0d0,f)

      cte = b ** 2 / a ** 2

      do i = 1,nballs
          do j = i + 1,nballs
c             ui = x(Q('u',i))
              call ivli(y(P(Q('u',i))),ui)
c             vi = x(Q('v',i))
              call ivli(y(P(Q('v',i))),vi)
c             si = x(Q('s',i))
              call ivli(y(P(Q('s',i))),si)

c             xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
              call scladd(-1.0d0,si,tmp)
              call sclmlt(cte,tmp,tmp)
              call scladd(1.0d0,tmp,tmp)
              call mult(tmp,ui,xxi)

c             yyi = si * vi
              call mult(si,vi,yyi)

c             uj = x(Q('u',j))
              call ivli(y(P(Q('u',j))),uj)
c             vj = x(Q('v',j))
              call ivli(y(P(Q('v',j))),vj)
c             sj = x(Q('s',j))
              call ivli(y(P(Q('s',j))),sj)

c             xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
              call scladd(-1.0d0,sj,tmp)
              call sclmlt(cte,tmp,tmp)
              call scladd(1.0d0,tmp,tmp)
              call mult(tmp,uj,xxj)

c             yyj = sj * vj
              call mult(sj,vj,yyj)

c             dist = ( xxi - xxj ) ** 2 + ( yyi - yyj ) ** 2
c             f = f - dist
              call sub(xxi,xxj,tmp)
              call mult(tmp,tmp,diffx)
              call sub(f,diffx,f)
              call sub(yyi,yyj,tmp)
              call mult(tmp,tmp,diffy)
              call sub(f,diffy,f)
          end do
      end do              

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalhi(n,y,hlin,hcol,hval,hnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n,hnnz

C     ARRAY ARGUMENTS
      integer hcol(*),hlin(*)
      double precision hval(*),y(2*n)

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     LOCAL SCALARS
      integer i,j
      double precision cte

C     LOCAL ARRAYS
      double precision diffx(2),diffy(2),dxxidui(2),dxxjduj(2),ui(2),
     +        uj(2),vi(2),vj(2),si(2),sj(2),xxi(2),xxj(2),yyi(2),yyj(2),
     +        tmp(2)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer P,Q
      external P,Q

      flag = 0

      cte = b ** 2 / a ** 2

      hnnz = 0

      do i = 1,nballs
          do j = i + 1,nballs
c             ui = x(Q('u',i))
              call ivli(y(P(Q('u',i))),ui)
c             vi = x(Q('v',i))
              call ivli(y(P(Q('v',i))),vi)
c             si = x(Q('s',i))
              call ivli(y(P(Q('s',i))),si)

c             xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
              call scladd(-1.0d0,si,tmp)
              call sclmlt(cte,tmp,tmp)
              call scladd(1.0d0,tmp,tmp)
              call mult(tmp,ui,xxi)

c             yyi = si * vi
              call mult(si,vi,yyi)

c             uj = x(Q('u',j))
              call ivli(y(P(Q('u',j))),uj)
c             vj = x(Q('v',j))
              call ivli(y(P(Q('v',j))),vj)
c             sj = x(Q('s',j))
              call ivli(y(P(Q('s',j))),sj)

c             xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
              call scladd(-1.0d0,sj,tmp)
              call sclmlt(cte,tmp,tmp)
              call scladd(1.0d0,tmp,tmp)
              call mult(tmp,uj,xxj)

c             yyj = sj * vj
              call mult(sj,vj,yyj)

c             diffx = 2.0d0 * ( xxi - xxj )
              call sub(xxi,xxj,diffx)
              call sclmlt(2.0d0,diffx,diffx)

c             diffy = 2.0d0 * ( yyi - yyj )
              call sub(yyi,yyj,diffy)
              call sclmlt(2.0d0,diffy,diffy)

c             dxxidui = ( 1.0d0 + ( si - 1.0d0 ) * cte )
              call scladd(-1.0d0,si,dxxidui)
              call sclmlt(cte,dxxidui,dxxidui)
              call scladd(1.0d0,dxxidui,dxxidui)

c             dxxjduj = ( 1.0d0 + ( sj - 1.0d0 ) * cte )
              call scladd(-1.0d0,sj,dxxjduj)
              call sclmlt(cte,dxxjduj,dxxjduj)
              call scladd(1.0d0,dxxjduj,dxxjduj)

              hlin(hnnz+1) = Q('u',i)
              hcol(hnnz+1) = Q('u',i)
c             hval(hnnz+1) = - 2.0d0 * dxxidui ** 2
              call mult(dxxidui,dxxidui,hval(P(hnnz+1)))
              call sclmlt(-2.0d0,hval(P(hnnz+1)),hval(P(hnnz+1)))

              hlin(hnnz+2) = Q('s',i)
              hcol(hnnz+2) = Q('u',i)
c             hval(hnnz+2) = - 2.0d0 * cte * ui * dxxidui - diffx * cte
              call mult(dxxidui,ui,hval(P(hnnz+2)))
              call sclmlt(-2.0d0*cte,hval(P(hnnz+2)),hval(P(hnnz+2)))
              call sclmlt(cte,diffx,tmp)
              call sub(hval(P(hnnz+2)),tmp,hval(P(hnnz+2)))

              hlin(hnnz+3) = Q('u',j)
              hcol(hnnz+3) = Q('u',i)
c             hval(hnnz+3) = 2.0d0 * dxxjduj * dxxidui
              call mult(dxxjduj,dxxidui,hval(P(hnnz+3)))
              call sclmlt(2.0d0,hval(P(hnnz+3)),hval(P(hnnz+3)))

              hlin(hnnz+4) = Q('s',j)
              hcol(hnnz+4) = Q('u',i)
c             hval(hnnz+4) = 2.0d0 * cte * uj * dxxidui 
              call mult(dxxidui,uj,hval(P(hnnz+4)))
              call sclmlt(2.0d0*cte,hval(P(hnnz+4)),hval(P(hnnz+4)))

              hlin(hnnz+5) = Q('v',i)
              hcol(hnnz+5) = Q('v',i)
c             hval(hnnz+5) = - 2.0d0 * si ** 2
              call mult(si,si,hval(P(hnnz+5)))
              call sclmlt(-2.0d0,hval(P(hnnz+5)),hval(P(hnnz+5)))

              hlin(hnnz+6) = Q('s',i)
              hcol(hnnz+6) = Q('v',i)
c             hval(hnnz+6) = - 4.0d0 * yyi +  2.0d0 * yyj
              call sclmlt(-4.0d0,yyi,hval(P(hnnz+6)))
              call sclmlt(2.0d0,yyj,tmp)
              call add(hval(P(hnnz+6)),tmp,hval(P(hnnz+6)))

              hlin(hnnz+7) = Q('v',j)
              hcol(hnnz+7) = Q('v',i)
c             hval(hnnz+7) = 2.0d0 * sj * si
              call mult(si,sj,hval(P(hnnz+7)))
              call sclmlt(2.0d0,hval(P(hnnz+7)),hval(P(hnnz+7)))

              hlin(hnnz+8) = Q('s',j)
              hcol(hnnz+8) = Q('v',i)
c             hval(hnnz+8) = 2.0d0 * vj * si
              call mult(vj,si,hval(P(hnnz+8)))
              call sclmlt(2.0d0,hval(P(hnnz+8)),hval(P(hnnz+8)))

              hlin(hnnz+9) = Q('s',i)
              hcol(hnnz+9) = Q('s',i)
c             hval(hnnz+9) = - 2.0d0 * ( cte ** 2 * ui ** 2 + vi ** 2 )
              call mult(vi,vi,hval(P(hnnz+9)))
              call mult(ui,ui,tmp)
              call sclmlt(cte**2,tmp,tmp)
              call add(hval(P(hnnz+9)),tmp,hval(P(hnnz+9)))
              call sclmlt(-2.0d0,hval(P(hnnz+9)),hval(P(hnnz+9)))

              hlin(hnnz+10) = Q('s',j)
              hcol(hnnz+10) = Q('s',i)
c             hval(hnnz+10) = 2.0d0 * ( cte ** 2 * ui * uj + vj * vi )
              call mult(vj,vi,hval(P(hnnz+10)))
              call mult(uj,ui,tmp)
              call sclmlt(cte**2,tmp,tmp)
              call add(hval(P(hnnz+10)),tmp,hval(P(hnnz+10)))
              call sclmlt(2.0d0,hval(P(hnnz+10)),hval(P(hnnz+10)))

              hlin(hnnz+11) = Q('s',i)
              hcol(hnnz+11) = Q('u',j)
c             hval(hnnz+11) = 2.0d0 * cte * ui * dxxjduj
              call mult(ui,dxxjduj,hval(P(hnnz+11)))
              call sclmlt(2.0d0*cte,hval(P(hnnz+11)),hval(P(hnnz+11)))

              hlin(hnnz+12) = Q('u',j)
              hcol(hnnz+12) = Q('u',j)
c             hval(hnnz+12) = - 2.0d0 * dxxjduj ** 2 
              call mult(dxxjduj,dxxjduj,hval(P(hnnz+12)))
              call sclmlt(-2.0d0*cte,hval(P(hnnz+12)),hval(P(hnnz+12)))

              hlin(hnnz+13) = Q('s',j)
              hcol(hnnz+13) = Q('u',j)
c             hval(hnnz+13) = - 2.0d0 * cte * uj * dxxjduj + diffx * cte
              call sclmlt(cte,diffx,hval(P(hnnz+13)))
              call mult(uj,dxxjduj,tmp)
              call sclmlt(-2.0d0*cte,tmp,tmp)
              call add(hval(P(hnnz+13)),tmp,hval(P(hnnz+13)))

              hlin(hnnz+14) = Q('s',i)
              hcol(hnnz+14) = Q('v',j)
c             hval(hnnz+14) = 2.0d0 * vi * sj
              call mult(vi,sj,hval(P(hnnz+14)))
              call sclmlt(2.0d0,hval(P(hnnz+14)),hval(P(hnnz+14)))

              hlin(hnnz+15) = Q('v',j)
              hcol(hnnz+15) = Q('v',j)
c             hval(hnnz+15) = - 2.0d0 * sj ** 2
              call mult(sj,sj,hval(P(hnnz+15)))
              call sclmlt(-2.0d0,hval(P(hnnz+15)),hval(P(hnnz+15)))

              hlin(hnnz+16) = Q('s',j)
              hcol(hnnz+16) = Q('v',j)
c             hval(hnnz+16) = 2.0d0 * yyi - 4.0d0 * yyj
              call sclmlt(-4.0d0,yyj,hval(P(hnnz+16)))
              call sclmlt(2.0d0,yyi,tmp)
              call add(hval(P(hnnz+16)),tmp,hval(P(hnnz+16)))

              hlin(hnnz+17) = Q('s',j)
              hcol(hnnz+17) = Q('s',j)
c             hval(hnnz+17) = - 2.0d0 * ( cte ** 2 * uj ** 2 + vj ** 2 )
              call mult(vj,vj,hval(P(hnnz+17)))
              call mult(uj,uj,tmp)
              call sclmlt(cte**2,tmp,tmp)
              call add(hval(P(hnnz+17)),tmp,hval(P(hnnz+17)))
              call sclmlt(-2.0d0,hval(P(hnnz+17)),hval(P(hnnz+17)))

              hnnz = hnnz + 17
          end do
      end do              

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalci(n,y,ind,c,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer ind,flag,n
      double precision c(2)

C     ARRAY ARGUMENTS
      double precision y(2*n)

C     LOCAL SCALARS
      integer i,j,k
      double precision cte

C     LOCAL ARRAYS
      double precision diffx(2),diffy(2),ui(2),uj(2),vi(2),vj(2),si(2),
     +        sj(2),xxi(2),xxj(2),yyi(2),yyj(2),tmp(2),tmp2(2),tmp3(2)

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer P,Q
      external P,Q

      flag = 0

      if ( ind .le. nballs ) then
          i  = ind

c         ui = x(Q('u',i))
          call ivli(y(P(Q('u',i))),ui)
c         vi = x(Q('v',i))
          call ivli(y(P(Q('v',i))),vi)

c         c = ( ui / a ) ** 2 + ( vi / b ) ** 2 - 1.0d0
          call sclmlt(1.0d0/a,ui,tmp)
          call mult(tmp,tmp,c)
          call sclmlt(1.0d0/b,vi,tmp)
          call mult(tmp,tmp,tmp2)
          call add(c,tmp2,c)
          call scladd(-1.0d0,c,c)

      else if ( ind .le. 2 * nballs ) then
          i  = ind - nballs

          cte = ( b ** 2 / a ** 2 ) ** 2

c         ui = x(Q('u',i))
          call ivli(y(P(Q('u',i))),ui)
c         vi = x(Q('v',i))
          call ivli(y(P(Q('v',i))),vi)
c         si = x(Q('s',i))
          call ivli(y(P(Q('s',i))),si)

c         c = r ** 2 - ( si - 1.0d0 ) ** 2 * ( cte * ui ** 2 + vi ** 2 )
          call mult(vi,vi,tmp)
          call mult(ui,ui,tmp2)
          call sclmlt(cte,tmp2,tmp2)
          call add(tmp2,tmp,tmp)
          call ivli(si,tmp)
          call scladd(-1.0d0,tmp,tmp)
          call mult(tmp,tmp,tmp3)
          call mult(tmp3,tmp2,tmp)
          call ivl1(r**2,c)
          call sub(c,tmp,c)

      else
          k = ind - 2 * nballs
          i = pairs(1,k)
          j = pairs(2,k)

          cte = b ** 2 / a ** 2

c         ui = x(Q('u',i))
          call ivli(y(P(Q('u',i))),ui)
c         vi = x(Q('v',i))
          call ivli(y(P(Q('v',i))),vi)
c         si = x(Q('s',i))
          call ivli(y(P(Q('s',i))),si)

c         xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
          call ivli(si,tmp)
          call scladd(-1.0d0,tmp,tmp)
          call sclmlt(cte,tmp,tmp)
          call scladd(1.0d0,tmp,tmp)
          call mult(tmp,ui,xxi)

c         yyi = si * vi
          call mult(si,vi,yyi)

c         uj = x(Q('u',j))
          call ivli(y(P(Q('u',j))),uj)
c         vj = x(Q('v',j))
          call ivli(y(P(Q('v',j))),vj)
c         sj = x(Q('s',j))
          call ivli(y(P(Q('s',j))),sj)

c         xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
          call ivli(sj,tmp)
          call scladd(-1.0d0,tmp,tmp)
          call sclmlt(cte,tmp,tmp)
          call scladd(1.0d0,tmp,tmp)
          call mult(tmp,uj,xxj)

c         yyj = sj * vj
          call mult(sj,vj,yyj)

c         c = 4 * r ** 2 - ( xxi - xxj ) ** 2 - ( yyi - yyj ) ** 2
          call ivl1(4.0d0*r**2,c)
          call ivli(xxi,tmp)
          call sub(tmp,xxj,tmp)
          call mult(tmp,tmp,diffx)
          call sub(c,diffx,c)
          call ivli(yyi,tmp)
          call sub(tmp,yyj,tmp)
          call mult(tmp,tmp,diffy)
          call sub(c,diffy,c)
      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evaljaci(n,y,ind,jcvar,jcval,jcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,n,jcnnz

C     ARRAY ARGUMENTS
      integer jcvar(n)
      double precision y(2*n),jcval(2*n)

C     LOCAL SCALARS
      integer i,j,k
      double precision cte

C     LOCAL ARRAYS
      double precision diffx(2),diffy(2),ui(2),uj(2),vi(2),vj(2),si(2),
     +        sj(2),xxi(2),xxj(2),yyi(2),yyj(2),tmp(2),tmp2(2)

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer P,Q
      external P,Q

      flag = 0

      if ( ind .le. nballs ) then
          i  = ind

c         ui = x(Q('u',i))
          call ivli(y(P(Q('u',i))),ui)
c         vi = x(Q('v',i))
          call ivli(y(P(Q('v',i))),vi)

          jcnnz = 2

          jcvar(1) = Q('u',i)
c         jcval(1) = 2.0d0 * ui / a ** 2
          call sclmlt(2.0d0/a**2,ui,jcval(P(1)))

          jcvar(2) = Q('v',i)
c         jcval(2) = 2.0d0 * vi / b ** 2
          call sclmlt(2.0d0/b**2,vi,jcval(P(2)))

      else if ( ind .le. 2 * nballs ) then
          i  = ind - nballs

          cte = ( b ** 2 / a ** 2 ) ** 2

c         ui = x(Q('u',i))
          call ivli(y(P(Q('u',i))),ui)
c         vi = x(Q('v',i))
          call ivli(y(P(Q('v',i))),vi)
c         si = x(Q('s',i))
          call ivli(y(P(Q('s',i))),si)

          jcnnz = 3

          jcvar(1) = Q('u',i)
c         jcval(1) = - ( si - 1.0d0 ) ** 2 * 2.0d0 * cte * ui
          call scladd(-1.0d0,si,tmp)
          call mult(tmp,tmp,jcval(P(1)))
          call sclmlt(-2.0d0*cte,ui,tmp)
          call mult(jcval(P(1)),tmp,jcval(P(1)))

          jcvar(2) = Q('v',i)
c         jcval(2) = - ( si - 1.0d0 ) ** 2 * 2.0d0 * vi
          call scladd(-1.0d0,si,tmp)
          call mult(tmp,tmp,jcval(P(2)))
          call sclmlt(-2.0d0,vi,tmp)
          call mult(jcval(P(2)),tmp,jcval(P(2)))

          jcvar(3) = Q('s',i)
c         jcval(3) = - 2.0d0 * ( si - 1.0d0 ) * 
c    +                 ( cte * ui ** 2 + vi ** 2 )
          call mult(vi,vi,tmp)
          call mult(ui,ui,tmp2)
          call sclmlt(cte,tmp2,tmp2)
          call add(tmp,tmp2,tmp)
          call scladd(-1.0d0,si,tmp2)
          call sclmlt(-2.0d0,tmp2,tmp2)
          call mult(tmp,tmp2,jcval(P(3)))

      else
          k = ind - 2 * nballs
          i = pairs(1,k)
          j = pairs(2,k)

          cte = b ** 2 / a ** 2

c         ui = x(Q('u',i))
          call ivli(y(P(Q('u',i))),ui)
c         vi = x(Q('v',i))
          call ivli(y(P(Q('v',i))),vi)
c         si = x(Q('s',i))
          call ivli(y(P(Q('s',i))),si)

c         xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
          call scladd(-1.0d0,si,tmp)
          call sclmlt(cte,tmp,tmp)
          call scladd(1.0d0,tmp,tmp)
          call mult(tmp,ui,xxi)

c         yyi = si * vi
          call mult(si,vi,yyi)

c         uj = x(Q('u',j))
          call ivli(y(P(Q('u',j))),uj)
c         vj = x(Q('v',j))
          call ivli(y(P(Q('v',j))),vj)
c         sj = x(Q('s',j))
          call ivli(y(P(Q('s',j))),sj)

c         xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
          call scladd(-1.0d0,sj,tmp)
          call sclmlt(cte,tmp,tmp)
          call scladd(1.0d0,tmp,tmp)
          call mult(tmp,uj,xxj)

c         yyj = sj * vj
          call mult(sj,vj,yyj)

c         diffx = 2.0d0 * ( xxi - xxj )
          call sub(xxi,xxj,diffx)
          call sclmlt(2.0d0,diffx,diffx)

c         diffy = 2.0d0 * ( yyi - yyj )
          call sub(yyi,yyj,diffy)
          call sclmlt(2.0d0,diffy,diffy)

          jcnnz = 6

          jcvar(1) = Q('u',i)
c         jcval(1) = - diffx * ( 1.0d0 + ( si - 1.0d0 ) * cte )
          call scladd(-1.0d0,si,tmp)
          call sclmlt(cte,tmp,tmp)
          call scladd(1.0d0,tmp,tmp)
          call mult(diffx,tmp,jcval(P(1)))
          call sclmlt(-1.0d0,jcval(P(1)),jcval(P(1)))

          jcvar(2) = Q('v',i)
c         jcval(2) = - diffy * si
          call mult(diffy,si,jcval(P(2)))
          call sclmlt(-1.0d0,jcval(P(2)),jcval(P(2)))

          jcvar(3) = Q('s',i)
c         jcval(3) = - diffx * cte * ui - diffy * vi
          call mult(diffy,vi,jcval(P(3)))
          call mult(diffx,ui,tmp)
          call sclmlt(cte,tmp,tmp)
          call add(jcval(P(3)),tmp,jcval(P(3)))
          call sclmlt(-1.0d0,jcval(P(3)),jcval(P(3)))

          jcvar(4) = Q('u',j)
c         jcval(4) = diffx * ( 1.0d0 + ( sj - 1.0d0 ) * cte )
          call scladd(-1.0d0,sj,tmp)
          call sclmlt(cte,tmp,tmp)
          call scladd(1.0d0,tmp,tmp)
          call mult(diffx,tmp,jcval(P(4)))

          jcvar(5) = Q('v',j)
c         jcval(5) = diffy * sj
          call mult(diffy,sj,jcval(P(5)))

          jcvar(6) = Q('s',j)
c         jcval(6) = diffx * cte * uj + diffy * vj
          call mult(diffy,vj,jcval(P(6)))
          call mult(diffx,uj,tmp)
          call sclmlt(cte,tmp,tmp)
          call add(jcval(P(6)),tmp,jcval(P(6)))

      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evalhci(n,y,ind,hclin,hccol,hcval,hcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,n,hcnnz

C     ARRAY ARGUMENTS
      integer hccol(*),hclin(*)
      double precision hcval(*),y(2*n)

C     LOCAL SCALARS
      integer i,j,k
      double precision cte

C     LOCAL ARRAYS
      double precision diffx(2),diffy(2),dxxidui(2),dxxjduj(2),ui(2),
     +        uj(2),vi(2),vj(2),si(2),sj(2),xxi(2),xxj(2),yyi(2),yyj(2),
     +        tmp(2)

C     PARAMETERS
      integer nballsmax,npairsmax
      parameter ( nballsmax = 100 )
      parameter ( npairsmax = nballsmax * ( nballsmax - 1 ) / 2 )

C     COMMON SCALARS
      integer nballs
      double precision a,b,r

C     COMMON ARRAYS
      integer pairs(2,npairsmax)

C     COMMON BLOCKS
      common /elldata/ pairs,a,b,r,nballs

C     EXTERNAL FUNCTIONS
      integer P,Q
      external P,Q

      flag = 0

      if ( ind .le. nballs ) then
          i  = ind

c         ui = x(Q('u',i))
          call ivli(y(P(Q('u',i))),ui)
c         vi = x(Q('v',i))
          call ivli(y(P(Q('v',i))),vi)

          hcnnz = 2

          hclin(1) = Q('u',i)
          hccol(1) = Q('u',i)
c         hcval(1) = 2.0d0 / a ** 2
          call ivl1(2.0d0/a**2,hcval(P(1)))

          hclin(2) = Q('v',i)
          hccol(2) = Q('v',i)
c         hcval(2) = 2.0d0 / b ** 2
          call ivl1(2.0d0/b**2,hcval(P(2)))

      else if ( ind .le. 2 * nballs ) then
          i  = ind - nballs

          cte = ( b ** 2 / a ** 2 ) ** 2

c         ui = x(Q('u',i))
          call ivli(y(P(Q('u',i))),ui)
c         vi = x(Q('v',i))
          call ivli(y(P(Q('v',i))),vi)
c         si = x(Q('s',i))
          call ivli(y(P(Q('s',i))),si)

          hcnnz = 5

          hclin(1) = Q('u',i)
          hccol(1) = Q('u',i)
c         hcval(1) = - ( si - 1.0d0 ) ** 2 * 2.0d0 * cte
          call scladd(-1.0d0,si,tmp)
          call mult(tmp,tmp,hcval(P(1)))
          call sclmlt(-2.0d0*cte,hcval(P(1)),hcval(P(1)))

          hclin(2) = Q('s',i)
          hccol(2) = Q('u',i)
c         hcval(2) = - 4.0d0 * ( si - 1.0d0 ) * cte * ui
          call scladd(-1.0d0,si,tmp)
          call sclmlt(-4.0d0*cte,tmp,tmp)
          call mult(tmp,ui,hcval(P(2)))

          hclin(3) = Q('v',i)
          hccol(3) = Q('v',i)
c         hcval(3) = - ( si - 1.0d0 ) ** 2 * 2.0d0
          call scladd(-1.0d0,si,tmp)
          call mult(tmp,tmp,hcval(P(3)))
          call sclmlt(-2.0d0,hcval(P(3)),hcval(P(3)))

          hclin(4) = Q('s',i)
          hccol(4) = Q('v',i)
c         hcval(4) = - 4.0d0 * ( si - 1.0d0 ) * vi
          call scladd(-1.0d0,si,tmp)
          call sclmlt(-4.0d0,tmp,tmp)
          call mult(tmp,vi,hcval(P(4)))

          hclin(5) = Q('s',i)
          hccol(5) = Q('s',i)
c         hcval(5) = - 2.0d0 * ( cte * ui ** 2 + vi ** 2 )
          call mult(vi,vi,hcval(P(5)))
          call mult(ui,ui,tmp)
          call sclmlt(cte,tmp,tmp)
          call add(hcval(P(5)),tmp,hcval(P(5)))
          call sclmlt(-2.0d0,hcval(P(5)),hcval(P(5)))

      else
          k = ind - 2 * nballs
          i = pairs(1,k)
          j = pairs(2,k)

          cte = b ** 2 / a ** 2

c         ui = x(Q('u',i))
          call ivli(y(P(Q('u',i))),ui)
c         vi = x(Q('v',i))
          call ivli(y(P(Q('v',i))),vi)
c         si = x(Q('s',i))
          call ivli(y(P(Q('s',i))),si)

c         xxi = ( 1.0d0 + ( si - 1.0d0 ) * cte ) * ui
          call scladd(-1.0d0,si,tmp)
          call sclmlt(cte,tmp,tmp)
          call scladd(1.0d0,tmp,tmp)
          call mult(tmp,ui,xxi)

c         yyi = si * vi
          call mult(si,vi,yyi)

c         uj = x(Q('u',j))
          call ivli(y(P(Q('u',j))),uj)
c         vj = x(Q('v',j))
          call ivli(y(P(Q('v',j))),vj)
c         sj = x(Q('s',j))
          call ivli(y(P(Q('s',j))),sj)

c         xxj = ( 1.0d0 + ( sj - 1.0d0 ) * cte ) * uj
          call scladd(-1.0d0,sj,tmp)
          call sclmlt(cte,tmp,tmp)
          call scladd(1.0d0,tmp,tmp)
          call mult(tmp,uj,xxj)

c         yyj = sj * vj
          call mult(sj,vj,yyj)

c         diffx = 2.0d0 * ( xxi - xxj )
          call sub(xxi,xxj,diffx)
          call sclmlt(2.0d0,diffx,diffx)

c         diffy = 2.0d0 * ( yyi - yyj )
          call sub(yyi,yyj,diffy)
          call sclmlt(2.0d0,diffy,diffy)

c         dxxidui = ( 1.0d0 + ( si - 1.0d0 ) * cte )
          call scladd(-1.0d0,si,dxxidui)
          call sclmlt(cte,dxxidui,dxxidui)
          call scladd(1.0d0,dxxidui,dxxidui)

c         dxxjduj = ( 1.0d0 + ( sj - 1.0d0 ) * cte )
          call scladd(-1.0d0,sj,dxxjduj)
          call sclmlt(cte,dxxjduj,dxxjduj)
          call scladd(1.0d0,dxxjduj,dxxjduj)

          hcnnz = 17

          hclin(1) = Q('u',i)
          hccol(1) = Q('u',i)
c         hcval(1) = - 2.0d0 * dxxidui ** 2
          call mult(dxxidui,dxxidui,hcval(P(1)))
          call sclmlt(-2.0d0,hcval(P(1)),hcval(P(1)))

          hclin(2) = Q('s',i)
          hccol(2) = Q('u',i)
c         hcval(2) = - 2.0d0 * cte * ui * dxxidui - diffx * cte
          call mult(dxxidui,ui,hcval(P(2)))
          call sclmlt(-2.0d0*cte,hcval(P(2)),hcval(P(2)))
          call sclmlt(cte,diffx,tmp)
          call sub(hcval(P(2)),tmp,hcval(P(2)))

          hclin(3) = Q('u',j)
          hccol(3) = Q('u',i)
c         hcval(3) = 2.0d0 * dxxjduj * dxxidui
          call mult(dxxjduj,dxxidui,hcval(P(3)))
          call sclmlt(2.0d0,hcval(P(3)),hcval(P(3)))

          hclin(4) = Q('s',j)
          hccol(4) = Q('u',i)
c         hcval(4) = 2.0d0 * cte * uj * dxxidui 
          call mult(dxxidui,uj,hcval(P(4)))
          call sclmlt(2.0d0*cte,hcval(P(4)),hcval(P(4)))

          hclin(5) = Q('v',i)
          hccol(5) = Q('v',i)
c         hcval(5) = - 2.0d0 * si ** 2
          call mult(si,si,hcval(P(5)))
          call sclmlt(-2.0d0,hcval(P(5)),hcval(P(5)))

          hclin(6) = Q('s',i)
          hccol(6) = Q('v',i)
c         hcval(6) = - 4.0d0 * yyi +  2.0d0 * yyj
          call sclmlt(-4.0d0,yyi,hcval(P(6)))
          call sclmlt(2.0d0,yyj,tmp)
          call add(hcval(P(6)),tmp,hcval(P(6)))

          hclin(7) = Q('v',j)
          hccol(7) = Q('v',i)
c         hcval(7) = 2.0d0 * sj * si
          call mult(si,sj,hcval(P(7)))
          call sclmlt(2.0d0,hcval(P(7)),hcval(P(7)))

          hclin(8) = Q('s',j)
          hccol(8) = Q('v',i)
c         hcval(8) = 2.0d0 * vj * si
          call mult(vj,si,hcval(P(8)))
          call sclmlt(2.0d0,hcval(P(8)),hcval(P(8)))

          hclin(9) = Q('s',i)
          hccol(9) = Q('s',i)
c         hcval(9) = - 2.0d0 * ( cte ** 2 * ui ** 2 + vi ** 2 )
          call mult(vi,vi,hcval(P(9)))
          call mult(ui,ui,tmp)
          call sclmlt(cte**2,tmp,tmp)
          call add(hcval(P(9)),tmp,hcval(P(9)))
          call sclmlt(-2.0d0,hcval(P(9)),hcval(P(9)))

          hclin(10) = Q('s',j)
          hccol(10) = Q('s',i)
c         hcval(10) = 2.0d0 * ( cte ** 2 * ui * uj + vj * vi )
          call mult(vj,vi,hcval(P(10)))
          call mult(uj,ui,tmp)
          call sclmlt(cte**2,tmp,tmp)
          call add(hcval(P(10)),tmp,hcval(P(10)))
          call sclmlt(2.0d0,hcval(P(10)),hcval(P(10)))

          hclin(11) = Q('s',i)
          hccol(11) = Q('u',j)
c         hcval(11) = 2.0d0 * cte * ui * dxxjduj
          call mult(ui,dxxjduj,hcval(P(11)))
          call sclmlt(2.0d0*cte,hcval(P(11)),hcval(P(11)))

          hclin(12) = Q('u',j)
          hccol(12) = Q('u',j)
c         hcval(12) = - 2.0d0 * dxxjduj ** 2 
          call mult(dxxjduj,dxxjduj,hcval(P(12)))
          call sclmlt(-2.0d0*cte,hcval(P(12)),hcval(P(12)))

          hclin(13) = Q('s',j)
          hccol(13) = Q('u',j)
c         hcval(13) = - 2.0d0 * cte * uj * dxxjduj + diffx * cte
          call sclmlt(cte,diffx,hcval(P(13)))
          call mult(uj,dxxjduj,tmp)
          call sclmlt(-2.0d0*cte,tmp,tmp)
          call add(hcval(P(13)),tmp,hcval(P(13)))

          hclin(14) = Q('s',i)
          hccol(14) = Q('v',j)
c         hcval(14) = 2.0d0 * vi * sj
          call mult(vi,sj,hcval(P(14)))
          call sclmlt(2.0d0,hcval(P(14)),hcval(P(14)))

          hclin(15) = Q('v',j)
          hccol(15) = Q('v',j)
c         hcval(15) = - 2.0d0 * sj ** 2
          call mult(sj,sj,hcval(P(15)))
          call sclmlt(-2.0d0,hcval(P(15)),hcval(P(15)))

          hclin(16) = Q('s',j)
          hccol(16) = Q('v',j)
c         hcval(16) = 2.0d0 * yyi - 4.0d0 * yyj
          call sclmlt(-4.0d0,yyj,hcval(P(16)))
          call sclmlt(2.0d0,yyi,tmp)
          call add(hcval(P(16)),tmp,hcval(P(16)))

          hclin(17) = Q('s',j)
          hccol(17) = Q('s',j)
c         hcval(17) = - 2.0d0 * ( cte ** 2 * uj ** 2 + vj ** 2 )
          call mult(vj,vj,hcval(P(17)))
          call mult(uj,uj,tmp)
          call sclmlt(cte**2,tmp,tmp)
          call add(hcval(P(17)),tmp,hcval(P(17)))
          call sclmlt(-2.0d0,hcval(P(17)),hcval(P(17)))
      end if

      end