C     =================================================================
C     File: problem02b-without-reform.f
C     =================================================================

C     =================================================================
C     Module: Subroutines that define the problem
C     =================================================================

C     Last update of any of the component of this module: 
 
C     November 6, 2006.

C     =================================================================

C     Problem 2b from E. G. Birgin, C. A. Floudas and J. M. Martinez,
C     Global minimization using an Augmented Lagrangian method with
C     variable lower-level constraints, submitted.

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

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

      implicit none

C     SCALAR ARGUMENTS
      integer annz,m,ml,meq,n

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

C     This subroutine must set some problem data. For achieving this 
C     objective YOU MUST MODIFY it according to your problem. See below 
C     where your modifications must be inserted.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     This subroutine has no input parameters.
C
C     On Return
C
C     n        integer,
C              number of variables,
C
C     x        double precision x(n),
C              initial point,
C
C     l        double precision l(n),
C              lower bounds on x,
C
C     u        double precision u(n),
C              upper bounds on x,
C
C     m        integer,
C              number of constraints (excluding the bounds),
C
C     lambda   double precision lambda(m),
C              initial estimation of the Lagrange multipliers,
C
C     rho      double precision rho(m),
C              initial penalty parameters.
C
C     equatn   logical equatn(m)
C              for each constraint j, set equatn(j) = .true. if it is an 
C              equality constraint of the form c_j(x) = 0, and set 
C              equatn(j) = .false. if it is an inequality constraint of 
C              the form c_j(x) <= 0,
C
C     linear   logical linear(m)
C              for each constraint j, set linear(j) = .true. if it is a 
C              linear constraint, and set linear(j) = .false. if it is a
C              nonlinear constraint.

C     COMMON SCALARS
      double precision c1,c2

C     LOCAL SCALARS
      integer i
      double precision seed

C     COMMON BLOCKS
      common /probdata/ c1,c2

C     EXTERNAL FUNCTIONS
      double precision drand

C     Set problem data
C     write(*,*) 'Enter with c1 and c2'
C     read(*,*) c1,c2
      c1 = 16.0d0
      c2 = 600.0d0

C     write(*,*) 'Seed for the random initial point: '
C     read(*,*) seed
      seed = 123456.0d0

C     Number of variables
      n = 9

C     Lower and upper bounds
      do i = 1,n
          l(i) = 0.0d0
      end do

      u(1) = c2
      u(2) = 200.0d0
      do i = 3,n
          u(i) = 500.0d0
      end do

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

C     Number of constraints (equalities plus inequalities)

      m = 3

C     Lagrange multipliers approximation.

      do i = 1,m
          lambda(i) = 0.0d0
      end do

C     Initial penalty parameters

      do i = 1,m
          rho(i) = 10.0d0
      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.

      equatn(1) = .false.
      equatn(2) = .false.
      equatn(3) = .true.

C     For each constraint i, set linear(i) = .true. if it is a linear
C     constraint, otherwise set linear(i) = .false.

      linear(1) = .false.
      linear(2) = .false.
      linear(3) = .false.

C     Set the sparse matrix and the right hand side of the linear
C     constraints of the form ai^t x = bi for i = 1, ..., meq and
C     ai^t x >= bi for i = meq + 1,ml. These linear constraints will
C     be lower-level constraints. All the previous constraints will
C     be upper-lever constraints. The suggestions is to consider the
C     linear constraints as lower-level constraints, but it is not
C     mandatory.

      ml  = 3
      meq = 3

      annz = 10

C     x(8) + x(9) - x(3) - x(4) = 0

      alin(1) = 1
      acol(1) = 8
      aval(1) =   1.0d0

      alin(2) = 1
      acol(2) = 9
      aval(2) =   1.0d0

      alin(3) = 1
      acol(3) = 3
      aval(3) = - 1.0d0

      alin(4) = 1
      acol(4) = 4
      aval(4) = - 1.0d0

      b(1) = 0.0d0

C     x(1) - x(8) - x(5) = 0

      alin(5) = 2
      acol(5) = 1
      aval(5) =   1.0d0

      alin(6) = 2
      acol(6) = 8
      aval(6) = - 1.0d0

      alin(7) = 2
      acol(7) = 5
      aval(7) = - 1.0d0

      b(2) = 0.0d0

