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

NAG FL Interface Introduction
Example description
!   E04SRF Example Program Text
!   Mark 30.0 Release. NAG Copyright 2024.

!   NLP example: Nonlinear objective, linear constraint and two nonlinear
!   constraints.

    Module e04srfe_mod
!     .. Use Statements ..
      Use nag_precisions, Only: wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: confun, congrd, objfun, objgrd
    Contains
      Subroutine objfun(nvar,x,fx,inform,iuser,ruser,cpuser)
!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: cpuser
        Real (Kind=wp), Intent (Out)   :: fx
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: nvar
!       .. Array Arguments ..
        Real (Kind=wp), Intent (Inout) :: ruser(*)
        Real (Kind=wp), Intent (In)    :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: cos
!       .. Executable Statements ..

        Continue

        fx = (x(1)+x(2)+x(3))**2 + 3.0_wp*x(3) + 5.0_wp*x(4) +                 &
          cos(0.01_wp*x(1)) - 1.0_wp

        inform = 0

      End Subroutine objfun

      Subroutine objgrd(nvar,x,nnzfd,fdx,inform,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: cpuser
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: nnzfd, nvar
!       .. Array Arguments ..
        Real (Kind=wp), Intent (Inout) :: fdx(nnzfd), ruser(*)
        Real (Kind=wp), Intent (In)    :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=wp)                 :: s
!       .. Intrinsic Procedures ..
        Intrinsic                      :: sin
!       .. Executable Statements ..

        Continue

        s = 2.0_wp*(x(1)+x(2)+x(3))

        fdx(1) = s - 0.01*sin(x(1))
        fdx(2) = s
        fdx(3) = s + 3.0_wp
        fdx(4) = 5.0_wp

        inform = 0

      End Subroutine objgrd

      Subroutine confun(nvar,x,ng,gx,inform,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: cpuser
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: ng, nvar
!       .. Array Arguments ..
        Real (Kind=wp), Intent (Out)   :: gx(ng)
        Real (Kind=wp), Intent (Inout) :: ruser(*)
        Real (Kind=wp), Intent (In)    :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Executable Statements ..

        Continue

        gx(1) = x(1)**2 + x(2)**2 + x(3)
        gx(2) = x(2)**4 + x(4)

        inform = 0

      End Subroutine confun

      Subroutine congrd(nvar,x,nnzfd,gdx,inform,iuser,ruser,cpuser)

!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Implicit None Statement ..
        Implicit None
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (In)      :: cpuser
        Integer, Intent (Inout)        :: inform
        Integer, Intent (In)           :: nnzfd, nvar
!       .. Array Arguments ..
        Real (Kind=wp), Intent (Inout) :: gdx(nnzfd), ruser(*)
        Real (Kind=wp), Intent (In)    :: x(nvar)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Executable Statements ..

        Continue

        gdx(1) = 2.0_wp*x(1)
        gdx(2) = 2.0_wp*x(2)
        gdx(3) = 1.0_wp
        gdx(4) = 4.0_wp*x(2)**3
        gdx(5) = 1.0_wp

        inform = 0

      End Subroutine congrd

    End Module e04srfe_mod

    Program e04srfe

!     .. Use Statements ..
      Use e04srfe_mod, Only: confun, congrd, objfun, objgrd
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nag_library, Only: e04raf, e04rcf, e04rgf, e04rhf, e04rjf, e04rkf,   &
                             e04rzf, e04srf, e04sru, e04srz, e04zmf
      Use nag_precisions, Only: wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=wp), Parameter        :: bigbnd = 1.0E20_wp
      Integer, Parameter               :: nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: cpuser, handle
      Integer                          :: i, idlc, ifail, inform, nclin,       &
                                          ncnln, nnzb, nnzfd, nnzgd, nnzu,     &
                                          nvar
      Character (6)                    :: rlmtest
!     .. Local Arrays ..
      Real (Kind=wp), Allocatable      :: b(:), blx(:), bux(:), fdx(:), fl(:), &
                                          fu(:), gdx(:), lambda(:), lcbl(:),   &
                                          lcbu(:), lmtest(:), u(:), x(:)
      Real (Kind=wp)                   :: rinfo(100), ruser(0), stats(100)
      Integer, Allocatable             :: icolb(:), icolgd(:), iidx(:),        &
                                          irowb(:), irowgd(:)
      Integer                          :: iuser(0)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, int, merge
!     .. Executable Statements ..

      ifail = 0
      cpuser = c_null_ptr

      Write (nout,Fmt=99999) 'E04SRF Example Program Results'

!     Problem size
      nvar = 4
!     Counter for Lagrange multipliers
      nnzu = 0
      Allocate (x(nvar),blx(nvar),bux(nvar))

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

