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

      subroutine testi(n,m,l,u,inform)

C     SCALAR ARGUMENTS
      integer inform,m,n

C     ARRAY ARGUMENTS
      double precision l(n),u(n)

C     This subroutine tests the intervalar subroutines evalf, evalgi, 
C     evalhi, evalci, evaljaci and evalhci.

C     PARAMETERS
      include "dim.f"

C     LOCAL SCALARS
      logical ok
      integer hnnz,hnnzint,i,ind,j,jcnnz,jcnnzint,k,kmax,flag
      double precision c,f,lt,seed,ut

C     LOCAL ARRAYS
      integer hcol(hnnzmax),hcolint(hnnzmax),hlin(hnnzmax),
     +        hlinint(hnnzmax),indjac(nmax),indjacint(nmax)
      double precision cint(2),fint(2),hval(hnnzmax),hvalint(2*hnnzmax),
     +        valjac(nmax),valjacint(2*nmax),x(nmax),y(2*nmax)

C     EXTERNAL FUNCTION
      integer P
      double precision drand,iinf,isup

      inform = 0

      kmax = 1000

      do i = 1,n
          y(2*i-1) = l(i)
          y(2*i  ) = u(i)
      end do

C     ------------------------------------------------------------------
C     TESTING INTERVALAR OBJECTIVE FUNCTION
C     ------------------------------------------------------------------

      call evalfi(n,y,fint,flag)
      if ( flag .ne. 0 ) then
          inform = - 80
          return
      end if

C     write(*,*) iinf(fint),isup(fint)

      ok = .true.

      k = 1
 5    if ( k .le. kmax .and. ok ) then

C         Set a random point in the box
          seed = k

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

C         Compute objective function
          call evalf(n,x,f,flag)
          if ( flag .ne. 0 ) then
              inform = - 90
              return
          end if

C         Test wether it is inside its interval
          lt = iinf(fint)
          ut = isup(fint)
          if ( f .lt. lt .or. f .gt. ut ) then
              ok = .false.
              write(*,*) 'Wrong interval for objective function'
              write(*,*) 'Interval: ',lt,ut,' value: ',f
          end if

          k = k + 1
          go to 5
      end if

      if ( ok ) then
          write(*,*) 'Intervalar objective function is ok'
      else
          write(*,*) 'Point: ',(x(j),j=1,n)
          write(*,*) 'Lower bounds: ',(l(i),i=1,n)
          write(*,*) 'Upper bounds: ',(u(i),i=1,n)
          inform = - 99
          return
      end if

C     ------------------------------------------------------------------
C     TESTING INTERVALAR HESSIAN OF THE OBJECTIVE FUNCTION
C     ------------------------------------------------------------------

      call evalhi(n,y,hlinint,hcolint,hvalint,hnnzint,flag)
      if ( flag .ne. 0 ) then
          inform = - 84
          return
      end if

      ok = .true.

      k = 1
 20   if ( k .le. kmax .and. ok ) then

C         Set a random point in the box
          seed = k

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

C         Compute the Hessian
          call evalh(n,x,hlin,hcol,hval,hnnz,flag)
          if ( flag .ne. 0 ) then
              inform = - 94
              return
          end if

          if ( hnnz .ne. hnnzint ) then

              ok = .false.
              write(*,*) 'Different number of elements for Hessian'
              write(*,*) 'hnnz: ',hnnz,' hnnzint: ',hnnzint

          else
C             Test wether each component is inside its interval
              do i = 1,hnnz
                  if ( hlin(i) .ne. hlinint(i) .or.
     +               hcol(i) .ne. hcolint(i) ) then
                    write(*,*) 'Different non-null elements'
                    write(*,*) 'i: ',i,' hlin(i): ',hlin(i),
     +              ' hcol(i): ',hcol(i),' hlinint(i): ',hlinint(i),
     +              ' hcolint(i): ',hcolint(i)
                  else
                    lt = iinf(hvalint(P(i)))
                    ut = isup(hvalint(P(i)))
                    if ( hval(i) .lt. lt .or. hval(i) .gt. ut ) then
                        ok = .false.
                        write(*,*) 'Wrong interval for Hessian ',
     +                  'component ',i,hlin(i),hcol(i)
                        write(*,*) 'Interval: ',lt,ut,' value: ',hval(i)
                    end if
                 end if
              end do
          end if

          k = k + 1
          go to 20
      end if

      if ( ok ) then
          write(*,*) 'Intervalar Hessian is ok'
      else
          write(*,*) 'Point: ',(x(j),j=1,n)
          inform = - 99
          return
      end if

C     ------------------------------------------------------------------
C     TESTING INTERVALAR CONSTRAINTS
C     ------------------------------------------------------------------

      do ind = 1,m

          call evalci(n,y,ind,cint,flag)
          if ( flag .ne. 0 ) then
              inform = - 81
              return
          end if

