C     =================================================================
C     File: prodpl0reflin.f
C     =================================================================

C     A production planning problem in the computer industry, 
C     from CUTEr (Source: L. Escudero, private communication, 1991.)

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

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

      implicit none

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

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

      include "prodpl0reflin.par"

C     LOCAL SCALARS
      integer i,ind,j,k

C     LOCAL ARRAYS
      double precision bb(imax,jmax),cc(jmax)

C     EXTERNAL FUNCTIONS
      integer II

C     Number of variables
      n = 3 * ( imax * jmax ) + 2 * ( jmax - 1 )

C     Lower and upper bounds
      do i = 1,imax
          do j = 1,jmax
              k = II('x',i,j)
              l(k) =  0.0d0
              u(k) = 10.0d0
              k = II('y',i,j)
              l(k) =  0.0d0
              u(k) = 10.0d0
              k = II('i',i,j)
              l(k) =  0.0d0
              u(k) = 10.0d0
          end do
      end do

      do j = 1,jmax - 1
          k = II('w',1,j)
          l(k) = -10.0d0
          u(k) =  30.0d0
          k = II('w',2,j)
          l(k) =   0.0d0
          u(k) =  20.0d0
      end do

C     Initial point
      do i = 1,n
          x(i) = 0.0d0
      end do

C     Number of constraints
      m   = jmax - 1
      ml  = 2 * ( jmax - 1 ) + imax * jmax + jmax
      meq = 2 * ( jmax - 1 ) + imax * jmax

C     Lagrange multipliers approximation and initial penalty parameter. 
      do i = 1,m
          lambda(i) =  0.0d0
          rho(i)    = 10.0d0
      end do

C     Lower-level constraints are all nonlinear inequalities

      do i = 1,m
          equatn(i) = .false.
          linear(i) = .false.
      end do

C     Linear constraints

      ind  = 0
      annz = 0

C     Linear equality constraints from reformulation

      do j = 1,jmax - 1
          ind = ind + 1

          b(ind) = 0.0d0

          alin(annz+1) = ind
          acol(annz+1) = II('x',3,j+1)
          aval(annz+1) = 1.0d0      

          alin(annz+2) = ind
          acol(annz+2) = II('x',4,j+1)
          aval(annz+2) = 1.0d0      

          alin(annz+3) = ind
          acol(annz+3) = II('x',3,j)
          aval(annz+3) = - 1.0d0   

          alin(annz+4) = ind
          acol(annz+4) = II('x',4,j)
          aval(annz+4) = 1.0d0      

          alin(annz+5) = ind
          acol(annz+5) = II('w',1,j)
          aval(annz+5) = - 1.0d0      

          annz = annz + 5
      end do

      do j = 1,jmax - 1
          ind = ind + 1

          b(ind) = 0.0d0

          alin(annz+1) = ind
          acol(annz+1) = II('x',3,j)
          aval(annz+1) = 1.0d0      

          alin(annz+2) = ind
          acol(annz+2) = II('x',4,j)
          aval(annz+2) = 1.0d0      

          alin(annz+3) = ind
          acol(annz+3) = II('w',2,j)
          aval(annz+3) = - 1.0d0      

          annz = annz + 3
      end do

C     Linear equality constraints

      bb(1,1) = - 1.0d0
      bb(2,1) = - 1.0d0
      bb(3,1) = - 1.0d0
      bb(4,1) = - 1.0d0
      bb(1,2) = - 2.667d0
      bb(2,2) = - 1.667d0
      bb(3,2) = - 2.667d0
      bb(4,2) = - 3.333d0
      bb(1,3) = - 2.667d0
      bb(2,3) = - 2.0d0
      bb(3,3) = - 3.0d0
      bb(4,3) = - 3.0d0
      bb(1,4) = - 2.667d0
      bb(2,4) = - 2.667d0
      bb(3,4) = - 2.667d0
      bb(4,4) = - 2.667d0
      bb(1,5) = - 2.667d0
      bb(2,5) = - 2.333d0
      bb(3,5) = - 2.333d0
      bb(4,5) = - 2.333d0

      do i = 1,imax
          do j = 1,jmax
              ind = ind + 1

              b(ind) = bb(i,j)

              alin(annz+1) = ind
              acol(annz+1) = II('x',i,j)
              aval(annz+1) = - 1.0d0      

              alin(annz+2) = ind
              acol(annz+2) = II('i',i,j)
              aval(annz+2) = 1.0d0      

              alin(annz+3) = ind
              acol(annz+3) = II('y',i,j)
              aval(annz+3) = - 1.0d0      

              annz = annz + 3

              if ( j .ne. 1 ) then
                  alin(annz+1) = ind
                  acol(annz+1) = II('i',i,j-1)
                  aval(annz+1) = - 1.0d0      

                  annz = annz + 1
              end if
          end do
      end do

