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

NAG FL Interface Introduction
Example description
!   H02DD Example Program Text

!   Mark 30.3 Release. nAG Copyright 2024.

    Module h02ddfe_mod

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: confun, congrd, objfun, objgrd

    Contains

      Subroutine objfun(nvar,x,fx,inform,iuser,ruser,cpuser)

!       .. 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(*)
!       .. Executable Statements ..

        fx = x(1)*(4.0_nag_wp*x(1)+3.0_nag_wp*x(2)-x(3)) +                     &
          x(2)*(3.0_nag_wp*x(1)+6.0_nag_wp*x(2)+x(3)) +                        &
          x(3)*(x(2)-x(1)+10.0_nag_wp*x(3))

      End Subroutine objfun

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

!       .. 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(*)
!       .. Executable Statements ..

!           Objective gradients for continuous variables
        fdx(1) = 8.0_nag_wp*x(1) + 6.0_nag_wp*x(2) - 2.0_nag_wp*x(3)
        fdx(2) = 6.0_nag_wp*x(1) + 12.0_nag_wp*x(2) + 2.0_nag_wp*x(3)
        fdx(3) = 2.0_nag_wp*(x(2)-x(1)) + 20.0_nag_wp*x(3)
        fdx(4) = 0.0_nag_wp

      End Subroutine objgrd

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

