Example description
!   E04STF Example Program Text
!   Mark 27.1 Release. NAG Copyright 2020.

!   NLP example: Linear objective + Linear constraint + Nonlinear constraint
!   For illustrative purposes, the linear objective will the trated as a
!   nonlinear function to show the usage of objfun and objgrd user
!   callbacks.

    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, monit, objfun, &
                                          objgrd
!     .. Derived Type Definitions ..
      Type, Public                     :: my_data
        Integer                        :: mon
        Real (Kind=nag_wp), Allocatable :: coeffs(:)
      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)
        inform = 0
        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 (Inout) :: fdx(nnzfd), 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)
        inform = 0
        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
        inform = 0
        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 (Inout) :: gdx(nnzgd), 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
        inform = 0
        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 monit(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
!       .. Intrinsic Procedures ..
        Intrinsic                      :: int
!       .. Executable Statements ..
        Continue
        Call c_f_pointer(cptr=cpuser,fptr=data)
        Write (data%mon,99999)
        Write (data%mon,99998)
        Write (data%mon,99997)(i,x(i),i=1,4)
        Write (data%mon,99995)
        Write (data%mon,99996)(i,u(2*i-1),u(2*i),i=1,7)
        Write (data%mon,99994)
        Write (data%mon,99993) 'Objective function value f(x)              ',  &
          rinfo(1)
        Write (data%mon,99993) 'Constraint violation (primal infeasibility)',  &
          rinfo(2)
        Write (data%mon,99993) 'Dual infeasibility                         ',  &
          rinfo(3)
        Write (data%mon,99993) 'Complementarity                            ',  &
          rinfo(4)
        Write (data%mon,99993) 'Karush-Kuhn-Tucker error                   ',  &
          rinfo(5)
        Write (data%mon,99991)
        Write (data%mon,99992) 'Number of the iteration                    ',  &
          int(stats(1))
        Write (data%mon,99992) 'Number of backtracking trial steps         ',  &
          int(stats(3))
        Write (data%mon,99990)

99999   Format (/,4('-'),2X,'Begin monitoring info',2X,32('-'))
99998   Format (/,'   Current point  X(:)',/,'  idx            value')
99997   Format (3X,I2,3X,E14.6)
99996   Format (3X,I2,3X,E14.6,3X,E14.6)
99995   Format (/,'   Lagrange mult. U(:)',/,'const.         L value',10X,     &
          'U value')
99994   Format (/,'   Rinfo:')
99993   Format (3X,A40,3X,'=',E14.6)
99992   Format (3X,A40,3X,'=',I14)
99991   Format (/,'   Stats:')
99990   Format (/,'----  End monitoring info  ',34('-'))
        Return
      End Subroutine monit
    End Module e04stfe_mod

    Program e04stfe
!     .. Use Statements ..
      Use e04stfe_mod, Only: confun, congrd, hess, monit, my_data, objfun,     &
                             objgrd
      Use, Intrinsic                   :: iso_c_binding, Only: c_loc,          &
                                          c_null_ptr, c_ptr
      Use nag_library, Only: e04raf, e04rgf, e04rhf, e04rjf, e04rkf, e04rlf,   &
                             e04rzf, e04stf, e04zmf, nag_wp, x04acf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: bigbnd = 1.0E20_nag_wp
      Integer, Parameter               :: log_file = 67, mon_file = 68,        &
                                          nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: cpuser, handle
      Type (my_data), Pointer          :: data
      Real (Kind=nag_wp)               :: lmtest
      Integer                          :: i, idlc, idx, ifail, j, nclin,       &
                                          ncnln, nnzfd, nnzgd, nnzu, nvar
      Character (60)                   :: opt_s
      Character (6)                    :: rlmtest
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: b(8), bl(4), bu(4), lambda(7),       &
                                          linbl(2), linbu(2), nlnbl(1),        &
                                          nlnbu(1), rinfo(32), ruser(1),       &
                                          stats(32), x(4)
      Real (Kind=nag_wp), Allocatable  :: fdx(:), gdx(:), u(:)
      Integer                          :: icolb(8), icolgd(4), icolh(10),      &
                                          idxfd(4), irowb(8), irowgd(4),       &
                                          irowh(10), iuser(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, int, merge
!     .. Executable Statements ..
      Continue

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

      ifail = 0
!     Specify the file to save user call-back monitoring information
      Call x04acf(mon_file,'e04stf.mon',1,ifail)
!     Specify the file to save verbose iteration log
      Call x04acf(log_file,'e04stf.log',1,ifail)

      Allocate (data)

      handle = c_null_ptr
!     Problem size
      nvar = 4
!     Counter for Lagrange multipliers
      nnzu = 0
!     Objective gradient nonzero elements quantity
      nnzfd = 4
!     Constraint jacobian nonzero elements quantity
      nnzgd = 4

!     Initialize handle
      Call e04raf(handle,nvar,ifail)

!     Add linear function as a nonlinear objective
      Allocate (data%coeffs(nvar),fdx(nvar),gdx(nnzgd))
      cpuser = c_loc(data)
      data%mon = mon_file
      data%coeffs(:) = (/24.55_nag_wp,26.75_nag_wp,39.00_nag_wp,40.50_nag_wp/)
      idxfd(1:nnzfd) = (/1,2,3,4/)
      Call e04rgf(handle,nnzfd,idxfd,ifail)

!     Add simple bounds (x>=0).
      bl(1:4) = 0.0_nag_wp
      bu(1:4) = bigbnd
      Call e04rhf(handle,nvar,bl,bu,ifail)
!     Specify the amount of Lagrange mult. required
      nnzu = nnzu + 2*nvar

!     Add two linear constraints
      nclin = 2
      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)
!     Update the Lagrange mult. count
      nnzu = nnzu + 2*nclin

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

!     Update the Lagrange mult. count and allocate space for storage
      nnzu = nnzu + 2*ncnln
      Allocate (u(nnzu))

!     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)
      Call e04rlf(handle,0,idx-1,irowh,icolh,ifail)

