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

      subroutine evalun(n,x,m,lambda,rho,equatn,linear,un,inform)

      implicit none

C     SCALAR ARGUMENTS
      double precision un
      integer inform,m,n

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

C     This subroutine computes the underestimator of the Augmented 
C     Lagrangian function.

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      logical under
      character * 5 utype

C     COMMON ARRAYS
      integer nsub(nmax)
      double precision alpha(nmax,nsubmax),beta(nmax,nsubmax),
     +        gamma(nmax,nsubmax),y(2*nmax)

C     LOCAL SCALARS
      integer i,j
      double precision l,q,tmp,u,yl,ys,yu

C     COMMON BLOCKS
      common /undata/ alpha,beta,gamma,y,nsub,under,utype

C     EXTERNAL FUNCTIONS
      integer P
      double precision iinf,isup

      call evalal(n,x,m,lambda,rho,equatn,linear,un,inform)
      if ( inform .lt. 0 ) then
          return
      end if

      if ( under ) then

C         PLAIN alpha-BB
          if ( utype .eq. 'PLAIN' ) then

              q = 0.0d0
              do i = 1,n
                  l = iinf(y(P(i)))
                  u = isup(y(P(i)))
                  tmp = max( 0.0d0, u - x(i) ) * max( 0.0d0, x(i) - l )
                  q = q + tmp * alpha(i,1)
              end do
              un = un - q

C         PIECEWISE alphaBB 
          else ! if ( utype .eq. 'PWISE' ) then

              q = 0.0d0
              do i = 1,n
                  yl = iinf(y(P(i)))
                  yu = isup(y(P(i)))
                  if ( nsub(i) .eq. 1 ) then
                      j  = 1
                      l  = yl
                      u  = yu
                  else
                      ys = ( yu - yl ) / nsub(i)
                      j  = min( int( ( x(i) - yl ) / ys ) + 1, nsub(i) )
                      l  = yl + ( j - 1 ) * ys
                      u  = yl +   j       * ys
                  end if
                  tmp = max( 0.0d0, u - x(i) ) * max( 0.0d0, x(i) - l )
                  q = 
     +            q + tmp * alpha(i,j) + x(i) * beta(i,j) + gamma(i,j) 
              end do
              un = un - q

          end if

      end if

      end

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

      subroutine evalnun(n,x,m,lambda,rho,equatn,linear,nun,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer inform,m,n

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

C     This subroutine computes the gradient of the underestimator of 
C     the Augmented Lagrangian function.

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      logical under
      character * 5 utype

C     COMMON ARRAYS
      integer nsub(nmax)
      double precision alpha(nmax,nsubmax),beta(nmax,nsubmax),
     +        gamma(nmax,nsubmax),y(2*nmax)

C     LOCAL SCALARS
      integer i,j
      double precision l,tmp,u,yl,ys,yu

C     COMMON BLOCKS
      common /undata/ alpha,beta,gamma,y,nsub,under,utype

C     EXTERNAL FUNCTIONS
      integer P
      double precision iinf,isup

      call evalnal(n,x,m,lambda,rho,equatn,linear,nun,inform)
      if ( inform .lt. 0 ) then
          return
      end if

      if ( under ) then

C         PLAIN alpha-BB
          if ( utype .eq. 'PLAIN' ) then

              do i = 1,n
                  l = iinf(y(P(i)))
                  u = isup(y(P(i)))
                  tmp = u + l - 2.0d0 * x(i)
                  nun(i) = nun(i) - tmp * alpha(i,1)
              end do

C         PIECEWISE alphaBB 
          else ! if ( utype .eq. 'PWISE' ) then

              do i = 1,n
                  yl = iinf(y(P(i)))
                  yu = isup(y(P(i)))
                  if ( nsub(i) .eq. 1 ) then
                      j  = 1
                      l  = yl
                      u  = yu
                  else
                      ys = ( yu - yl ) / nsub(i)
                      j  = min( int( ( x(i) - yl ) / ys ) + 1, nsub(i) )
                      l  = yl + ( j - 1 ) * ys
                      u  = yl +   j       * ys
                  end if
                  tmp = u + l - 2.0d0 * x(i)
                  nun(i) = nun(i) - tmp * alpha(i,j) - beta(i,j)
              end do

          end if

      end if

      end

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

      subroutine evalhun(n,x,m,lambda,rho,equatn,linear,ldh,H,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer inform,ldh,m,n

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

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      logical under
      character * 5 utype

C     COMMON ARRAYS
      integer nsub(nmax)
      double precision alpha(nmax,nsubmax),beta(nmax,nsubmax),
     +        gamma(nmax,nsubmax),y(2*nmax)

C     LOCAL SCALARS
      integer i,j
      double precision yl,ys,yu

C     COMMON BLOCKS
      common /undata/ alpha,beta,gamma,y,nsub,under,utype

C     EXTERNAL FUNCTIONS
      integer P
      double precision iinf,isup

      call evalhal(n,x,m,lambda,rho,equatn,linear,ldh,H,inform)
      if ( inform .lt. 0 ) then
          return
      end if

      if ( under ) then

C         PLAIN alpha-BB
          if ( utype .eq. 'PLAIN' ) then

              do i = 1,n
                  H(i,i) = H(i,i) + 2.0d0 * alpha(i,1)
              end do

C         PIECEWISE alphaBB 
          else ! if ( utype .eq. 'PWISE' ) then

              do i = 1,n
                  if ( nsub(i) .eq. 1 ) then
                      j  = 1
                  else
                      yl = iinf(y(P(i)))
                      yu = isup(y(P(i)))
                      ys = ( yu - yl ) / nsub(i)
                      j  = min( int( ( x(i) - yl ) / ys ) + 1, nsub(i) )
                  end if
                  H(i,i) = H(i,i) + 2.0d0 * alpha(i,j)
              end do

          end if

      end if

      end