!       .. 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(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: rho
        Integer                        :: p
!       .. Intrinsic Procedures ..
        Intrinsic                      :: real
!       .. Executable Statements ..

!       Constraints
        p = iuser(1)
        rho = ruser(1)

!       Mean return rho:
        gx(1) = 8.0_nag_wp*x(1) + 9.0_nag_wp*x(2) + 12.0_nag_wp*x(3) +         &
          7.0_nag_wp*x(4) - rho
!       Maximum of p assets in portfolio:
        gx(2) = real(p,kind=nag_wp) - x(5) - x(6) - x(7) - x(8)

      End Subroutine confun

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

!       .. 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(*)
!       .. Executable Statements ..

!       Jacobian
        gdx(1:4) = (/8.0_nag_wp,9.0_nag_wp,12.0_nag_wp,7.0_nag_wp/)

!       (Constraint 2 does not include continuous variables, hence their
!       derivatives are not required)

      End Subroutine congrd

    End Module h02ddfe_mod

    Program h02ddfe

!     .. Use Statements ..
      Use h02ddfe_mod, Only: confun, congrd, objfun, objgrd
      Use iso_c_binding, Only: c_null_ptr, c_ptr
      Use nag_library, Only: e04raf, e04rcf, e04rgf, e04rhf, e04rjf, e04rkf,   &
                             e04rzf, e04zmf, h02ddf, h02ddu, nag_wp, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: bigbnd = 1.0E6_nag_wp
      Real (Kind=nag_wp), Parameter    :: infbnd = 1.0E20_nag_wp
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter               :: nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: cpuser, handle
      Integer                          :: i, idlc, ifail, lgroup, nclin,       &
                                          ncnln, nnzb, nnzfd, nnzgd, nvar
      Character (32)                   :: ptype
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: b(:), linbl(:), linbu(:), nlnbl(:),  &
                                          nlnbu(:), ruser(:), varbl(:),        &
                                          varbu(:), x(:)
      Real (Kind=nag_wp)               :: rinfo(100), stats(100)
      Integer, Allocatable             :: group(:), icolb(:), icolgd(:),       &
                                          idxfd(:), irowb(:), irowgd(:),       &
                                          iuser(:)
!     .. Executable Statements ..

      Write (nout,*) 'H02DDF Example Program Results'
      Write (nout,*)

!     Initialise handle
      nvar = 8
      ifail = 0
      Call e04raf(handle,nvar,ifail)

!     Set variable types
!     1:4 continuous, 5:8 binary
      ptype = 'BINARY'
      lgroup = 4
      Allocate (group(lgroup))
      group(1:lgroup) = (/(i,i=5,8)/)
      ifail = 0
      Call e04rcf(handle,ptype,lgroup,group,ifail)

!     Set non-linear objective
      nnzfd = 4
      Allocate (idxfd(nnzfd))
      idxfd(1:nnzfd) = (/(i,i=1,4)/)
      ifail = 0
      Call e04rgf(handle,nnzfd,idxfd,ifail)

!     Set variable bounds
      Allocate (varbl(nvar),varbu(nvar))
!     xi >= 0
      varbl(1:4) = zero
      varbu(1:4) = bigbnd
!     yi already defined as binary
      varbl(5:8) = zero
      varbu(5:8) = one
      ifail = 0
      Call e04rhf(handle,nvar,varbl,varbu,ifail)

!     Construct linear constraints
      nclin = 5
      nnzb = 12
      Allocate (b(nnzb),irowb(nnzb),icolb(nnzb),linbl(nclin),linbu(nclin))
!     sum(x1:x4) = 1; i.e. 1 <= x1+x2+x3+x4 <= 1
      linbl(1) = one
      linbu(1) = one
      irowb(1:4) = (/1,1,1,1/)
      icolb(1:4) = (/1,2,3,4/)
      b(1:4) = (/one,one,one,one/)
!     x1 <= x5; i.e. 0 <= x5-x1 <= inf
!     x2 <= x6; i.e. 0 <= x6-x2 <= inf
!     x3 <= x7; i.e. 0 <= x7-x3 <= inf
!     x4 <= x8; i.e. 0 <= x8-x4 <= inf
      linbl(2:5) = zero
      linbu(2:5) = infbnd
      irowb(5:12) = (/2,2,3,3,4,4,5,5/)
      icolb(5:12) = (/1,5,2,6,3,7,4,8/)
      b(5:12) = (/-one,one,-one,one,-one,one,-one,one/)

!     Add linear constraints
      idlc = 0
      ifail = 0
      Call e04rjf(handle,nclin,linbl,linbu,nnzb,irowb,icolb,b,idlc,ifail)

!     Construct non-linear constraints
      ncnln = 2
      Allocate (nlnbl(ncnln),nlnbu(ncnln))
!     r*x = rho:  g1(x) = r*x - rho;  0 <= g1(x) <= 0
      nlnbl(1) = zero
      nlnbu(1) = zero
!     sum(x5:x8) <= p:  g2(x) = p - sum(x5:x8);  0 <= g2(x) <= inf
      nlnbl(2) = zero
      nlnbu(2) = infbnd

!     Specify sparsity pattern of Jacobian from congrd
!     - congrd returns non-zero values for g1, x1->x4 only
!     - x5:x8 are binary and g2 only uses binary variables
      nnzgd = 4
      Allocate (icolgd(nnzgd),irowgd(nnzgd))
      irowgd(1:4) = 1
      icolgd(1:4) = (/(i,i=1,4)/)

!     Add non-linear constraints
      ifail = 0
      Call e04rkf(handle,ncnln,nlnbl,nlnbu,nnzgd,irowgd,icolgd,ifail)

!     Initial estimate (binary variables need not be given)
      Allocate (x(nvar))
      x(1:4) = one
      x(5:8) = zero

!     Set up user workspace
      Allocate (iuser(1),ruser(1))
      cpuser = c_null_ptr

!     Add portfolio parameters p=3 and rho=10 for confun to access
      iuser(1) = 3
      ruser(1) = 10.0_nag_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
      ifail = -1
      Call h02ddf(handle,objfun,objgrd,confun,congrd,h02ddu,nvar,x,rinfo,      &
        stats,iuser,ruser,cpuser,ifail)

      If (ifail==0) Then
!       Print final estimate
        Call x04caf('G','N',nvar,1,x,nvar,'Final estimate:',ifail)
        Write (nout,'(/1x,a,1x,g12.4)') 'Optimised value:', rinfo(1)
      Else
!       Report failure
        Write (nout,'(/1x,a,i4/)') 'h02ddf returns ifail = ', ifail
      End If

!     Clean up
      Deallocate (b,group,icolb,icolgd,idxfd,irowb,irowgd,iuser,linbl,linbu,   &
        nlnbl,nlnbu,ruser,varbl,varbu,x)

!     Free the handle memory
      ifail = 0
      Call e04rzf(handle,ifail)

    End Program h02ddfe