!     Specify initial starting point
      x(1:nvar) = (/1.0_nag_wp,1.0_nag_wp,1.0_nag_wp,1.0_nag_wp/)

!     Specify to solver where to store the monitoring iteration log
      Write (opt_s,99978) log_file
      Call e04zmf(handle,opt_s,ifail)
!     Set  print level for monitoring iteration output
      Call e04zmf(handle,'Monitoring Level = 2',ifail)
!     Set  print level for iteration output
      Call e04zmf(handle,'Print Level = 0',ifail)
      Call e04zmf(handle,'Outer Iteration Limit = 26',ifail)
      Call e04zmf(handle,'Stop Tolerance 1 = 2.5e-8',ifail)
      Call e04zmf(handle,'Time Limit = 10',ifail)

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

      If (ifail==0) Then
        Write (nout,99998)
        Write (nout,99993)(i,x(i),i=1,nvar)
        Write (nout,99997)
        Write (nout,99990)(i,u(2*i-1),i,u(2*i),i=1,nvar)
        Write (nout,99995)
        Write (nout,99991)(i,u(2*i-1+2*nvar),i,u(2*i+2*nvar),i=1,nclin)
        Write (nout,99994)
        Write (nout,99991)(i,u(2*i-1+2*nvar+2*nclin),i,u(2*i+2*nvar+2*nclin),  &
          i=1,ncnln)

!       Check of Lagrange multipliers (complementariety)
!       Gradient of the objective (also available in data%coeffs(1:nnzfd))
        Call objgrd(nvar,x,nnzfd,fdx,ifail,iuser,ruser,cpuser)
!       Gradient of the linear constraints is already available in
!       irowb(1:8), icolb(1:8) and b(1:8)

!       Gradient of the nonlinear constraint
        Call congrd(nvar,x,nnzgd,gdx,ifail,iuser,ruser,cpuser)

!       Obtain the multipliers with correct sign
!       4 bound constraints, 2 linear constraints, and 1 nonlinear constraint
        Do i = 1, 7
          lambda(i) = u(2*i) - u(2*i-1)
        End Do
!       lambda (1..4) mult. related to bounds
!       lambda (5..6) mult. related to linear constraints
!       lambda (7) mult. related to the nonlinear constraint
        Write (nout,99996)
        Do i = 1, 4
          lmtest = fdx(i) + lambda(i) + lambda(5)*b(i) + lambda(6)*b(4+i) +    &
            lambda(7)*gdx(i)
          rlmtest = merge('Ok    ','Failed',abs(lmtest)<2.5E-8_nag_wp)
          Write (nout,99992) i, lmtest, rlmtest
        End Do

        Write (nout,99989) rinfo(1)
        Write (nout,99988) rinfo(2)
        Write (nout,99987) rinfo(3)
        Write (nout,99986) rinfo(4)
        Write (nout,99985) rinfo(5)
        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))

        Write (nout,99999)                                                     &
          'See e04stf.mon file for monitoring information, and',               &
          '    e04stf.log file for a detailed iteration log.'
      End If

      Flush (nout)

!     Release all memory
      Call e04rzf(handle,ifail)
      Deallocate (data%coeffs,u,fdx,gdx)
      Deallocate (data)

99999 Format (/,2(A,/))
99998 Format (/,'Variables')
99997 Format ('Variable bound Lagrange multipliers')
99996 Format ('Stationarity test for Lagrange multipliers')
99995 Format ('Linear constraints Lagrange multipliers')
99994 Format ('Nonlinear constraints Lagrange multipliers')
99993 Format (5X,'x(',I10,')',17X,'=',1P,E16.2)
99992 Format (4X,'lx(',I10,')',17X,'=',1P,E16.2,5X,A6)
99991 Format (4X,'zL(',I10,')',17X,'=',1P,E16.2,/,4X,'zU(',I10,')',17X,'=',1P, &
        E16.2)
99990 Format (4X,'zL(',I10,')',17X,'=',1P,E16.2,/,4X,'zU(',I10,')',17X,'=',1P, &
        E16.2)
99989 Format (/,'At solution, Objective minimum     =',1P,E16.4)
99988 Format ('             Constraint violation  =',1P,6X,E10.2)
99987 Format ('             Dual infeasibility    =',1P,6X,E10.2)
99986 Format ('             Complementarity       =',1P,6X,E10.2)
99985 Format ('             KKT error             =',1P,6X,E10.2)
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)
    End Program e04stfe