C         write(*,*) iinf(cint),isup(cint)

          ok = .true.

          k = 1
 30       if ( k .le. kmax .and. ok ) then

C             Set a random point in the box
              seed = k

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

C             Compute ind-th constraint
              call evalc(n,x,ind,c,flag)
              if ( flag .ne. 0 ) then
                  inform = - 91
                  return
              end if

C             Test wether it is inside its interval
              lt = iinf(cint)
              ut = isup(cint)
              if ( c .lt. lt .or. c .gt. ut ) then
                  ok = .false.
                  write(*,*) 'Wrong interval for constraint ',ind
                  write(*,*) 'Interval: ',lt,ut,' value: ',c
              end if

              k = k + 1
              go to 30
          end if

          if ( ok ) then
              write(*,*) 'Intervalar constraint ',ind,' is ok'
          else
              write(*,*) 'Point: ',(x(j),j=1,n)
              inform = - 99
              return
          end if

      end do

C     ------------------------------------------------------------------
C     TESTING INTERVALAR GRADIENT OF CONSTRAINTS
C     ------------------------------------------------------------------

      do ind = 1,m

          call evaljaci(n,y,ind,indjacint,valjacint,jcnnzint,flag)
          if ( flag .ne. 0 ) then
              inform = - 83
              return
          end if

C         do i = 1,jcnnzint
C             write(*,*) indjacint(i),iinf(valjacint(P(i))),
C    +                                isup(valjacint(P(i)))
C         end do

          ok = .true.

          k = 1
 40       if ( k .le. kmax .and. ok ) then

C             Set a random point in the box
              seed = k

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

C             Compute the gradient
              call evaljac(n,x,ind,indjac,valjac,jcnnz,flag)
              if ( flag .ne. 0 ) then
                  inform = - 93
                  return
              end if

              if ( jcnnz .ne. jcnnzint ) then

                  ok = .false.
                  write(*,*) 'Different number of elements in ',
     +            'gradient of constraint ',ind
                  write(*,*) 'jcnnz: ',jcnnz,' jcnnzint: ',jcnnzint

              else
C                 Test wether each component is inside its interval
                  do i = 1,jcnnz
                      if ( indjac(i) .ne. indjacint(i) ) then
                        write(*,*) 'Different non-null elements in ',
     +                  'of constraint ',ind
                        write(*,*) 'i: ',i,' indjac(i): ',indjac(i),
     +                  ' indjacint(i): ',indjacint(i)
                      else
                        lt = iinf(valjacint(P(i)))
                        ut = isup(valjacint(P(i)))
                        if ( valjac(i) .lt. lt .or. 
     +                       valjac(i) .gt. ut ) then
                            ok = .false.
                            write(*,*) 'Wrong interval for gradient ',
     +                      'component ',i,indjac(i),' of ',
     +                      'constraint ',ind
                            write(*,*) 'Int: ',lt,ut,' val: ',valjac(i)
                        end if
                     end if
                  end do
              end if

              k = k + 1
              go to 40
          end if

          if ( ok ) then
              write(*,*) 'Intervalar gradient of constraint ',ind,
     +        ' is ok'
          else
              write(*,*) 'Point: ',(x(j),j=1,n)
              inform = - 99
              return
          end if

      end do

C     ------------------------------------------------------------------
C     TESTING INTERVALAR HESSIAN OF CONSTRAINTS
C     ------------------------------------------------------------------

      do ind = 1,m

          call evalhci(n,y,ind,hlinint,hcolint,hvalint,hnnzint,flag)
          if ( flag .ne. 0 ) then
              inform = - 85
              return
          end if

          ok = .true.

          k = 1
 50       if ( k .le. kmax .and. ok ) then

C             Set a random point in the box
              seed = k

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

C             Compute the Hessian
              call evalhc(n,x,ind,hlin,hcol,hval,hnnz,flag)
              if ( flag .ne. 0 ) then
                  inform = - 95
                  return
              end if

              if ( hnnz .ne. hnnzint ) then

                  ok = .false.
                  write(*,*) 'Different number of elements in Hessian ',
     +            'of constraint ',ind
                  write(*,*) 'hnnz: ',hnnz,' hnnzint: ',hnnzint

              else
