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

      subroutine evalali(n,y,m,lambda,rho,equatn,linear,al,inform)

      implicit none

C     SCALAR ARGUMENTS
      double precision al(2)
      integer inform,m,n

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      double precision lambda(m),rho(m),y(2*n)

C     This subroutine computes the Augmented Lagrangian function.

C     LOCAL SCALARS
      integer flag,i
      double precision c(2),f(2),p(2)

      inform = 0

C     COMPUTE OBJECTIVE FUNTION

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

C     COMPUTES AL(x) = f(x) + \sum_j=1^m P(c(j),rho(j),lambda(j))

C     al = f
      call ivli(f,al)

      do i = 1,m

C         COMPUTE i-TH CONSTRAINT
          call evalci(n,y,i,c,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

C         ADD P(c(i),rho(i),lambda(i))
          call evalpi(c,rho(i),lambda(i),equatn(i),p)

C         al = al + p
          call add(al,p,al)

      end do

      end

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

      subroutine evalhali(n,y,m,lambda,rho,equatn,linear,H,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer inform,m,n

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      double precision H(2*n**2),lambda(m),rho(m),y(m)

C     This subroutine cumputes the interval Hessian of the Augmented
C     Lagrangian.

C     PARAMETERS
      include "dim.f"

C     LOCAL SCALARS
      logical addterm
      integer flag,hnnz,i,i1,i2,ind,ind1,ind2,j,jcnnz

C     LOCAL ARRAYS
      integer hcol(hnnzmax),hlin(hnnzmax),indjac(jcnnzmax)
      double precision c(2),dpdc(2),hval(hnnzmax),tmp(2),
     +        valjac(jcnnzmax),zeroi(2)

C     EXTERNAL FUNCTIONS
      integer P,P2
      double precision iinf,isup

      inform = 0

C     SET CONSTANT INTERVALS

      call ivl1(0.0d0,zeroi)

C     INITIALIZE DENSE AUGMENTED LAGRANGIAN HESSIAN

      do j = 1,n
          do i = j,n
C             H(i,j) = 0.0d0
              ind = P2(i,j,n)
              call ivl1(0.0d0,H(ind))
          end do
      end do

C     COMPUTE THE INTERVALAR HESSIAN OF THE OBJECTIVE FUNCTION

      call evalhi(n,y,hlin,hcol,hval,hnnz,flag)
      if ( flag .ne. 0 ) then
          inform = - 84
          return
      end if

      do i = 1,hnnz
C         Ignore upper triangle elements
          if ( hlin(i) .ge. hcol(i) ) then
              ind = P2(hlin(i),hcol(i),n)
              call add(H(ind),hval(P(i)),H(ind))
          end if
      end do

C     COMPUTE THE INTERVALAR GRADIENTS AND HESSIANS OF THE CONSTRAINTS

      do j = 1,m
C     if ( .not. linear(j) ) then

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

C         dpdc = lambda + rho * c
          call sclmlt(rho(j),c,dpdc)
          call scladd(lambda(j),dpdc,dpdc)

C         FIRST ORDER TERMS
 
          call evaljaci(n,y,j,indjac,valjac,jcnnz,flag)
          if ( flag .ne. 0 ) then
              inform = - 83
              return
          end if

          do i1 = 1,jcnnz
              ind1 = indjac(i1)
              do i2 = 1,jcnnz
                  ind2 = indjac(i2)
                  if ( ind1 .ge. ind2 ) then
                      call mult(valjac(P(i1)),valjac(P(i2)),tmp)
                      call sclmlt(rho(j),tmp,tmp)
                      if ( equatn(j) ) then
C                         Add
                          addterm = .true.
                      else
                          if ( iinf(dpdc) .gt. 0.0d0 ) then
C                             Add
                              addterm = .true.
                          else if ( isup(dpdc) .le. 0.0d0 ) then
C                             Do not add
                              addterm = .false.
                          else 
C                             Add but also include zero into the interval
                              call ihull(tmp,zeroi,tmp)
                              addterm = .true.
                          end if
                      end if
                      if ( addterm ) then
                          ind = P2(ind1,ind2,n)
                          call add(H(ind),tmp,H(ind))
                      end if
                  end if
              end do
          end do

C         SECOND ORDER TERMS

          call evalhci(n,y,j,hlin,hcol,hval,hnnz,flag)
          if ( flag .ne. 0 ) then
              inform = - 85
              return
          end if

          do i = 1,hnnz
              ind1 = hlin(i)
              ind2 = hcol(i)
              if ( ind1 .ge. ind2 ) then
                  if ( equatn(j) ) then
C                     Add multiplied by dpdc
                      call ivli(dpdc,tmp)
                      addterm = .true.
                  else
                      if ( iinf(dpdc) .gt. 0.0d0 ) then
C                         Add multiplied by dpdc
                          call ivli(dpdc,tmp)
                          addterm = .true.
                      else if ( isup(dpdc) .le. 0.0d0 ) then
C                         Do not add because max(0, dpdc) = 0
                          addterm = .false.
                      else 
C                         Add but multiplied by max(0,dpdc)
                          call ivl2(0.0d0,isup(dpdc),tmp)
                          addterm = .true.
                      end if
                  end if
                  if ( addterm ) then
                      call mult(tmp,hval(P(i)),hval(P(i)))
                      ind = P2(hlin(i),hcol(i),n)
                      call add(H(ind),hval(P(i)),H(ind))
                  end if
              end if
          end do

C     end if
      end do

      end

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

      subroutine evalpi(c,rho,lambda,equatn,p)

      implicit none

C     SCALAR ARGUMENTS
      logical equatn
      double precision lambda,p(2),rho,c(2)

C     LOCAL SCALARS
      double precision scltmp

C     LOCAL ARRAYS
      double precision tmp(2)

C     EXTERNAL FUNCTIONS
      double precision iinf,isup

      if ( equatn ) then
C         p = c * ( lambda + 0.5d0 * rho * c ) 
C           + 0.5d0 * lambda ** 2 / rho
          scltmp = 0.5d0 * rho
          call sclmlt(scltmp,c,p)
          call scladd(lambda,p,p)
          call mult(c,p,p)
          scltmp = 0.5d0 * lambda ** 2 / rho
          call scladd(scltmp,p,p)
      else
C         tmp = lambda + rho * c
          call sclmlt(rho,c,tmp)
          call scladd(lambda,tmp,tmp)

          if ( iinf(tmp) .ge. 0.0d0 ) then
C             p = c * ( lambda + 0.5d0 * rho * c ) 
C               + 0.5d0 * lambda ** 2 / rho
              scltmp = 0.5d0 * rho
              call sclmlt(scltmp,c,p)
              call scladd(lambda,p,p)
              call mult(c,p,p)
              scltmp = 0.5d0 * lambda ** 2 / rho
              call scladd(scltmp,p,p)

          else if (isup(tmp) .le. 0.0d0 ) then
C             p = 0.0d0
              call ivl1(0.0d0,p)

          else
C             It can be one or the other

C             p1 = c * ( lambda + 0.5d0 * rho * c )
C                + 0.5d0 * lambda ** 2 / rho
              scltmp = 0.5d0 * rho
              call sclmlt(scltmp,c,p)
              call scladd(lambda,p,p)
              call mult(c,p,p)
              scltmp = 0.5d0 * lambda ** 2 / rho
              call scladd(scltmp,p,p)

C             p2 = 0.0d0
              call ivl1(0.0d0,tmp)

C             p = ihull(p1,p2)
              call ihull(tmp,p,p)
          end if
      end if

      return

      end

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

      subroutine evalpi2(c,rho,lambda,equatn,p)

      implicit none

C     SCALAR ARGUMENTS
      logical equatn
      double precision lambda,p(2),rho,c(2)

C     LOCAL SCALARS
      double precision scltmp

C     LOCAL ARRAYS
      double precision tmp(2)

C     EXTERNAL FUNCTIONS
      double precision iinf,isup

      if ( equatn ) then
C         p  = c * ( lambda + 0.5d0 * rho * c )
          scltmp = 0.5d0 * rho
          call sclmlt(scltmp,c,p)
          call scladd(lambda,p,p)
          call mult(c,p,p)
      else
C         tmp = lambda + rho * c
          call sclmlt(rho,c,tmp)
          call scladd(lambda,tmp,tmp)

          if ( iinf(tmp) .ge. 0.0d0 ) then
C             p  = c * ( lambda + 0.5d0 * rho * c )
              scltmp = 0.5d0 * rho
              call sclmlt(scltmp,c,p)
              call scladd(lambda,p,p)
              call mult(c,p,p)

          else if (isup(tmp) .le. 0.0d0 ) then
C             p = - 0.5d0 * lambda ** 2 / rho
              scltmp = - 0.5d0 * lambda ** 2 / rho
              call ivl1(scltmp,p)

          else
C             It can be one or the other

C             p1  = c * ( lambda + 0.5d0 * rho * c )
              scltmp = 0.5d0 * rho
              call sclmlt(scltmp,c,p)
              call scladd(lambda,p,p)
              call mult(c,p,p)
c             write(*,*) 'parte 1: ',p(1),p(2)

C             p2 = - 0.5d0 * lambda ** 2 / rho
              scltmp = - 0.5d0 * lambda ** 2 / rho
              call ivl1(scltmp,tmp)
c             write(*,*) 'parte 2: ',tmp(1),tmp(2)

C             p = ihull(p1,p2)
              call ihull(tmp,p,p)
c             write(*,*) 'ihull: ',p(1),p(2)
          end if
      end if

      return

      end