!     Add simple box bounds on x
      blx(1:nvar) = (/-bigbnd,-bigbnd,0.0_wp,0.0_wp/)
      bux(1:nvar) = bigbnd
      ifail = 0
      Call e04rhf(handle,nvar,blx,bux,ifail)
!     Specify the amount of Lagrange mult. required
      nnzu = 2*nvar

!     Add nonlinear objective gradient information
      nnzfd = nvar
      Allocate (iidx(nnzfd),fdx(nnzfd))
      iidx(:) = (/(i,i=1,nvar)/)
      ifail = 0
      Call e04rgf(handle,nnzfd,iidx,ifail)

!     Add two nonlinear constraints
      ncnln = 2
      Allocate (fl(ncnln),fu(ncnln))
      fl(:) = (/2.0_wp,4.0_wp/)
      fu(:) = (/2.0_wp,4.0_wp/)
!     Number of nonzero elements in the constraint Jacobian
      nnzgd = 5
      Allocate (irowgd(nnzgd),icolgd(nnzgd),gdx(nnzgd))
      irowgd(:) = (/1,1,1,2,2/)
      icolgd(:) = (/1,2,3,2,4/)
!     Add nonlinear constraint information
      ifail = 0
      Call e04rkf(handle,ncnln,fl,fu,nnzgd,irowgd,icolgd,ifail)
!     Update the Lagrange mult. count
      nnzu = nnzu + 2*ncnln

!     Add one linear constraint
      nclin = 1
!     Number of nonzero elements in the constraint
      nnzb = 2
      Allocate (lcbl(nclin),lcbu(nclin),irowb(nnzb),icolb(nnzb),b(nnzb))
      lcbl(1) = 0.0_wp
      lcbu(1) = bigbnd
      irowb(:) = (/1,1/)
      icolb(:) = (/1,2/)
      b(:) = (/2.0_wp,4.0_wp/)
      idlc = 0
      ifail = 0
      Call e04rjf(handle,nclin,lcbl,lcbu,nnzb,irowb,icolb,b,idlc,ifail)
!     Update the Lagrange mult. count
      nnzu = nnzu + 2*nclin

      Allocate (u(nnzu))

!     Optionally, define variable x(4) to be linear
!     Hinting which variables are linear in problems with many
!     variables can speed up performance
      Call e04rcf(handle,'linear',1,(/4/),ifail)

!     Define initial guess point
      x(1:nvar) = (/1.0_wp,2.0_wp,3.0_wp,4.0_wp/)

!     Add options
      ifail = 0
!     Disable printing
      Call e04zmf(handle,'Print Level = 0',ifail)
!     Do not print options
      Call e04zmf(handle,'Print Options = No',ifail)
!     It is recommended on new problems to always verify the derivatives
      Call e04zmf(handle,'Verify Derivatives = Yes',ifail)
!     Do not print solution, x and f(x) will be printed afterwards
      Call e04zmf(handle,'Print Solution = No',ifail)

!     Solve the problem
      ifail = -1
      Call e04srf(handle,objfun,objgrd,confun,congrd,e04srz,e04sru,nvar,x,     &
        nnzu,u,rinfo,stats,iuser,ruser,cpuser,ifail)

      Allocate (lambda(nnzu/2),lmtest(nnzfd))
      If (ifail==0) Then
!       Print solution point and objective function value
        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)

!       Gradient of the objective
        inform = 0
        Call objgrd(nvar,x,nnzfd,fdx,inform,iuser,ruser,cpuser)

!       Gradient of the linear constraints is already available in
!       irowb(1:2), icolb(1:2) and b(1:2)

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

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

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

      Flush (nout)

!     Clean up
      ifail = 0
      Call e04rzf(handle,ifail)

      Deallocate (b,blx,bux,fl,fu,lcbl,lcbu,u,x,lambda,fdx,gdx,lmtest,icolb,   &
        icolgd,iidx,irowb,irowgd)

99999 Format (A)
99998 Format (/,'Variables')
99997 Format ('Variable bound Lagrange multipliers')
99996 Format ('Stationarity test for Lagrange multipliers')
99995 Format ('Linear constraint 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 (/,'Solution found. Objective minimum  =',1P,E16.4)
99988 Format ('             Constraint violation  =',1P,6X,E10.2)
99987 Format ('             Dual infeasibility    =',1P,6X,E10.2)
99986 Format ('   Iterations                      =',5X,I11)
99985 Format ('   Objective evaluations           =',5X,I11)
99984 Format ('   Objective gradient evaluations  =',5X,I11)
99983 Format ('   Constraint evaluations          =',5X,I11)
99982 Format ('   Constraint gradient evaluations =',5X,I11)
99981 Format ('   Hessian evaluations             =',5X,I11)

    End Program e04srfe