C                 Test wether each component is inside its interval
                  do i = 1,hnnz
                      if ( hlin(i) .ne. hlinint(i) .or.
     +                   hcol(i) .ne. hcolint(i) ) then
                        write(*,*) 'Different non-null elements in ',
     +                  'Hessian of constraint ',ind
                        write(*,*) 'i: ',i,' hlin(i): ',hlin(i),
     +                  ' hcol(i): ',hcol(i),' hlinint(i): ',
     +                  hlinint(i),' hcolint(i): ',hcolint(i)
                      else
                        lt = iinf(hvalint(P(i)))
                        ut = isup(hvalint(P(i)))
                        if ( hval(i) .lt. lt .or. hval(i) .gt. ut ) then
                            ok = .false.
                            write(*,*) 'Wrong interval for Hessian ',
     +                      'component ',i,hlin(i),hcol(i),' of ',
     +                      'constraint ',ind
                            write(*,*) 'Interval: ',lt,ut,' value: ',
     +                      hval(i)
                          end if
                     end if
                  end do
              end if

              k = k + 1
              go to 50
          end if

          if ( ok ) then
              write(*,*) 'Intervalar Hessian of constraint ',ind,
     +        ' is ok'
          else
              write(*,*) 'Point: ',(x(j),j=1,n)
              inform = - 99
              return
          end if

      end do

      end 

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

      subroutine testali(n,l,u,m,lambda,rho,equatn,linear,inform)

C     SCALAR ARGUMENTS
      integer inform,m,n

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

C     This subroutine tests evalali and evalhali.

C     PARAMETERS
      include "dim.f"

C     LOCAL SCALARS
      logical ok
      integer i,j,k,kmax
      double precision al,lt,seed,ut

C     LOCAL ARRAYS
      double precision alint(2),H(nmax,nmax),Hint(2*nmax**2),x(nmax),
     +        y(2*nmax)

C     EXTERNAL FUNCTION
      integer P2
      double precision drand,iinf,isup

      kmax = 1000

      do i = 1,n
          y(2*i-1) = l(i)
          y(2*i  ) = u(i)
      end do

C     ------------------------------------------------------------------
C     TESTING INTERVALAR AUGMENTED LAGRANGIAN
C     ------------------------------------------------------------------

      call evalali(n,y,m,lambda,rho,equatn,linear,alint,inform)
      if ( inform .lt. 0 ) then
          return
      end if

C     write(*,*) iinf(alint),isup(alint)

      ok = .true.

      k = 1
 5    if ( k .le. kmax .and. ok ) then

C         Set a random point in the box
          seed = k

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

C         Compute the Augmented Lagrangian
          call evalal(n,x,m,lambda,rho,equatn,linear,al,inform)
          if ( inform .lt. 0 ) then
              return
          end if

           lt = iinf(alint)
           ut = isup(alint)
           if ( al .lt. lt .or. al .gt. ut ) then
               ok = .false.
               write(*,*) 'Wrong interval Augmented Lagrangian'
               write(*,*) 'Interval: ',lt,ut,' value: ',al
          end if

          k = k + 1
          go to 5
      end if

      if ( ok ) then
          write(*,*) 'Intervalar Augmented Lagrangian function is ok'
      else
          write(*,*) 'Point: ',(x(j),j=1,n)
          inform = - 99
          return
      end if

C     ------------------------------------------------------------------
C     TESTING INTERVALAR HESSIAN OF THE AUGMENTED LAGRANGIAN
C     ------------------------------------------------------------------

      call evalhali(n,y,m,lambda,rho,equatn,linear,Hint,inform)
      if ( inform .lt. 0 ) then
          return
      end if

      ok = .true.

      k = 1
 10   if ( k .le. kmax .and. ok ) then

C         Set a random point in the box
          seed = k

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

