!   E04STF Example Program Text
!   Mark 26 Release. NAG Copyright 2016.

!   NLP example : Linear objective + Linear constraint + Non-Linear constraint

    Module e04stfe_mod

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: confun, congrd, hess, mon, objfun,   &
                                          objgrd
    Contains
      Subroutine objfun(nvar,x,fx,inform,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
        Use nag_library, Only: f06eaf
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: cpuser
        Real (Kind=nag_wp), Intent (Out) :: fx
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: nvar
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Executable Statements ..
        fx = f06eaf(4,x,1,ruser,1)
        Return
      End Subroutine objfun
      Subroutine objgrd(nvar,x,nnzfd,fdx,inform,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: cpuser
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: nnzfd, nvar
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: fdx(nnzfd)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Executable Statements ..
        Continue
        fdx(1:nnzfd) = ruser(1:nnzfd)
        Return
      End Subroutine objgrd
      Subroutine confun(nvar,x,ncnln,gx,inform,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: cpuser
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: ncnln, nvar
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: gx(ncnln)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sqrt
!       .. Executable Statements ..
        Continue
        gx(1) = 12.0_nag_wp*x(1) + 11.9_nag_wp*x(2) + 41.8_nag_wp*x(3) +       &
          52.1_nag_wp*x(4) - 1.645_nag_wp*sqrt(.28_nag_wp*x(1)**2+.19_nag_wp*x &
          (2)**2+20.5_nag_wp*x(3)**2+.62_nag_wp*x(4)**2)
        Return
      End Subroutine confun
      Subroutine congrd(nvar,x,nnzgd,gdx,inform,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: cpuser
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: nnzgd, nvar
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: gdx(nnzgd)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: tmp
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sqrt
!       .. Executable Statements ..
        Continue
        tmp = sqrt(0.62_nag_wp*x(4)**2+20.5_nag_wp*x(3)**2+                    &
          0.19_nag_wp*x(2)**2+0.28_nag_wp*x(1)**2)
        gdx(1) = (12.0_nag_wp*tmp-0.4606_nag_wp*x(1))/tmp
        gdx(2) = (11.9_nag_wp*tmp-0.31255_nag_wp*x(2))/tmp
        gdx(3) = (41.8_nag_wp*tmp-33.7225_nag_wp*x(3))/tmp
        gdx(4) = (52.1_nag_wp*tmp-1.0199_nag_wp*x(4))/tmp
        Return
      End Subroutine congrd
      Subroutine hess(nvar,x,ncnln,idf,sigma,lambda,nnzh,hx,inform,iuser,      &
        ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: cpuser
        Real (Kind=nag_wp), Intent (In) :: sigma
        Integer, Intent (In)           :: idf, ncnln, nnzh, nvar
        Integer, Intent (Inout)        :: inform
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: hx(nnzh)
        Real (Kind=nag_wp), Intent (In) :: lambda(ncnln), x(nvar)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: tmp
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sqrt
!       .. Executable Statements ..
        inform = -1
        hx = 0.0_nag_wp
        Select Case (idf)
        Case (0)
          inform = 0
        Case (1,-1)
          tmp = sqrt(0.62_nag_wp*x(4)**2+20.5_nag_wp*x(3)**2+                  &
            0.19_nag_wp*x(2)**2+0.28_nag_wp*x(1)**2)
          tmp = tmp*(x(4)**2+33.064516129032258064_nag_wp*x(3)**2+             &
            0.30645161290322580645_nag_wp*x(2)**2+                             &
            0.45161290322580645161_nag_wp*x(1)**2)
!         1,1..4
          hx(1) = (-0.4606_nag_wp*x(4)**2-15.229516129032258064_nag_wp*x(3)**2 &
            -0.14115161290322580645_nag_wp*x(2)**2)/tmp
          hx(2) = (0.14115161290322580645_nag_wp*x(1)*x(2))/tmp
          hx(3) = (15.229516129032258064_nag_wp*x(1)*x(3))/tmp
          hx(4) = (0.4606_nag_wp*x(1)*x(4))/tmp
!         2,2..4
          hx(5) = (-0.31255_nag_wp*x(4)**2-10.334314516129032258_nag_wp*x(3)** &
            2-0.14115161290322580645_nag_wp*x(1)**2)/tmp
          hx(6) = (10.334314516129032258_nag_wp*x(2)*x(3))/tmp
          hx(7) = (0.31255_nag_wp*x(2)*x(4))/tmp
!         3,3..4
          hx(8) = (-33.7225_nag_wp*x(4)**2-10.334314516129032258_nag_wp*x(2)** &
            2-15.229516129032258065_nag_wp*x(1)**2)/tmp
          hx(9) = (33.7225_nag_wp*x(3)*x(4))/tmp
!         4,4
          hx(10) = (-33.7225_nag_wp*x(3)**2-0.31255_nag_wp*x(2)**2-            &
            0.4606_nag_wp*x(1)**2)/tmp
          If (idf==-1) Then
            hx = lambda(1)*hx
          End If
          inform = 0
        Case Default
        End Select
        Return
      End Subroutine hess
      Subroutine mon(nvar,x,nnzu,u,inform,rinfo,stats,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: cpuser
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: nnzu, nvar
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: rinfo(32), stats(32), u(nnzu),      &
                                          x(nvar)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Integer                        :: i, io_unit
!       .. Executable Statements ..
        Continue
        io_unit = iuser(2)
        Write (io_unit,99999)
        Write (io_unit,99998) x
        Write (io_unit,99997)
        Write (io_unit,99998) u
        Write (io_unit,99996)
        Write (io_unit,99995)(i,rinfo(i),i=1,32)
        Write (io_unit,99994)
        Write (io_unit,99995)(i,stats(i),i=1,32)
99999   Format ('Monitoring...',/,'   X(*)')
99998   Format (1P,E14.6)
99997   Format ('   U(*)')
99996   Format ('   RINFO(32)')
99995   Format (4(I2,1P,E14.6,1X))
99994   Format ('   STATS(32)')
        Return
      End Subroutine mon
    End Module e04stfe_mod

    Program e04stfe

!     .. Use Statements ..
      Use nag_library, Only: e04raf, e04ref, e04rgf, e04rhf, e04rjf, e04rkf,   &
                             e04rlf, e04ryf, e04rzf, e04stf, e04zmf, nag_wp,   &
                             x04acf
      Use e04stfe_mod, Only: confun, congrd, hess, mon, objfun, objgrd
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: solve_timeout = 5._nag_wp
      Integer, Parameter               :: nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: cpuser, handle
      Real (Kind=nag_wp)               :: bigbnd
      Integer                          :: i, idlc, idx, ifail, ilinear, j,     &
                                          nclin, ncnln, nnzu, nvar, ry_ifail
      Character (60)                   :: opt_s
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: b(8), bl(4), bu(4), linbl(2),        &
                                          linbu(2), nlnbl(1), nlnbu(1),        &
                                          rinfo(32), ruser(4), stats(32), x(4)
      Real (Kind=nag_wp), Allocatable  :: u(:)
      Integer                          :: icolb(8), icolgd(4), icolh(10),      &
                                          idxfd(4), irowb(8), irowgd(4),       &
                                          irowh(10), iuser(2)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: int
!     .. Executable Statements ..

      Write (nout,*) 'E04STF Example Program Results'
      Write (nout,*)
      Flush (nout)

      ifail = 0
      Call x04acf(66,'e04stf.out',1,ifail)
      Call x04acf(67,'e04stf.mon',1,ifail)
      Call x04acf(68,'e04stf.umon',1,ifail)

      Do ilinear = 0, 1
        handle = c_null_ptr
        bigbnd = 1.0E40_nag_wp
        nvar = 4
        nnzu = 0
        Call e04raf(handle,nvar,ifail)

        Write (opt_s,99992) bigbnd
        Call e04zmf(handle,opt_s,ifail)

!       Add simple bounds (x_i>=0).
        bl(1:4) = 0.0_nag_wp
        bu(1:4) = bigbnd
        nnzu = nnzu + 2*nvar
        Call e04rhf(handle,nvar,bl,bu,ifail)

        iuser(1) = ilinear

!       Add linear objective
        ruser(1:4) = (/24.55_nag_wp,26.75_nag_wp,39.00_nag_wp,40.50_nag_wp/)
        If (ilinear==1) Then
          Call e04ref(handle,4,ruser(1:4),ifail)
        Else
          idxfd(1:4) = (/1,2,3,4/)
          Call e04rgf(handle,4,idxfd,ifail)
        End If

!       Add two linear constraints
        nclin = 2
        nnzu = nnzu + 2*nclin
        linbl(1:2) = (/5.0_nag_wp,1.0_nag_wp/)
        linbu(1:2) = (/bigbnd,1.0_nag_wp/)
        irowb(1:8) = (/1,1,1,1,2,2,2,2/)
        icolb(1:8) = (/1,2,3,4,1,2,3,4/)
        b(1:8) = (/2.3_nag_wp,5.6_nag_wp,11.1_nag_wp,1.3_nag_wp,1.0_nag_wp,    &
          1.0_nag_wp,1.0_nag_wp,1.0_nag_wp/)
        idlc = 0
        Call e04rjf(handle,nclin,linbl,linbu,nclin*nvar,irowb,icolb,b,idlc,    &
          ifail)

!       Add one nonlinear constraint
        ncnln = 1
        nnzu = nnzu + 2*ncnln
        nlnbl(1:1) = (/21.0_nag_wp/)
        nlnbu(1:1) = (/bigbnd/)
        irowgd(1:4) = (/1,1,1,1/)
        icolgd(1:4) = (/1,2,3,4/)
        Call e04rkf(handle,ncnln,nlnbl,nlnbu,4,irowgd,icolgd,ifail)

!       Define dense structure of the Hessian
        idx = 1
        Do i = 1, nvar
          Do j = i, nvar
            icolh(idx) = j
            irowh(idx) = i
            idx = idx + 1
          End Do
        End Do

!       Hessian of nonlinear constraint
        Call e04rlf(handle,1,idx-1,irowh,icolh,ifail)
        If (ilinear/=1) Then
          Call e04rlf(handle,0,idx-1,irowh,icolh,ifail)
        End If

        Call e04ryf(handle,nout,'Overview',ifail)

!       call solver
        x = (/1.0_nag_wp,1.0_nag_wp,1.0_nag_wp,1.0_nag_wp/)
        Allocate (u(nnzu))

        iuser(2) = 68
        Call e04zmf(handle,'Monitoring File = 67',ifail)
        Call e04zmf(handle,'Monitoring Level = 5',ifail)
        Call e04zmf(handle,'Outer Iteration Limit = 26',ifail)
        Call e04zmf(handle,'Print File = 66',ifail)
        Call e04zmf(handle,'Print Level = 2',ifail)
        Call e04zmf(handle,'Stop Tolerance 1 = 2.5e-8',ifail)
        Call e04zmf(handle,'Time Limit = 60',ifail)

        ifail = -1
        Call e04stf(handle,objfun,objgrd,confun,congrd,hess,mon,nvar,x,nnzu,u, &
          rinfo,stats,iuser,ruser,cpuser,ifail)

        ry_ifail = 0
        Call e04ryf(handle,nout,'Options',ry_ifail)

        If (ifail==0) Then
          Write (nout,99999)
          Write (nout,99995)(i,x(i),i=1,nvar)
          Write (nout,99998)
          Write (nout,99993)(i,u(2*i-1),i,u(2*i),i=1,nvar)
          Write (nout,99997)
          Write (nout,99994)(i,u(2*i-1+2*nvar),i,u(2*i+2*nvar),i=1,nclin)
          Write (nout,99996)
          Write (nout,99994)(i,u(2*i-1+2*nvar+2*nclin),i,                      &
            u(2*i+2*nvar+2*nclin),i=1,ncnln)
          Write (nout,99991) rinfo(1)
          Write (nout,99990) rinfo(2)
          Write (nout,99989) rinfo(3)
          Write (nout,99988) rinfo(4)
          Write (nout,99987) rinfo(5)
          If (stats(8)<solve_timeout) Then
            Write (nout,99986)
          Else
            Write (nout,99985) stats(8)
          End If
          Write (nout,99984) int(stats(1))
          Write (nout,99983) int(stats(19))
          Write (nout,99982) int(stats(5))
          Write (nout,99981) int(stats(20))
          Write (nout,99980) int(stats(21))
          Write (nout,99979) int(stats(4))
        End If

        ifail = 0
        Call e04ryf(handle,nout,'Overview',ifail)

99999   Format (/,'Variables')
99998   Format ('Variable bound Lagrange multipliers')
99997   Format ('Linear constraints Lagrange multipliers')
99996   Format ('Nonlinear constraints Lagrange multipliers')
99995   Format (5X,'x(',I10,')',17X,'=',1P,E20.6)
99994   Format (5X,'l+(',I10,')',16X,'=',1P,E20.6,/,5X,'l-(',I10,')',16X,'=',  &
          1P,E20.6)
99993   Format (5X,'zL(',I10,')',16X,'=',1P,E20.6,/,5X,'zU(',I10,')',16X,'=',  &
          1P,E20.6)
99992   Format ('Infinite Bound Size               = ',1P,E14.6)
99991   Format ('At solution, Objective minimum     =',1P,E20.7)
99990   Format ('             Constraint violation  =',1P,E20.2)
99989   Format ('             Dual infeasibility    =',1P,E20.2)
99988   Format ('             Complementarity       =',1P,E20.2)
99987   Format ('             KKT error             =',1P,E20.2)
99986   Format ('Solved in allotted time limit')
99985   Format ('Solution took ',E10.3,' sec, which is longer than expected')
99984   Format ('    after iterations                        :',I11)
99983   Format ('    after objective evaluations             :',I11)
99982   Format ('    after objective gradient evaluations    :',I11)
99981   Format ('    after constraint evaluations            :',I11)
99980   Format ('    after constraint gradient evaluations   :',I11)
99979   Format ('    after hessian evaluations               :',I11)
        Deallocate (u)
!       release all memory held in the handle
        Call e04rzf(handle,ifail)
        Write (nout,*) '---------------------------------------------'
      End Do
    End Program e04stfe