C     Linear inequality constraints

      cc(1) = -    3.0d0
      cc(2) = -    6.0d0
      cc(3) = -   10.0d0
      cc(4) = - 2000.0d0
      cc(5) = -   18.0d0

      do j = 1,jmax
          ind = ind + 1

          b(ind) = cc(j)

          do i = 1,imax
              alin(annz+i) = ind
              acol(annz+i) = II('x',i,j)
              aval(annz+i) = - 1.0d0
          end do

          annz = annz + imax
      end do

      fknown = 6.091923647619048d+01

      end

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

      subroutine evalf(n,x,f,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n
      double precision f

C     ARRAY ARGUMENTS
      double precision x(n)

C     This subroutine must compute the objective function. For achieving 
C     this objective YOU MUST MODIFY it according to your problem. See 
C     below where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     On Return
C
C     f        double precision,
C              objective function value at x,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the objective 
C              function. (For example, trying to compute the square root 
C              of a negative number, dividing by zero or a very small 
C              number, etc.) If everything was o.k. you must set it 
C              equal to zero.

      include "prodpl0reflin.par"

C     LOCAL SCALARS
      integer i,j

C     EXTERNAL FUNCTIONS
      integer II

      flag = 0

      f = 0.0d0

      do i = 1,imax
          do j = 1,jmax
              f = f + float(i) * x(II('i',i,j))
          end do
      end do

      do j = 1,jmax
          f = f + 2.0d0 * x(II('y',1,j))
          f = f + 3.0d0 * x(II('y',2,j))
          f = f + 2.0d0 * x(II('y',3,j))
          f = f + 5.0d0 * x(II('y',4,j))
      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     This subroutine must compute the gradient vector of the objective 
C     function. For achieving these objective YOU MUST MODIFY it in the 
C     way specified below. However, if you decide to use numerical 
C     derivatives (we dont encourage this option at all!) you dont need
C     to modify evalg.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     On Return
C
C     g        double precision g(n),
C              gradient vector of the objective function evaluated at x,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of any component 
C              of the gradient vector. (For example, trying to compute 
C              the square root of a negative number, dividing by zero or 
C              a very small number, etc.) If everything was o.k. you 
C              must set it equal to zero.

      include "prodpl0reflin.par"

C     LOCAL SCALARS
      integer i,j

C     EXTERNAL FUNCTIONS
      integer II

      flag = 0

      do i = 1,imax
          do j = 1,jmax
              g(II('x',i,j)) = 0.0d0
          end do
      end do

      do i = 1,imax
          do j = 1,jmax
              g(II('i',i,j)) = float(i)
          end do
      end do

      do j = 1,jmax
          g(II('y',1,j)) = 2.0d0
          g(II('y',2,j)) = 3.0d0
          g(II('y',3,j)) = 2.0d0
          g(II('y',4,j)) = 5.0d0
      end do

      do j = 1,jmax - 1
          g(II('w',1,j)) = 0.0d0
          g(II('w',2,j)) = 0.0d0
      end do

      end

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

      implicit none

C     SCALAR ARGUMENTS
      integer flag,n,hnnz

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

C     This subroutine might compute the Hessian matrix of the objective 
C     function. For achieving this objective YOU MAY MODIFY it according 
C     to your problem. To modify this subroutine IS NOT MANDATORY. See 
C     below where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     On Return
C
C     hnnz     integer,
C              number of perhaps-non-null elements of the computed 
C              Hessian,
C
C     hlin     integer hlin(hnnz),
C              see below,
C
C     hcol     integer hcol(hnnz),
C              see below,
C
C     hval     double precision hval(hnnz),
C              the non-null value of the (hlin(k),hcol(k)) position 
C              of the Hessian matrix of the objective function must 
C              be saved at hval(k). Just the lower triangular part of
C              Hessian matrix must be computed,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the Hessian
C              matrix of the objective funtion. (For example, trying 
C              to compute the square root of a negative number, 
C              dividing by zero or a very small number, etc.) If 
C              everything was o.k. you must set it equal to zero.

      flag = 0

      hnnz = 0

      end

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

      implicit none

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

C     ARRAY ARGUMENTS
      double precision x(n)

C     This subroutine must compute the ind-th constraint of your 
C     problem. For achieving this objective YOU MUST MOFIFY it 
C     according to your problem. See below the places where your 
C     modifications must be inserted.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     ind      integer,
C              index of the constraint to be computed,
C
C     On Return
C
C     c        double precision,
C              ind-th constraint evaluated at x,
C
C     flag     integer
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the 
C              constraint. (For example, trying to compute the square 
C              root of a negative number, dividing by zero or a very 
C              small number, etc.) If everything was o.k. you must set 
C              it equal to zero.

C     EXTERNAL FUNCTIONS
      integer II

      flag = 0

      c = x(II('w',1,ind)) ** 2 - 0.01d0 * x(II('w',2,ind)) ** 2

      end

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

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,n,jcnnz

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

C     This subroutine must compute the gradient of the constraint ind. 
C     For achieving these objective YOU MUST MODIFY it in the way 
C     specified below.
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     ind      integer,
C              index of the constraint whose gradient will be computed,
C
C     On Return
C
C     jcnnz   integer,
C              number of perhaps-non-null elements of the computed 
C              gradient,
C
C     jcvar   integer jcvar(jcnnz),
C              see below,
C
C     jcval   double precision jcval(jcnnz),
C              the non-null value of the partial derivative of the 
C              ind-th constraint with respect to the jcvar(k)-th 
C              variable must be saved at jcval(k).
C
C     flag     integer
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the 
C              constraint. (For example, trying to compute the square 
C              root of a negative number, dividing by zero or a very 
C              small number, etc.) If everything was o.k. you must set 
C              it equal to zero.

C     LOCAL SCALARS
      integer i1,i2

C     EXTERNAL FUNCTIONS
      integer II

      flag = 0

      i1 = II('w',1,ind)
      i2 = II('w',2,ind)

      jcnnz = 2

      jcvar(1) = i1
      jcval(1) = 2.0d0 * x(i1)

      jcvar(2) = i2
      jcval(2) = - 0.02d0 * x(i2)

      end

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

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,n,hcnnz

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

C     This subroutine might compute the Hessian matrix of the ind-th
C     constraint. For achieving this objective YOU MAY MODIFY it 
C     according to your problem. To modify this subroutine IS NOT 
C     MANDATORY. See below where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     ind      integer,
C              index of the constraint whose Hessian will be computed,
C
C     On Return
C
C     hcnnz    integer,
C              number of perhaps-non-null elements of the computed 
C              Hessian,
C
C     hclin    integer hclin(hcnnz),
C              see below,
C
C     hccol    integer hccol(hcnnz),
C              see below,
C
C     hcval    double precision hcval(hcnnz),
C              the non-null value of the (hclin(k),hccol(k)) position 
C              of the Hessian matrix of the ind-th constraint must 
C              be saved at hcval(k). Just the lower triangular part of
C              Hessian matrix must be computed,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if 
C              some error ocurred during the evaluation of the Hessian
C              matrix of the ind-th constraint. (For example, trying 
C              to compute the square root of a negative number, 
C              dividing by zero or a very small number, etc.) If 
C              everything was o.k. you must set it equal to zero.

C     LOCAL SCALARS
      integer i1,i2

C     EXTERNAL FUNCTIONS
      integer II

      flag = 0

      i1 = II('w',1,ind)
      i2 = II('w',2,ind)

      hcnnz = 2

      hclin(1) = i1
      hccol(1) = i1
      hcval(1) = 2.0d0

      hclin(2) = i2
      hccol(2) = i2
      hcval(2) = - 0.02d0

      end

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

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

      implicit none

C     This subroutine might compute the product of the Hessian of the
C     Lagrangian times vector p (just the Hessian of the objective
C     function in the unconstrained or bound-constrained case).
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              current point,
C
C     m        integer,
C              number of constraints,
C
C     lambda   double precision lambda(m),
C              vector of Lagrange multipliers,
C
C     p        double precision p(n),
C              vector of the matrix-vector product,
C
C     goth     logical,
C              can be used to indicate if the Hessian matrices were
C              computed at the current point. It is set to .false.
C              by the optimization method every time the current
C              point is modified. Sugestion: if its value is .false.
C              then compute the Hessians, save them in a common
C              structure and set goth to .true.. Otherwise, just use
C              the Hessians saved in the common block structure,
C
C     On Return
C
C     hp       double precision hp(n),
C              Hessian-vector product,
C
C     goth     logical,
C              see above,
C
C     flag     integer,
C              You must set it to any number different of 0 (zero) if
C              some error ocurred during the evaluation of the
C              Hessian-vector product. (For example, trying to compute
C              the square root of a negative number, dividing by zero
C              or a very small number, etc.) If everything was o.k. you
C              must set it equal to zero.

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,equatn,linear)

      implicit none

