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

      subroutine evalal(n,x,m,lambda,rho,equatn,linear,al,inform)

      implicit none

C     SCALAR ARGUMENTS
      double precision al
      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 Augmented Lagrangian function.

C     LOCAL SCALARS
      integer flag,i
      double precision c,p

      inform = 0

C     COMPUTE OBJECTIVE FUNTION

      call evalf(n,x,al,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))

      do i = 1,m

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

C         ADD P(c(i),rho(i),lambda(i))
          call evalp(c,rho(i),lambda(i),equatn(i),p)
          al = al + p

      end do

      end

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

      subroutine evalnal(n,x,m,lambda,rho,equatn,linear,nal,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer inform,m,n

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

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

C     PARAMETERS
      include "dim.f"

C     COMMON ARRAYS
      integer jcvar(jcnnzmax),jcsta(mmax),jclen(mmax)
      double precision dpdc(mmax),jcval(jcnnzmax)

C     LOCAL SCALARS
      integer flag,i,ind,j
      double precision c

C     COMMON BLOCKS
      common /jcdata/ dpdc,jcval,jcvar,jcsta,jclen

      inform = 0

C     COMPUTE THE GRADIENT OF THE OBJECTIVE FUNCTION

      call evalg(n,x,nal,flag)
      if ( flag .ne. 0 ) then
          inform = - 92
          return
      end if

C     COMPUTE \nabla L(x) = \nabla f(x) + \sum_{j=1}^m dPdcj * dcjdx

      ind = 0

      do j = 1,m

C         COMPUTE THE j-TH CONSTRAINT
          call evalc(n,x,j,c,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

C         COMPUTE dP/dc
          call evaldpdy(c,rho(j),lambda(j),equatn(j),dpdc(j))

          if ( equatn(j) .or. dpdc(j) .gt. 0.0d0 ) then

              jcsta(j) = ind + 1

C             COMPUTE THE GRADIENT OF THE j-TH CONSTRAINT
              call evaljac(n,x,j,jcvar(ind+1),jcval(ind+1),jclen(j),
     +        flag)
              if ( flag .ne. 0 ) then
                  inform = - 93
                  return
              end if

              ind = ind + jclen(j)

C             ADD dPdc * dcdx
              do i = jcsta(j),jcsta(j) + jclen(j) - 1
                  nal(jcvar(i)) = nal(jcvar(i)) + dpdc(j) * jcval(i)
              end do

          end if

      end do

      end

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

      subroutine evalhal(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     This subroutine cumputes the Hessian of the Augmented Lagrangian.

C     PARAMETERS
      include "dim.f"

C     COMMON ARRAYS
      integer jcvar(jcnnzmax),jcsta(mmax),jclen(mmax)
      double precision dpdc(mmax),jcval(jcnnzmax)

C     LOCAL SCALARS
      integer flag,hnnz,i,i1,i2,ind,ind1,ind2,j
      double precision c

C     LOCAL ARRAYS
      integer hcol(hnnzmax),hlin(hnnzmax)
      double precision hval(hnnzmax)

C     COMMON BLOCKS
      common /jcdata/ dpdc,jcval,jcvar,jcsta,jclen

      inform = 0

C     INITIALIZE DENSE AUGMENTED LAGRANGIAN HESSIAN

      do j = 1,n
          do i = j,n
              H(i,j) = 0.0d0
          end do
      end do

C     COMPUTE THE HESSIAN OF THE OBJECTIVE FUNCTION

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

      do i = 1,hnnz
          ind1 = hlin(i)
          ind2 = hcol(i)
          if ( ind1 .ge. ind2 ) then
              H(ind1,ind2) = H(ind1,ind2) + hval(i)
          end if
      end do

C     COMPUTE GRADIENTS AND HESSIANS OF THE CONSTRAINTS

      ind = 0

      do j = 1,m

          call evalc(n,x,j,c,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

C         COMPUTE dP/dc
          call evaldpdy(c,rho(j),lambda(j),equatn(j),dpdc(j))

          if ( equatn(j) .or. dpdc(j) .gt. 0.0d0 ) then

              jcsta(j) = ind + 1

C             COMPUTE THE GRADIENT OF THE j-TH CONSTRAINT
              call evaljac(n,x,j,jcvar(ind+1),jcval(ind+1),jclen(j),
     +        flag)
              if ( flag .ne. 0 ) then
                  inform = - 93
                  return
              end if

              ind = ind + jclen(j)

C             ADD FIRST-ORDER TERM
              do i1 = jcsta(j),jcsta(j) + jclen(j) - 1
                  ind1 = jcvar(i1)
                  do i2 = jcsta(j),jcsta(j) + jclen(j) - 1
                      ind2 = jcvar(i2)
                      if ( ind1 .ge. ind2 ) then
                         H(ind1,ind2) = 
     +                   H(ind1,ind2) + rho(j) * jcval(i1) * jcval(i2)
                      end if
                  end do
              end do

C             COMPUTE THE HESSIAN OF THE j-TH CONSTRAINT
              call evalhc(n,x,j,hlin,hcol,hval,hnnz,flag)
              if ( flag .ne. 0 ) then
                  inform = - 95
                  return
              end if

C             ADD SECOND-ORDER TERM
              do i = 1,hnnz
                  ind1 = hlin(i)
                  ind2 = hcol(i)
                  if ( ind1 .ge. ind2 ) then
                      H(ind1,ind2) = H(ind1,ind2) + dpdc(j) * hval(i)
                  end if
              end do

          end if

      end do

      end

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

      subroutine evalhal2(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     This subroutine cumputes the Hessian of the Augmented Lagrangian.

C     PARAMETERS
      include "dim.f"

C     COMMON ARRAYS
      integer jcvar(jcnnzmax),jcsta(mmax),jclen(mmax)
      double precision dpdc(mmax),jcval(jcnnzmax)

C     LOCAL SCALARS
      integer flag,hnnz,i,i1,i2,ind1,ind2,j

C     LOCAL ARRAYS
      integer hcol(hnnzmax),hlin(hnnzmax)
      double precision hval(hnnzmax)

C     COMMON BLOCKS
      common /jcdata/ dpdc,jcval,jcvar,jcsta,jclen

      inform = 0

C     INITIALIZE DENSE AUGMENTED LAGRANGIAN HESSIAN

      do j = 1,n
          do i = j,n
              H(i,j) = 0.0d0
          end do
      end do

C     COMPUTE THE HESSIAN OF THE OBJECTIVE FUNCTION

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

      do i = 1,hnnz
          ind1 = hlin(i)
          ind2 = hcol(i)
          if ( ind1 .ge. ind2 ) then
              H(ind1,ind2) = H(ind1,ind2) + hval(i)
          end if
      end do

C     COMPUTE GRADIENTS AND HESSIANS OF THE CONSTRAINTS

      do j = 1,m

          if ( equatn(j) .or. dpdc(j) .gt. 0.0d0 ) then

C             ADD FIRST-ORDER TERMS
              do i1 = jcsta(j),jcsta(j) + jclen(j) - 1
                  ind1 = jcvar(i1)
                  do i2 = jcsta(j),jcsta(j) + jclen(j) - 1
                      ind2 = jcvar(i2)
                      if ( ind1 .ge. ind2 ) then
                         H(ind1,ind2) = 
     +                   H(ind1,ind2) + rho(j) * jcval(i1) * jcval(i2)
                      end if
                  end do
              end do

C             COMPUTE THE HESSIAN OF THE j-TH CONSTRAINT
              call evalhc(n,x,j,hlin,hcol,hval,hnnz,flag)
              if ( flag .ne. 0 ) then
                  inform = - 95
                  return
              end if

C             ADD SECOND-ORDER TERMS
              do i = 1,hnnz
                  ind1 = hlin(i)
                  ind2 = hcol(i)
                  if ( ind1 .ge. ind2 ) then
                      H(ind1,ind2) = H(ind1,ind2) + dpdc(j) * hval(i)
                  end if
              end do

          end if

      end do

      end

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

      subroutine evalobjc(n,x,f,m,c,inform)

      implicit none

C     This subroutine computes the objective function and the 
C     constraints.

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

C     ARRAY ARGUMENTS
      double precision c(m),x(n)

C     Parameters of the subroutine:
C     =============================
C
C     On Entry:
C     =========
C
C     N integer: Number of variables.
C     ---------
C
C     X double precision x(n): Current point.
C     -----------------------
C
C     M integer: Number of constraints.
C     ---------
C
C     On Return:
C     ==========
C
C     F double precision: Objective function value at x.
C     ------------------
C
C     C double precision c(m): Constraints at x.
C     -----------------------
C
C     INFORM integer
C     --------------
C
C     0 means that the evaluation was successfuly done.
C
C     Any other value means that some error occured during the 
C     evaluation.

C     LOCAL SCALARS
      integer flag,i

      inform = 0

C     COMPUTE OBJECTIVE FUNTION

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

C     COMPUTE CONSTRAINTS

      do i = 1,m

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

      end do

      end

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

      subroutine evalp(y,rho,lambda,equatn,p)

      implicit none

C     SCALAR ARGUMENTS
      logical equatn
      double precision lambda,p,rho,y

      if ( equatn ) then
          p  = y * ( lambda + 0.5d0 * rho * y )
      else
          if ( lambda + rho * y .ge. 0.0d0 ) then
              p = y * ( lambda + 0.5d0 * rho * y )
          else
              p = - 0.5d0 * lambda ** 2 / rho
          end if
      end if

C     Added October 5th, 2006 to have that f(x) <= L(x).
      p = p + 0.5d0 * lambda ** 2 / rho

      end

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

      subroutine evaldpdy(y,rho,lambda,equatn,dpdy)

      implicit none

C     SCALAR ARGUMENTS
      logical equatn
      double precision y,rho,lambda,dpdy

      if ( equatn ) then
          dpdy = lambda + rho * y
      else
          dpdy = max( 0.0d0, lambda + rho * y )
      end if

      return

      end