C     x(2) - x(9) - x(6) = 0

      alin(8) = 3
      acol(8) = 2
      aval(8) =   1.0d0

      alin(9) = 3
      acol(9) = 9
      aval(9) = - 1.0d0

      alin(10) = 3
      acol(10) = 6
      aval(10) = - 1.0d0

      b(3) = 0.0d0

      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.

C     COMMON SCALARS
      double precision c1,c2

C     COMMON BLOCKS
      common /probdata/ c1,c2

      flag = 0

      f = - 9.0d0 * x(1) - 15.0d0 * x(2) + 6.0d0 * x(3) + c1 * x(4)
     +    + 10.0d0 * ( x(5) + x(6) )

      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.

C     COMMON SCALARS
      double precision c1,c2

C     COMMON BLOCKS
      common /probdata/ c1,c2

      flag = 0

      g(1) = -  9.0d0
      g(2) = - 15.0d0
      g(3) =    6.0d0
      g(4) =       c1
      g(5) =   10.0d0
      g(6) =   10.0d0
      g(7) =    0.0d0
      g(8) =    0.0d0
      g(9) =    0.0d0

      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.
 
      flag = 0

      if ( ind .eq. 1 ) then
          c = x(7) * x(8) + 2.0d0 * x(5) - 2.5d0 * x(1)

      else if ( ind .eq. 2 ) then
          c = x(7) * x(9) + 2.0d0 * x(6) - 1.5d0 * x(2)

      else ! if ( ind .eq. 3 ) then
          c =  3.0d0 * x(3) + x(4) - x(7) * ( x(8) + x(9) )
      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evaljac(n,x,ind,indjac,valjac,jcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,n,jcnnz

C     ARRAY ARGUMENTS
      integer indjac(n)
      double precision x(n),valjac(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     indjac   integer indjac(jcnnz),
C              see below,
C
C     valjac   double precision valjac(jcnnz),
C              the non-null value of the partial derivative of the 
C              ind-th constraint with respect to the indjac(k)-th 
C              variable must be saved at valjac(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.

      flag = 0

      if ( ind .eq. 1 ) then

          jcnnz = 4

          indjac(1) = 7
          valjac(1) = x(8)

          indjac(2) = 8
          valjac(2) = x(7)

          indjac(3) = 5
          valjac(3) = 2.0d0

          indjac(4) = 1
          valjac(4) = - 2.5d0

      else if ( ind .eq. 2 ) then

          jcnnz = 4

          indjac(1) = 7
          valjac(1) = x(9)

          indjac(2) = 9
          valjac(2) = x(7)

          indjac(3) = 6
          valjac(3) = 2.0d0

          indjac(4) = 2
          valjac(4) = - 1.5d0

      else ! if ( ind .eq. 3 ) then

          jcnnz = 5

          indjac(1) = 3
          valjac(1) = 3.0d0

          indjac(2) = 4
          valjac(2) = 1.0d0

          indjac(3) = 7
          valjac(3) = - ( x(8) + x(9) )

          indjac(4) = 8
          valjac(4) = - x(7)

          indjac(5) = 9
          valjac(5) = - x(7) 

      end if

      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.

      flag = 0

      if ( ind .eq. 1 ) then
          hcnnz = 1

          hclin(1) = 8
          hccol(1) = 7
          hcval(1) = 1.0d0

      else if ( ind .eq. 2 ) then
          hcnnz = 1

          hclin(1) = 9
          hccol(1) = 7
          hcval(1) = 1.0d0

      else ! if ( ind .eq. 3 ) then
          hcnnz = 2

          hclin(1) = 8
          hccol(1) = 7
          hcval(1) = - 1.0d0

          hclin(2) = 9
          hccol(2) = 7
          hcval(2) = - 1.0d0
      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)

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.

      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)

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     ******************************************************************

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

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

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

C     COMMON SCALARS
      double precision c1,c2

C     COMMON BLOCKS
      common /probdata/ c1,c2

      mla   = ml
      annza = annz

C     x(7) * x(8) + 2.0d0 * x(5) - 2.5d0 * x(1) <= 0

C     Understimator of bilinear term:
C     l(7) * x(8) + l(8) * x(7) - l(8) * l(7) <= x(7) * x(8)

      ind(mla+1) = 1

      alin(annza+1) = mla + 1
      acol(annza+1) = 1
      aval(annza+1) = 2.5d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 5
      aval(annza+2) = - 2.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 7 
      aval(annza+3) = - l(8)

      alin(annza+4) = mla + 1
      acol(annza+4) = 8
      aval(annza+4) = - l(7)

      b(mla+1) = - l(8) * l(7)

      mla   = mla   + 1
      annza = annza + 4

C     Understimator of bilinear term:
C     u(7) * x(8) + u(8) * x(7) - u(8) * u(7) <= x(7) * x(8)

      ind(mla+1) = 1

      alin(annza+1) = mla + 1
      acol(annza+1) = 1
      aval(annza+1) = 2.5d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 5
      aval(annza+2) = - 2.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 7 
      aval(annza+3) = - u(8)

      alin(annza+4) = mla + 1
      acol(annza+4) = 8
      aval(annza+4) = - u(7)

      b(mla+1) = - u(8) * u(7)

      mla   = mla   + 1
      annza = annza + 4

C     x(7) * x(9) + 2.0d0 * x(6) - 1.5d0 * x(2) <= 0

C     Understimator of bilinear term:
C     l(7) * x(9) + l(9) * x(7) - l(9) * l(7) <= x(7) * x(9)

      ind(mla+1) = 2

      alin(annza+1) = mla + 1
      acol(annza+1) = 9
      aval(annza+1) = - l(7)

      alin(annza+2) = mla + 1
      acol(annza+2) = 6
      aval(annza+2) = - 2.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 2
      aval(annza+3) = 1.5d0

      alin(annza+4) = mla + 1
      acol(annza+4) = 7 
      aval(annza+4) = - l(9)

      b(mla+1) = - l(9) * l(7)

      mla   = mla   + 1
      annza = annza + 4

C     Understimator of bilinear term:
C     u(7) * x(9) + u(9) * x(7) - u(9) * u(7) <= x(7) * x(9)

      ind(mla+1) = 2

      alin(annza+1) = mla + 1
      acol(annza+1) = 9
      aval(annza+1) = - u(7)

      alin(annza+2) = mla + 1
      acol(annza+2) = 6
      aval(annza+2) = - 2.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 2
      aval(annza+3) = 1.5d0

      alin(annza+4) = mla + 1
      acol(annza+4) = 7 
      aval(annza+4) = - u(9)

      b(mla+1) = - u(9) * u(7)

      mla   = mla   + 1
      annza = annza + 4

C     3.0d0 * x(3) + x(4) - x(7) * x(8) - x(7) * x(9) = 0

C     3.0d0 * x(3) + x(4) - l(7) * x(8) - l(7) * x(9) - 
C     ( u(8) + u(9) ) * x(7) + u(8) * l(7) + u(9) * l(7) <= 0

      ind(mla+1) = 3

      alin(annza+1) = mla + 1
      acol(annza+1) = 3
      aval(annza+1) = - 3.0d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 4
      aval(annza+2) = - 1.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 8
      aval(annza+3) = l(7)

      alin(annza+4) = mla + 1
      acol(annza+4) = 9
      aval(annza+4) = l(7)

      alin(annza+5) = mla + 1
      acol(annza+5) = 7 
      aval(annza+5) = u(8) + u(9)

      b(mla+1) = u(8) * l(7) + u(9) * l(7)

      mla   = mla   + 1
      annza = annza + 5

C     3.0d0 * x(3) + x(4) - l(7) * x(8) - u(7) * x(9) - 
C     ( u(8) + l(9) ) * x(7) + u(8) * l(7) + l(9) * u(7) <= 0

      ind(mla+1) = 3

      alin(annza+1) = mla + 1
      acol(annza+1) = 3
      aval(annza+1) = - 3.0d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 4
      aval(annza+2) = - 1.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 8
      aval(annza+3) = l(7)

      alin(annza+4) = mla + 1
      acol(annza+4) = 9
      aval(annza+4) = u(7)

      alin(annza+5) = mla + 1
      acol(annza+5) = 7 
      aval(annza+5) = u(8) + l(9)

      b(mla+1) = u(8) * l(7) + l(9) * u(7)

      mla   = mla   + 1
      annza = annza + 5

C     3.0d0 * x(3) + x(4) - u(7) * x(8) - l(7) * x(9) - 
C     ( l(8) + u(9) ) * x(7) + l(8) * u(7) + u(9) * l(7) <= 0

      ind(mla+1) = 3

      alin(annza+1) = mla + 1
      acol(annza+1) = 3
      aval(annza+1) = - 3.0d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 4
      aval(annza+2) = - 1.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 8
      aval(annza+3) = u(7)

      alin(annza+4) = mla + 1
      acol(annza+4) = 9
      aval(annza+4) = l(7)

      alin(annza+5) = mla + 1
      acol(annza+5) = 7 
      aval(annza+5) = l(8) + u(9)

      b(mla+1) = l(8) * u(7) + u(9) * l(7)

      mla   = mla   + 1
      annza = annza + 5

C     3.0d0 * x(3) + x(4) - u(7) * x(8) - u(7) * x(9) - 
C     ( l(8) + l(9) ) * x(7) + l(8) * u(7) + l(9) * u(7) <= 0

      ind(mla+1) = 3

      alin(annza+1) = mla + 1
      acol(annza+1) = 3
      aval(annza+1) = - 3.0d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 4
      aval(annza+2) = - 1.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 8
      aval(annza+3) = u(7)

      alin(annza+4) = mla + 1
      acol(annza+4) = 9
      aval(annza+4) = u(7)

      alin(annza+5) = mla + 1
      acol(annza+5) = 7 
      aval(annza+5) = l(8) + l(9)

      b(mla+1) = l(8) * u(7) + l(9) * u(7)

      mla   = mla   + 1
      annza = annza + 5

C   - 3.0d0 * x(3) - x(4) + l(7) * x(8) + l(7) * x(9) + 
C     ( l(8) + l(9) ) * x(7) - l(8) * l(7) - l(9) * l(7) <= 0

      ind(mla+1) = - 3

      alin(annza+1) = mla + 1
      acol(annza+1) = 3
      aval(annza+1) = 3.0d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 4
      aval(annza+2) = 1.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 8
      aval(annza+3) = - l(7)

      alin(annza+4) = mla + 1
      acol(annza+4) = 9
      aval(annza+4) = - l(7)

      alin(annza+5) = mla + 1
      acol(annza+5) = 7 
      aval(annza+5) = - ( l(8) + l(9) )

      b(mla+1) = - ( l(8) * l(7) + l(9) * l(7) )

      mla   = mla   + 1
      annza = annza + 5

C   - 3.0d0 * x(3) - x(4) + l(7) * x(8) + u(7) * x(9) + 
C     ( l(8) + u(9) ) * x(7) - l(8) * l(7) - u(9) * u(7) <= 0

      ind(mla+1) = - 3

      alin(annza+1) = mla + 1
      acol(annza+1) = 3
      aval(annza+1) = 3.0d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 4
      aval(annza+2) = 1.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 8
      aval(annza+3) = - l(7)

      alin(annza+4) = mla + 1
      acol(annza+4) = 9
      aval(annza+4) = - u(7)

      alin(annza+5) = mla + 1
      acol(annza+5) = 7 
      aval(annza+5) = - ( l(8) + u(9) )

      b(mla+1) = - ( l(8) * l(7) + u(9) * u(7) )

      mla   = mla   + 1
      annza = annza + 5

C   - 3.0d0 * x(3) - x(4) + u(7) * x(8) + l(7) * x(9) + 
C     ( u(8) + l(9) ) * x(7) - u(8) * u(7) - l(9) * l(7) <= 0

      ind(mla+1) = - 3

      alin(annza+1) = mla + 1
      acol(annza+1) = 3
      aval(annza+1) = 3.0d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 4
      aval(annza+2) = 1.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 8
      aval(annza+3) = - u(7)

      alin(annza+4) = mla + 1
      acol(annza+4) = 9
      aval(annza+4) = - l(7)

      alin(annza+5) = mla + 1
      acol(annza+5) = 7 
      aval(annza+5) = - ( u(8) + l(9) )

      b(mla+1) = - ( u(8) * u(7) + l(9) * l(7) )

      mla   = mla   + 1
      annza = annza + 5

C   - 3.0d0 * x(3) - x(4) + u(7) * x(8) + u(7) * x(9) + 
C     ( u(8) + u(9) ) * x(7) - u(8) * u(7) - u(9) * u(7) <= 0

      ind(mla+1) = - 3

      alin(annza+1) = mla + 1
      acol(annza+1) = 3
      aval(annza+1) = 3.0d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 4
      aval(annza+2) = 1.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 8
      aval(annza+3) = - u(7)

      alin(annza+4) = mla + 1
      acol(annza+4) = 9
      aval(annza+4) = - u(7)

      alin(annza+5) = mla + 1
      acol(annza+5) = 7 
      aval(annza+5) = - ( u(8) + u(9) )

      b(mla+1) = - ( u(8) * u(7) + u(9) * u(7) )

      mla   = mla   + 1
      annza = annza + 5

C   - 9.0d0 * x(1) - 15.0d0 * x(2) + 6.0d0 * x(3) + c1 * x(4) + 
C     10.0d0 * ( x(5) + x(6) ) - f <= 0

      ind(mla+1) = 0

      alin(annza+1) = mla + 1
      acol(annza+1) = 1
      aval(annza+1) = 9.0d0

      alin(annza+2) = mla + 1
      acol(annza+2) = 2
      aval(annza+2) = 15.0d0

      alin(annza+3) = mla + 1
      acol(annza+3) = 3
      aval(annza+3) = - 6.0d0

      alin(annza+4) = mla + 1
      acol(annza+4) = 4
      aval(annza+4) = - c1

      alin(annza+5) = mla + 1
      acol(annza+5) = 5 
      aval(annza+5) = - 10.0d0

      alin(annza+6) = mla + 1
      acol(annza+6) = 6 
      aval(annza+6) = - 10.0d0

      b(mla+1) = - f

      mla   = mla   + 1
      annza = annza + 6

      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     This subroutine computes the intervalar objective function.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     y        double precision y(2*n),
C              bounds on x variables, y(2*i-1) <= x(i) <= y(2*i),
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.

C     COMMON SCALARS
      double precision c1,c2

C     COMMON BLOCKS
      common /probdata/ c1,c2

C     LOCAL ARRAYS
      double precision tmp(2)

C     EXTERNAL FUNCTIONS
      integer P

      flag = 0

C     f = - 9.0d0 * x(1) - 15.0d0 * x(2) + 6.0d0 * x(3) + c1 * x(4)
C    +    + 10.0d0 * ( x(5) + x(6) )

      call sclmlt(-9.0d0,y(P(1)),f)

      call sclmlt(-15.0d0,y(P(2)),tmp)
      call add(f,tmp,f)

      call sclmlt(6.0d0,y(P(3)),tmp)
      call add(f,tmp,f)

      call sclmlt(c1,y(P(4)),tmp)
      call add(f,tmp,f)

      call add(y(P(5)),y(P(6)),tmp)
      call sclmlt(10.0d0,tmp,tmp)
      call add(f,tmp,f)

      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     This subroutine computes the intervalar Hessian matrix of the 
C     objective function.
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     y        double precision y(2*n),
C              bounds on x variables, y(2*i-1) <= x(i) <= y(2*i),
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 interval of the non-null value of the (hlin(k),hcol(k)) 
C              position of the Hessian matrix of the objective function must 
C              be saved at hval(P(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 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     This subroutine computes the interval of the ind-th constraint of 
C     your problem. 
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     y        double precision y(2*n),
C              bounds on x variables, y(2*i-1) <= x(i) <= y(2*i),
C
C     ind      integer,
C              index of the constraint to be computed,
C
C     On Return
C
C     c        double precision,
C              interval of the ind-th constraint "evaluated at y",
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 ARRAYS
      double precision tmp(2)
 
C     EXTERNAL FUNCTION
      integer P

      flag = 0

      if ( ind .eq. 1 ) then
C         c = x(7) * x(8) + 2.0d0 * x(5) - 2.5d0 * x(1)
          call ivli(y(P(7)),c)
          call mult(c,y(P(8)),c)
          call sclmlt(2.0d0,y(P(5)),tmp)
          call add(c,tmp,c)
          call sclmlt(-2.5d0,y(P(1)),tmp)
          call add(c,tmp,c)

      else if ( ind .eq. 2 ) then
C         c = x(7) * x(9) + 2.0d0 * x(6) - 1.5d0 * x(2)
          call ivli(y(P(7)),c)
          call mult(c,y(P(9)),c)
          call sclmlt(2.0d0,y(P(6)),tmp)
          call add(c,tmp,c)
          call sclmlt(-1.5d0,y(P(2)),tmp)
          call add(c,tmp,c)

      else ! if ( ind .eq. 3 ) then
C         c = 3.0d0 * x(3) + x(4) - x(7) * ( x(8) + x(9) )
          call sclmlt(3.0d0,y(P(3)),c)
          call add(c,y(P(4)),c)
          call add(y(P(8)),y(P(9)),tmp)
          call mult(tmp,y(P(7)),tmp)
          call sub(c,tmp,c)
      end if

      end

C     ******************************************************************
C     ******************************************************************
 
      subroutine evaljaci(n,y,ind,indjac,valjac,jcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS
      integer flag,ind,n,jcnnz

C     ARRAY ARGUMENTS
      integer indjac(n)
      double precision y(2*n),valjac(2*n)

C     This subroutine must compute the intervalar gradient of the 
C     constraint ind. 
C
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     y        double precision y(2*n),
C              bounds on x variables, y(2*i-1) <= x(i) <= y(2*i),
C
C     ind      integer,
C              index of the constraint whose intervalar gradient will be 
C              computed,
C
C     On Return
C
C     jcnnz   integer,
C              number of perhaps-non-null elements of the computed 
C              gradient,
C
C     indjac   integer indjac(jcnnz),
C              see below,
C
C     valjac   double precision valjac(jcnnz),
C              the interval of the non-null value of the partial 
C              derivative of the ind-th constraint with respect to the 
C              indjac(k)-th variable must be saved at valjac(P(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     EXTERNAL FUNCTION
      integer P

      flag = 0

      if ( ind .eq. 1 ) then

          jcnnz = 4

          indjac(1) = 7
C         valjac(1) = x(8)
          call ivli(y(P(8)),valjac(P(1)))

          indjac(2) = 8
C         valjac(2) = x(7)
          call ivli(y(P(7)),valjac(P(2)))

          indjac(3) = 5
C         valjac(3) = 2.0d0
          call ivl1(2.0d0,valjac(P(3)))

          indjac(4) = 1
C         valjac(4) = - 2.5d0
          call ivl1(-2.5d0,valjac(P(4)))

      else if ( ind .eq. 2 ) then

          jcnnz = 4

          indjac(1) = 7
C         valjac(1) = x(9)
          call ivli(y(P(9)),valjac(P(1)))

          indjac(2) = 9
C         valjac(2) = x(7)
          call ivli(y(P(7)),valjac(P(2)))

          indjac(3) = 6
C         valjac(3) = 2.0d0
          call ivl1(2.0d0,valjac(P(3)))

          indjac(4) = 2
C         valjac(4) = - 1.5d0
          call ivl1(-1.5d0,valjac(P(4)))

      else ! if ( ind .eq. 3 ) then

          jcnnz = 5

          indjac(1) = 3
C         valjac(1) = 3.0d0
          call ivl1(3.0d0,valjac(P(1)))

          indjac(2) = 4
C         valjac(2) = 1.0d0
          call ivl1(1.0d0,valjac(P(2)))

          indjac(3) = 7
C         valjac(3) = - ( x(8) + x(9) )
          call add(y(P(8)),y(P(9)),valjac(P(3)))
          call sclmlt(-1.0d0,valjac(P(3)),valjac(P(3)))

          indjac(4) = 8
C         valjac(4) = - x(7)
          call sclmlt(-1.0d0,y(P(7)),valjac(P(4)))

          indjac(5) = 9
C         valjac(5) = - x(7) 
          call sclmlt(-1.0d0,y(P(7)),valjac(P(5)))

      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     This subroutine might compute the intervalar Hessian matrix of 
C     the ind-th constraint. 
C     
C     Parameters of the subroutine:
C
C     On Entry:
C
C     n        integer,
C              number of variables,
C
C     y        double precision y(2*n),
C              bounds on x variables, y(2*i-1) <= x(i) <= y(2*i),
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 interval of the non-null value of the (hclin(k),hccol(k)) 
C              position of the Hessian matrix of the ind-th constraint must 
C              be saved at hcval(P(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     EXTERNAL FUNCTION
      integer P

      flag = 0

      if ( ind .eq. 1 ) then
          hcnnz = 1

          hclin(1) = 8
          hccol(1) = 7
C         hcval(1) = 1.0d0
          call ivl1(1.0d0,hcval(P(1)))

      else if ( ind .eq. 2 ) then
          hcnnz = 1

          hclin(1) = 9
          hccol(1) = 7
C         hcval(1) = 1.0d0
          call ivl1(1.0d0,hcval(P(1)))

      else ! if ( ind .eq. 3 ) then
          hcnnz = 2

          hclin(1) = 8
          hccol(1) = 7
C         hcval(1) = - 1.0d0
          call ivl1(-1.0d0,hcval(P(1)))

          hclin(2) = 9
          hccol(2) = 7
C         hcval(2) = - 1.0d0
          call ivl1(-1.0d0,hcval(P(2)))
      end if

      end