C     SCALAR ARGUMENTS
      integer m,n

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

C     This subroutine can be used to do some extra job after the solver
C     has found the solution,like some extra statistics, or to save the
C     solution in some special format or to draw some graphical
C     representation of the solution. If the information given by the
C     solver is enough for you then leave the body of this subroutine
C     empty.
C     
C     Parameters of the subroutine:
C
C     The paraemters of this subroutine are the same parameters of
C     subroutine inip. But in this subroutine there are not output
C     parameter. All the parameters are input parameters.

      end

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

      integer function II(vtype,i,j)

C     SCALAR ARGUMENTS
      integer i,j
      character vtype
 
      include "prodpl0reflin.par"

C     Original variables

C     x00101 i00101 y00101
C     x00201 i00201 y00201
C     x00301 i00301 y00301
C     x00401 i00401 y00401
C     x00102 i00102 y00102
C     x00202 i00202 y00202
C     x00302 i00302 y00302
C     x00402 i00402 y00402
C     x00103 i00103 y00103
C     x00203 i00203 y00203
C     x00303 i00303 y00303
C     x00403 i00403 y00403
C     x00104 i00104 y00104
C     x00204 i00204 y00204
C     x00304 i00304 y00304
C     x00404 i00404 y00404
C     x00105 i00105 y00105
C     x00205 i00205 y00205
C     x00305 i00305 y00305
C     x00405 i00405 y00405

      if ( vtype .eq. 'x' .or. vtype .eq. 'X' ) then
          II = ( i - 1 ) * jmax + j

      else if ( vtype .eq. 'i' .or. vtype .eq. 'I' ) then
          II = imax * jmax + ( i - 1 ) * jmax + j

      else if ( vtype .eq. 'y' .or. vtype .eq. 'Y' ) then
          II = 2 * ( imax * jmax ) + ( i - 1 ) * jmax + j

      else if ( vtype .eq. 'w' .or. vtype .eq. 'W' ) then
          II = 3 * ( imax * jmax ) + ( i - 1 ) * ( jmax  - 1 ) + j
      end if

      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)

      include "prodpl0reflin.par"

