!   E04UDA Example Program Text
!   Mark 25 Release. NAG Copyright 2014.

    Module e04udae_mod

!     E04UDA Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                               :: confun, objfun
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Integer, Parameter, Public           :: inc1 = 1, iset = 1, lcwsav = 1,  &
                                              liwsav = 610, llwsav = 120,      &
                                              lrwsav = 475, nin = 5,           &
                                              ninopt = 7, nout = 6
    Contains
      Subroutine objfun(mode,n,x,objf,objgrd,nstate,iuser,ruser)
!       Routine to evaluate objective function and its 1st derivatives.

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: objf
        Integer, Intent (Inout)              :: mode
        Integer, Intent (In)                 :: n, nstate
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: objgrd(n), ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: x(n)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Executable Statements ..
        If (mode==0 .Or. mode==2) Then
          objf = x(1)*x(4)*(x(1)+x(2)+x(3)) + x(3)
        End If

        If (mode==1 .Or. mode==2) Then
          objgrd(1) = x(4)*(x(1)+x(1)+x(2)+x(3))
          objgrd(2) = x(1)*x(4)
          objgrd(3) = x(1)*x(4) + one
          objgrd(4) = x(1)*(x(1)+x(2)+x(3))
        End If

        Return

      End Subroutine objfun
      Subroutine confun(mode,ncnln,n,ldcj,needc,x,c,cjac,nstate,iuser,ruser)
!       Routine to evaluate the nonlinear constraints and their 1st
!       derivatives.

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: ldcj, n, ncnln, nstate
        Integer, Intent (Inout)              :: mode
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: c(ncnln)
        Real (Kind=nag_wp), Intent (Inout)   :: cjac(ldcj,n), ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: x(n)
        Integer, Intent (Inout)              :: iuser(*)
        Integer, Intent (In)                 :: needc(ncnln)
!       .. Executable Statements ..
        If (nstate==1) Then

!         First call to CONFUN.  Set all Jacobian elements to zero.
!         Note that this will only work when 'Derivative Level = 3'
!         (the default; see Section 11.2).

          cjac(1:ncnln,1:n) = zero
        End If

        If (needc(1)>0) Then

          If (mode==0 .Or. mode==2) Then
            c(1) = x(1)**2 + x(2)**2 + x(3)**2 + x(4)**2
          End If

          If (mode==1 .Or. mode==2) Then
            cjac(1,1) = x(1) + x(1)
            cjac(1,2) = x(2) + x(2)
            cjac(1,3) = x(3) + x(3)
            cjac(1,4) = x(4) + x(4)
          End If

        End If

        If (needc(2)>0) Then

          If (mode==0 .Or. mode==2) Then
            c(2) = x(1)*x(2)*x(3)*x(4)
          End If

          If (mode==1 .Or. mode==2) Then
            cjac(2,1) = x(2)*x(3)*x(4)
            cjac(2,2) = x(1)*x(3)*x(4)
            cjac(2,3) = x(1)*x(2)*x(4)
            cjac(2,4) = x(1)*x(2)*x(3)
          End If

        End If

        Return

      End Subroutine confun
    End Module e04udae_mod
    Program e04udae

!     E04UDA Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: dgemv, e04uca, e04uda, e04uea, e04wbf, nag_wp,    &
                             x04abf, x04acf, x04baf
      Use e04udae_mod, Only: confun, inc1, iset, lcwsav, liwsav, llwsav,       &
                             lrwsav, nin, ninopt, nout, objfun, one, zero
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Character (*), Parameter             :: fname = 'e04udae.opt'
!     .. Local Scalars ..
      Real (Kind=nag_wp)                   :: objf
      Integer                              :: i, ifail, inform, iter, j, lda,  &
                                              ldcj, ldr, liwork, lwork, mode,  &
                                              n, nclin, ncnln, outchn, sda,    &
                                              sdcjac
      Character (80)                       :: rec
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: a(:,:), bl(:), bu(:), c(:),      &
                                              cjac(:,:), clamda(:), objgrd(:), &
                                              r(:,:), work(:), x(:)
      Real (Kind=nag_wp)                   :: ruser(1), rwsav(lrwsav)
      Integer, Allocatable                 :: istate(:), iwork(:)
      Integer                              :: iuser(1), iwsav(liwsav)
      Logical                              :: lwsav(llwsav)
      Character (80)                       :: cwsav(lcwsav)
!     .. Intrinsic Procedures ..
      Intrinsic                            :: max
!     .. Executable Statements ..
      Write (rec,99991) 'E04UDA Example Program Results'
      Call x04baf(nout,rec)

