!   E04STF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2017.

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

    Module e04stfe_mod

!     Derived type my_data serves to demonstrate how to use the cpuser
!     communication argument in callbacks.

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: confun, congrd, hess, mon, objfun,   &
                                          objgrd
!     .. Derived Type Definitions ..
      Type, Public                     :: my_data
        Integer                        :: io_unit
        Real (Kind=nag_wp), Allocatable :: coeffs(:)
        Logical                        :: objfun_called = .False.
      End Type my_data
    Contains
      Subroutine objfun(nvar,x,fx,inform,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_f_pointer,    &
                                          c_ptr
        Use nag_library, Only: f06eaf
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: 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(*)
!       .. Local Scalars ..
        Type (my_data), Pointer        :: data
!       .. Executable Statements ..
        Call c_f_pointer(cptr=cpuser,fptr=data)
        fx = f06eaf(4,x,1,data%coeffs,1)
        data%objfun_called = .True.
        Return
      End Subroutine objfun
      Subroutine objgrd(nvar,x,nnzfd,fdx,inform,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_f_pointer,    &
                                          c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: 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(*)
!       .. Local Scalars ..
        Type (my_data), Pointer        :: data
!       .. Executable Statements ..
        Continue
        Call c_f_pointer(cptr=cpuser,fptr=data)
        fdx(1:nnzfd) = data%coeffs(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 (In)      :: 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 (In)      :: 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 (In)      :: 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_f_pointer,    &
                                          c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: 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 ..
        Type (my_data), Pointer        :: data
        Integer                        :: i
!       .. Executable Statements ..
        Continue
        Call c_f_pointer(cptr=cpuser,fptr=data)
        Write (data%io_unit,99999)
        Write (data%io_unit,99998) x
        Write (data%io_unit,99997)
        Write (data%io_unit,99998) u
        Write (data%io_unit,99996)
        Write (data%io_unit,99995)(i,rinfo(i),i=1,32)
        Write (data%io_unit,99994)
        Write (data%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 e04stfe_mod, Only: confun, congrd, hess, mon, my_data, objfun,       &
                             objgrd
      Use, Intrinsic                   :: iso_c_binding, Only: c_loc,          &
                                          c_null_ptr, c_ptr
      Use nag_library, Only: e04raf, e04ref, e04rgf, e04rhf, e04rjf, e04rkf,   &
                             e04rlf, e04ryf, e04rzf, e04stf, e04zmf, nag_wp,   &
                             x04acf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: solve_timeout = 10._nag_wp
      Integer, Parameter               :: nag_mon_file = 67, nout = 6,         &
                                          print_file = 66, user_mon_file = 68
!     .. Local Scalars ..
      Type (c_ptr)                     :: cpuser, handle
      Type (my_data), Pointer          :: data
      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(1), 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(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: int
!     .. Executable Statements ..

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

      ifail = 0
      Call x04acf(print_file,'e04stf.out',1,ifail)
      Call x04acf(nag_mon_file,'e04stf.mon',1,ifail)
      Call x04acf(user_mon_file,'e04stf.umon',1,ifail)

      Allocate (data)

      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)

!       Add linear objective
        Allocate (data%coeffs(nvar))
        data%coeffs(:) = (/24.55_nag_wp,26.75_nag_wp,39.00_nag_wp,             &
          40.50_nag_wp/)
        If (ilinear==1) Then
          Call e04ref(handle,4,data%coeffs,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
        Flush (nout)

        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))

        data%io_unit = user_mon_file
        Write (opt_s,99978) nag_mon_file
        Call e04zmf(handle,opt_s,ifail)
        Call e04zmf(handle,'Monitoring Level = 5',ifail)
        Call e04zmf(handle,'Outer Iteration Limit = 26',ifail)
        Write (opt_s,99977) print_file
        Call e04zmf(handle,opt_s,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)

        cpuser = c_loc(data)

        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
          If (.Not. data%objfun_called) Then
            Write (nout,*) 'Warning: expected a call to objfun'
          End If
          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
        Flush (nout)

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

        Deallocate (data%coeffs,u)
!       release all memory held in the handle
        Call e04rzf(handle,ifail)
        Write (nout,*) '---------------------------------------------'
      End Do

      Deallocate (data)

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)
99978 Format ('Monitoring File = ',I20)
99977 Format ('Print File = ',I20)
    End Program e04stfe