C     LOCAL SCALARS
      integer i,i1,i2,itimes,j
      double precision tmp1,tmp2,tmpl2,tmpu2

C     EXTERNAL FUNCTIONS
      integer II

      mla   = ml
      annza = annz

C     Remove the comment below to avoid using linear relaxations of 
C     the penalized constraints
c     return

C     f(x) - f <= 0

      mla = mla + 1

      ind(mla) = 0

      b(mla) = - min( f, fknown )

      do i = 1,imax
          do j = 1,jmax
              annza = annza + 1

              alin(annza) = mla
              acol(annza) = II('i',i,j)
              aval(annza) = - float(i)
          end do
      end do

      do j = 1,jmax
          annza = annza + 1

          alin(annza) = mla
          acol(annza) = II('y',1,j)
          aval(annza) = - 2.0d0

          annza = annza + 1

          alin(annza) = mla
          acol(annza) = II('y',2,j)
          aval(annza) = - 3.0d0

          annza = annza + 1

          alin(annza) = mla
          acol(annza) = II('y',3,j)
          aval(annza) = - 2.0d0

          annza = annza + 1

          alin(annza) = mla
          acol(annza) = II('y',4,j)
          aval(annza) = - 5.0d0
      end do

C     For j = 1,jmax - 1:

C     c = x(II('w',1,j)) ** 2 - 0.01d0 * x(II('w',2,j)) ** 2.

C     Several linear underestimators for x(II('w',1,ind)) ** 2 and the
C     and the linear understimation of concave functions for - 0.01d0 * 
C     x(II('w',2,ind)) ** 2.

      do j = 1,jmax - 1
          i1 = II('w',1,j)
          i2 = II('w',2,j)

