NAG Library Manual, Mark 30.2
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
!   E04STF Example Program Text
!   Mark 30.2 Release. NAG Copyright 2024.

!   NLP example: linear objective + linear constraint + nonlinear constraint.
!   For illustrative purposes, the linear objective will be treated as a
!   nonlinear function to show the usage of objfun and objgrd user call-backs.

    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, objfun, objgrd
!     .. Derived Type Definitions ..
      Type, Public                     :: my_data
        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 ..
        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 ..
        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 ..
        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 (Inout) :: hx(nnzh), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: lambda(ncnln), x(nvar)
        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
        End Select
        Return
      End Subroutine hess
    End Module e04stfe_mod

    Program e04stfe
!     .. Use Statements ..
      Use e04stfe_mod, Only: confun, congrd, hess, 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, e04stu, e04zmf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: bigbnd = 1.0E20_nag_wp
      Integer, Parameter               :: 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 (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(100), ruser(1),      &
                                          stats(100), 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 ..

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

      ifail = 0
      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)

      Allocate (data)
      Allocate (data%coeffs(nvar),fdx(nvar),gdx(nnzgd))
      cpuser = c_loc(data)
      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/)
!     Add linear function as a nonlinear objective
      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/)

!     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 solver
      ifail = -1
      Call e04stf(handle,objfun,objgrd,confun,congrd,hess,e04stu,nvar,x,nnzu,  &
        u,rinfo,stats,iuser,ruser,cpuser,ifail)

      If (ifail==0) Then
        Write (nout,99999)
        Write (nout,99994)(i,x(i),i=1,nvar)
        Write (nout,99998)
        Write (nout,99991)(i,u(2*i-1),i,u(2*i),i=1,nvar)
        Write (nout,99996)
        Write (nout,99992)(i,u(2*i-1+2*nvar),i,u(2*i+2*nvar),i=1,nclin)
        Write (nout,99995)
        Write (nout,99992)(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,99997)
        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,99993) i, lmtest, rlmtest
        End Do

        Write (nout,99990) rinfo(1)
        Write (nout,99989) rinfo(2)
        Write (nout,99988) rinfo(3)
        Write (nout,99987) rinfo(4)
        Write (nout,99986) rinfo(5)
        Write (nout,99985) int(stats(1))
        Write (nout,99984) int(stats(19))
        Write (nout,99983) int(stats(5))
        Write (nout,99982) int(stats(20))
        Write (nout,99981) int(stats(21))
        Write (nout,99980) int(stats(4))
      End If

      Flush (nout)

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

99999 Format (/,'Variables')
99998 Format ('Variable bound Lagrange multipliers')
99997 Format ('Stationarity test for Lagrange multipliers')
99996 Format ('Linear constraints Lagrange multipliers')
99995 Format ('Nonlinear constraints Lagrange multipliers')
99994 Format (5X,'x(',I10,')',17X,'=',1P,E16.2)
99993 Format (4X,'lx(',I10,')',17X,'=',1P,E16.2,5X,A6)
99992 Format (4X,'zL(',I10,')',17X,'=',1P,E16.2,/,4X,'zU(',I10,')',17X,'=',1P, &
        E16.2)
99991 Format (4X,'zL(',I10,')',17X,'=',1P,E16.2,/,4X,'zU(',I10,')',17X,'=',1P, &
        E16.2)
99990 Format (/,'At solution, Objective minimum     =',1P,E16.4)
99989 Format ('             Constraint violation  =',1P,6X,E10.2)
99988 Format ('             Dual infeasibility    =',1P,6X,E10.2)
99987 Format ('             Complementarity       =',1P,6X,E10.2)
99986 Format ('             KKT error             =',1P,6X,E10.2)
99985 Format ('    after iterations                        :',I11)
99984 Format ('    after objective evaluations             :',I11)
99983 Format ('    after objective gradient evaluations    :',I11)
99982 Format ('    after constraint evaluations            :',I11)
99981 Format ('    after constraint gradient evaluations   :',I11)
99980 Format ('    after hessian evaluations               :',I11)
    End Program e04stfe