C         Compute the Hessian
          call evalhal(n,x,m,lambda,rho,equatn,linear,nmax,H,inform)
          if ( inform .lt. 0 ) then
              return
          end if

          do i = 1,n
              do j = 1,i
                  lt = iinf(Hint(P2(i,j,n)))
                  ut = isup(Hint(P2(i,j,n)))
                  if ( H(i,j) .lt. lt .or. H(i,j) .gt. ut ) then
                      ok = .false.
                      write(*,*) 'Wrong interval for AL Hessian ',
     +                'component ',i,j
                      write(*,*) 'Int: ',lt,ut,' value: ',H(i,j)
                 end if
              end do
          end do

          k = k + 1
          go to 10
      end if

      if ( ok ) then
          write(*,*) 'Intervalar Hessian of AL is ok'
      else
          write(*,*) 'Point: ',(x(j),j=1,n)
          inform = - 99
          return
      end if

      end

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

      subroutine testun(n,l,u,m,lambda,rho,equatn,linear,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer inform,m,n

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

C     This subroutine tests the convexity of the underestimator of 
C     the Augmented Lagrangian function.

C     PARAMETERS
      include "dim.f"

C     LOCAL SCALARS
      logical ok
      integer j,k,kmax
      double precision seed,tmp,un1,un2,un3

C     LOCAL ARRAYS
      double precision x1(nmax),x2(nmax),x3(nmax)

C     EXTERNAL FUNCTIONS
      double precision drand

      kmax = 1000

C     ------------------------------------------------------------------
C     TESTING CONVEXITY
C     ------------------------------------------------------------------

      ok = .true.

      k = 1
 10   if ( k .le. kmax .and. ok ) then

          seed = k

C         Set two random points in the box

          do j = 1,n
              x1(j) = l(j) + ( u(j) - l(j) ) * drand(seed)
          end do

          do j = 1,n
              x2(j) = l(j) + ( u(j) - l(j) ) * drand(seed)
          end do

C         Set a random point in the segment joining x1 and x2

          tmp = drand(seed)
          do j = 1,n
              x3(j) = ( 1.0d0 - tmp ) * x1(j) + tmp * x2(j)
          end do

C         Compute un in x1, x2 and x3

          call evalun(n,x1,m,lambda,rho,equatn,linear,un1,inform)
          if ( inform .lt. 0 ) then
              return
          end if

          call evalun(n,x2,m,lambda,rho,equatn,linear,un2,inform)
          if ( inform .lt. 0 ) then
              return
          end if

          call evalun(n,x3,m,lambda,rho,equatn,linear,un3,inform)
          if ( inform .lt. 0 ) then
              return
          end if

C         Verify convexity
          if ( un3 .gt. ( 1.0d0 - tmp ) * un1 + tmp * un2 + 1.0d-08 ) 
     +    then
              ok = .false.
          end if

          k = k + 1
          go to 10
      end if

      if ( ok ) then
          write(*,*) 'The underestimator is convex!'
      else
          write(*,*) 'Underestimator failure!',un1,un2,un3,tmp,
     +    un1 + ( un2 - un1 ) * tmp
          inform = - 99
          return
      end if

      end

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

      subroutine testaddc(n,l,u,m,lambda,rho,equatn,linear,ml,meq,alin,
     +acol,aval,annz,b,mla,annza,ind,fbest,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer annz,annza,inform,m,meq,ml,mla,n
      double precision fbest

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      integer alin(annz),acol(annz),ind(ml)
      double precision aval(annz),b(ml),l(n),lambda(m),rho(m),u(n)

C     This subroutine tests that the addicional constraints are valid
C     constraints and that they are not reducing the feasible set.

C     PARAMETERS
      include "dim.f"

C     LOCAL SCALARS
      logical ok
      integer i,j,k,kmax,flag
      double precision al,atx,c,seed

C     LOCAL ARRAYS
      double precision x(nmax)

C     EXTERNAL FUNCTIONS
      double precision drand

      inform = 0

      kmax = 1000

C     ------------------------------------------------------------------
C     TESTING ADITIONAL CONSTRAINTS
C     ------------------------------------------------------------------

      ok = .true.

      k = 1
 10   if ( k .le. kmax .and. ok ) then

          seed = k

C         Set a random points in the box

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

C         Test the aditional constraints at that point

          do j = ml + 1,mla

C             Compute the aditional constraint
              atx = b(j)
              do i = annz + 1,annza
                  if ( alin(i) .eq. j ) then
                      atx = atx - aval(i) * x(acol(i))
                  end if
              end do

              if ( ind(j) .ne. 0 ) then

C                 Compute the original constraint of which the aditional
C                 one is an underestimator
                  call evalc(n,x,abs(ind(j)),c,flag)
                  if ( flag .ne. 0 ) then
                      inform = - 91
                      return
                  end if

                  if ( ind(j) .lt. 0 ) then
                      c = - c
                  end if

C                 Compare them
                  if ( atx .gt. c .and. abs(atx - c) .gt. 1.0d-08 ) then
                      write(*,*) 'Aditional constraint ',j-ml,' is not',
     +                ' an underestimator of constraint ',abs(ind(j)),
     +                atx,c
                      ok = .false.
                  end if

              else

C                 Compute the Augmented Lagrangian
                  call evalal(n,x,m,lambda,rho,equatn,linear,al,inform)
                  if ( inform .lt. 0 ) then
                      return
                  end if

C                 Compare them
                  if ( atx .gt. al - fbest .and. 
     +                 abs(atx - ( al - fbest ) ) .gt. 1.0d-08 ) then
                      write(*,*) 'Aditional constraint ',j-ml,' is ',
     +                ' not an underestimator of the AL function',atx,
     +                al - fbest
                      ok = .false.
                  end if
                  
              end if

          end do

          k = k + 1
          go to 10
      end if

      if ( ok ) then
          write(*,*) 'The aditional constraints are valid!'
      else
          write(*,*) 'Aditional constraint failure!'
          write(*,*) 'Point      : ',(x(j),j=1,n)
          write(*,*) 'Lower bound: ',(l(j),j=1,n)
          write(*,*) 'Upper bound: ',(u(j),j=1,n)
          inform = - 99
          return
      end if

      end