C         tmpl2 + tmp2 * ( x(i2) - l(i2) ) <= - x(i2) ** 2

          tmpl2 = - l(i2) ** 2
          tmpu2 = - u(i2) ** 2
          tmp2  = ( tmpu2 - tmpl2 ) / ( u(i2) - l(i2) )

          itimes = 3

          do i = 1,itimes
              mla = mla + 1

              ind(mla) = j

              annza = annza + 1

              alin(annza) = mla
              acol(annza) = i2
              aval(annza) = - 0.01d0 * tmp2

C             - tmp1 ** 2 + 2.0d0 * tmp1 * x(i1) <= x(i1) ** 2

              tmp1 = l(i1) + ( u(i1) - l(i1) ) * ( i - 1 ) / (itimes-1)

              annza = annza + 1

              alin(annza) = mla
              acol(annza) = i1
              aval(annza) = - 2.0d0 * tmp1

              b(mla) = 0.01d0 * ( tmpl2 - tmp2 * l(i2) ) - tmp1 ** 2
          end do
      end do

      end

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

      subroutine setlv(n,linvar)

      implicit none

C     SCALAR ARGUMENTS
      integer n

C     ARRAY ARGUMENTS
      logical linvar(n)

      include "prodpl0reflin.par"

C     LOCAL SCALARS
      integer i,j

C     EXTERNAL FUNCTIONS
      integer II

      do i = 1,imax
          do j = 1,jmax
              linvar(II('x',i,j)) = .true.
              linvar(II('i',i,j)) = .true.
              linvar(II('y',i,j)) = .true.
          end do
      end do

      do j = 1,jmax - 1
          linvar(II('w',1,j)) = .false.
          linvar(II('w',2,j)) = .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)

      include "prodpl0reflin.par"

C     LOCAL SCALARS
      integer i,j
      double precision r

C     LOCAL ARRAYS
      double precision tmp(2)

C     EXTERNAL FUNCTIONS
      integer II,P

      flag = 0

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

      do i = 1,imax
          do j = 1,jmax
C             f = f + i * y(P(II('i',i,j)))
              r = float(i)
              call sclmlt(r,y(P(II('i',i,j))),tmp)
              call add(tmp,f,f)
          end do
      end do

      do j = 1,jmax
C         f = f + 2.0d0 * y(P(II('y',1,j)))
          call sclmlt(2.0d0,y(P(II('y',1,j))),tmp)
          call add(tmp,f,f)

C         f = f + 3.0d0 * y(P(II('y',2,j)))
          call sclmlt(3.0d0,y(P(II('y',2,j))),tmp)
          call add(tmp,f,f)

C         f = f + 2.0d0 * y(P(II('y',3,j)))
          call sclmlt(2.0d0,y(P(II('y',3,j))),tmp)
          call add(tmp,f,f)

C         f = f + 5.0d0 * y(P(II('y',4,j)))
          call sclmlt(5.0d0,y(P(II('y',4,j))),tmp)
          call add(tmp,f,f)
      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)

      flag = 0

      hnnz = 0

      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 ARRAYS
      double precision tmp1(2),tmp2(2)

C     EXTERNAL FUNCTIONS
      integer II,P

      flag = 0

C     c = x(II('w',1,ind)) ** 2 - 0.01d0 * x(II('w',2,ind)) ** 2

      call power(y(P(II('w',1,ind))),2,tmp1)
      call power(y(P(II('w',2,ind))),2,tmp2)
      call sclmlt(0.01d0,tmp2,tmp2)
      call sub(tmp1,tmp2,c)

      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 i1,i2

C     EXTERNAL FUNCTIONS
      integer II,P

      flag = 0

      i1 = II('w',1,ind)
      i2 = II('w',2,ind)

      jcnnz = 2

      jcvar(1) = i1
C     jcval(1) = 2.0d0 * x(i1) 
      call sclmlt(2.0d0,y(P(i1)),jcval(P(1)))

      jcvar(2) = i2
C     jcval(2) = - 0.02d0 * x(i2)
      call sclmlt(-0.02d0,y(P(i2)),jcval(P(2)))

      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 i1,i2

C     EXTERNAL FUNCTIONS
      integer II,P

      flag = 0

      i1 = II('w',1,ind)
      i2 = II('w',2,ind)

      hcnnz = 2

      hclin(1) = i1
      hccol(1) = i1
C     hcval(1) = 2.0d0
      call ivl1(2.0d0,hcval(P(1)))

      hclin(2) = i2
      hccol(2) = i2
C     hcval(2) = - 0.02d0
      call ivl1(-0.02d0,hcval(P(2)))

      end