!     Skip heading in data file
      Read (nin,*)

      Read (nin,*) n, nclin, ncnln
      liwork = 3*n + nclin + 2*ncnln
      lda = max(1,nclin)

      If (nclin>0) Then
        sda = n
      Else
        sda = 1
      End If

      ldcj = max(1,ncnln)

      If (ncnln>0) Then
        sdcjac = n
      Else
        sdcjac = 1
      End If

      ldr = n

      If (ncnln==0 .And. nclin>0) Then
        lwork = 2*n**2 + 20*n + 11*nclin
      Else If (ncnln>0 .And. nclin>=0) Then
        lwork = 2*n**2 + n*nclin + 2*n*ncnln + 20*n + 11*nclin + 21*ncnln
      Else
        lwork = 20*n
      End If

      Allocate (istate(n+nclin+ncnln),iwork(liwork),a(lda,sda), &
        bl(n+nclin+ncnln),bu(n+nclin+ncnln),c(max(1, &
        ncnln)),cjac(ldcj,sdcjac),clamda(n+nclin+ncnln),objgrd(n),r(ldr,n), &
        x(n),work(lwork))

      If (nclin>0) Then
        Read (nin,*)(a(i,1:sda),i=1,nclin)
      End If

      Read (nin,*) bl(1:(n+nclin+ncnln))
      Read (nin,*) bu(1:(n+nclin+ncnln))
      Read (nin,*) x(1:n)

!     Set the unit number for advisory messages to OUTCHN

      outchn = nout
      Call x04abf(iset,outchn)

!     Initialise E04UCA

      ifail = 0
      Call e04wbf('E04UCA',cwsav,lcwsav,lwsav,llwsav,iwsav,liwsav,rwsav, &
        lrwsav,ifail)

!     Set two options using E04UEA

      Call e04uea(' Infinite Bound Size = 1.0D+25 ',lwsav,iwsav,rwsav,inform)

      If (inform==0) Then

        Call e04uea(' Verify Level = -1 ',lwsav,iwsav,rwsav,inform)

      End If

      If (inform/=0) Then
        Write (rec,99999) 'E04UEA terminated with ' // 'INFORM = ', inform
        Call x04baf(nout,rec)
        Go To 100
      End If

!     Open the options file for reading

      mode = 0

      ifail = 0
      Call x04acf(ninopt,fname,mode,ifail)

!     Read the options file for the remaining options

      Call e04uda(ninopt,lwsav,iwsav,rwsav,inform)

      If (inform/=0) Then
        Write (rec,99999) 'E04UDA terminated with ' // 'INFORM = ', inform
        Call x04baf(nout,rec)
      End If

!     Solve the problem

      ifail = -1
      Call e04uca(n,nclin,ncnln,lda,ldcj,ldr,a,bl,bu,confun,objfun,iter, &
        istate,c,cjac,clamda,objf,objgrd,r,x,iwork,liwork,work,lwork,iuser, &
        ruser,lwsav,iwsav,rwsav,ifail)

      Select Case (ifail)
      Case (0:6,8)
        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,99998)
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)

        Do i = 1, n
          Write (rec,99997) i, istate(i), x(i), clamda(i)
          Call x04baf(nout,rec)
        End Do

        If (nclin>0) Then

!         A*x --> work.
!         The NAG name equivalent of dgemv is f06paf
          Call dgemv('N',nclin,n,one,a,lda,x,inc1,zero,work,inc1)

          Write (rec,'()')
          Call x04baf(nout,rec)
          Write (rec,'()')
          Call x04baf(nout,rec)
          Write (rec,99996)
          Call x04baf(nout,rec)
          Write (rec,'()')
          Call x04baf(nout,rec)

          Do i = n + 1, n + nclin
            j = i - n
            Write (rec,99995) j, istate(i), work(j), clamda(i)
            Call x04baf(nout,rec)
          End Do

        End If

        If (ncnln>0) Then
          Write (rec,'()')
          Call x04baf(nout,rec)
          Write (rec,'()')
          Call x04baf(nout,rec)
          Write (rec,99994)
          Call x04baf(nout,rec)
          Write (rec,'()')
          Call x04baf(nout,rec)

          Do i = n + nclin + 1, n + nclin + ncnln
            j = i - n - nclin
            Write (rec,99993) j, istate(i), c(j), clamda(i)
            Call x04baf(nout,rec)
          End Do

        End If

        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,'()')
        Call x04baf(nout,rec)
        Write (rec,99992) objf
        Call x04baf(nout,rec)
      End Select

100   Continue

99999 Format (1X,A,I5)
99998 Format (1X,'Varbl',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99997 Format (1X,'V',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99996 Format (1X,'L Con',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99995 Format (1X,'L',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99994 Format (1X,'N Con',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99993 Format (1X,'N',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99992 Format (1X,'Final objective value = ',1P,G15.7)
99991 Format (1X,A)
    End Program e04udae