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

NAG FL Interface Introduction
Example description
!   E05USF Example Program Text
!   Mark 30.3 Release. nAG Copyright 2024.

    Module e05usfe_mod

!     E05USF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: confun, objfun, start
    Contains
      Subroutine objfun(mode,m,n,ldfjsl,needfi,x,f,fjsl,nstate,iuser,ruser)
!       Evaluates the subfunctions and their 1st derivatives.

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: a(44) = (/8._nag_wp,8._nag_wp,       &
                                          10._nag_wp,10._nag_wp,10._nag_wp,    &
                                          10._nag_wp,12._nag_wp,12._nag_wp,    &
                                          12._nag_wp,12._nag_wp,14._nag_wp,    &
                                          14._nag_wp,14._nag_wp,16._nag_wp,    &
                                          16._nag_wp,16._nag_wp,18._nag_wp,    &
                                          18._nag_wp,20._nag_wp,20._nag_wp,    &
                                          20._nag_wp,22._nag_wp,22._nag_wp,    &
                                          22._nag_wp,24._nag_wp,24._nag_wp,    &
                                          24._nag_wp,26._nag_wp,26._nag_wp,    &
                                          26._nag_wp,28._nag_wp,28._nag_wp,    &
                                          30._nag_wp,30._nag_wp,30._nag_wp,    &
                                          32._nag_wp,32._nag_wp,34._nag_wp,    &
                                          36._nag_wp,36._nag_wp,38._nag_wp,    &
                                          38._nag_wp,40._nag_wp,42._nag_wp/)
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: ldfjsl, m, n, needfi, nstate
        Integer, Intent (Inout)        :: mode
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(m)
        Real (Kind=nag_wp), Intent (Inout) :: fjsl(ldfjsl,*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: temp
        Integer                        :: i
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..

!       This is a two-dimensional objective function.
!       As an example of using the mode mechanism,
!       terminate if any other problem size is supplied.

        If (n/=2) Then
          mode = -1
        End If

        If (nstate==1) Then
!         This is the first call.
!         Take any special action here if desired.
        End If

        If (mode==0 .And. needfi>0) Then
          f(needfi) = x(1) + (0.49_nag_wp-x(1))*exp(-x(2)*(a(needfi)-          &
            8.0_nag_wp))
        Else
          Do i = 1, m
            temp = exp(-x(2)*(a(i)-8._nag_wp))

            If (mode==0 .Or. mode==2) Then
              f(i) = x(1) + (0.49_nag_wp-x(1))*temp
            End If

            If (mode==1 .Or. mode==2) Then
              fjsl(i,1) = 1._nag_wp - temp
              fjsl(i,2) = -(0.49_nag_wp-x(1))*(a(i)-8._nag_wp)*temp
            End If

          End Do
        End If

        Return
      End Subroutine objfun
      Subroutine confun(mode,ncnln,n,ldcjsl,needc,x,c,cjsl,nstate,iuser,ruser)
!       Evaluates the nonlinear constraints and their 1st derivatives.

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: ldcjsl, n, ncnln, nstate
        Integer, Intent (Inout)        :: mode
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: c(ncnln)
        Real (Kind=nag_wp), Intent (Inout) :: cjsl(ldcjsl,*), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: needc(ncnln)
!       .. Executable Statements ..

!       This problem has only one constraint.
!       As an example of using the mode mechanism,
!       terminate if any other size is supplied.

        If (ncnln/=1) Then
          mode = -1
        End If

        If (nstate==1) Then

!         First call to CONFUN.  Set all Jacobian elements to zero.
!         Note that this will only work when 'Derivative Level = 3'
!         (the default; see Section 11.1 of the E04USA document).

          cjsl(1:ncnln,1:n) = 0._nag_wp
        End If

        If (needc(1)>0) Then

          If (mode==0 .Or. mode==2) Then
            c(1) = -0.09_nag_wp - x(1)*x(2) + 0.49_nag_wp*x(2)
          End If

          If (mode==1 .Or. mode==2) Then
            cjsl(1,1) = -x(2)
            cjsl(1,2) = -x(1) + 0.49_nag_wp
          End If

        End If

        Return
      End Subroutine confun
      Subroutine start(npts,quas,n,repeat1,bl,bu,iuser,ruser,mode)

!       Sets the initial points.
!       A typical user-defined start procedure.

!       .. Use Statements ..
        Use nag_library, Only: g05kgf, g05saf
!       .. Scalar Arguments ..
        Integer, Intent (Inout)        :: mode
        Integer, Intent (In)           :: n, npts
        Logical, Intent (In)           :: repeat1
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: bl(n), bu(n)
        Real (Kind=nag_wp), Intent (Inout) :: quas(n,npts), ruser(*)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Integer                        :: genid, i, ifail, lstate, subid
!       .. Local Arrays ..
        Integer, Allocatable           :: state(:)
!       .. Executable Statements ..
!       quas is pre-assigned to zero.
        If (repeat1) Then
          quas(1,1) = 0.4_nag_wp
          quas(2,2) = 1._nag_wp
        Else
!         Generate a non-repeatable spread of points between bl and bu.
          genid = 2
          subid = 53
          lstate = -1
          Allocate (state(lstate))
          ifail = 0
          Call g05kgf(genid,subid,state,lstate,ifail)
          Deallocate (state)
          Allocate (state(lstate))
          ifail = 0
          Call g05kgf(genid,subid,state,lstate,ifail)
          Do i = 1, npts
            ifail = 0
            Call g05saf(n,state,quas(1,i),ifail)
            quas(1:n,i) = bl(1:n) + (bu(1:n)-bl(1:n))*quas(1:n,i)
          End Do
          Deallocate (state)
        End If
!       Set mode negative to terminate execution for any reason.
        mode = 0
        Return
      End Subroutine start
    End Module e05usfe_mod
    Program e05usfe

!     E05USF Example Main Program

!     .. Use Statements ..
      Use e05usfe_mod, Only: confun, objfun, start
      Use nag_library, Only: dgemv, e05usf, e05zkf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: liopts = 740, lopts = 485, m = 44,   &
                                          n = 2, nb = 1, nclin = 1, ncnln = 1, &
                                          nin = 5, nout = 6, npts = 3
      Integer, Parameter               :: sdfjac = n
      Integer, Parameter               :: lda = nclin
      Integer, Parameter               :: ldc = ncnln
      Integer, Parameter               :: ldcjac = ncnln
      Integer, Parameter               :: ldclda = n + nclin + ncnln
      Integer, Parameter               :: ldfjac = m
      Integer, Parameter               :: ldx = n
      Integer, Parameter               :: listat = n + nclin + ncnln
      Logical, Parameter               :: repeat1 = .True.
!     .. Local Scalars ..
      Integer                          :: i, ifail, j, k, l, sda, sdcjac
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), bl(:), bu(:), c(:,:),        &
                                          cjac(:,:,:), clamda(:,:), f(:,:),    &
                                          fjac(:,:,:), work(:), x(:,:), y(:)
      Real (Kind=nag_wp)               :: objf(nb), opts(lopts), ruser(1)
      Integer                          :: info(nb), iopts(liopts), iter(nb),   &
                                          iuser(1)
      Integer, Allocatable             :: istate(:,:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..
      Write (nout,*) 'E05USF Example Program Results'
      Flush (nout)

!     Skip heading in data file.
      Read (nin,*)

      If (nclin>0) Then
        sda = n
      Else
        sda = 1
      End If

      If (ncnln>0) Then
        sdcjac = n
      Else
        sdcjac = 0
      End If

      Allocate (a(lda,sda),bl(n+nclin+ncnln),bu(n+nclin+ncnln),y(m),c(ldc,nb), &
        cjac(ldcjac,sdcjac,nb),f(m,nb),fjac(ldfjac,sdfjac,nb),                 &
        clamda(ldclda,nb),istate(listat,nb),x(ldx,nb),work(max(1,nclin)))

      If (nclin>0) Then
        Read (nin,*)(a(i,1:sda),i=1,nclin)
      End If

      Read (nin,*) y(1:m)
      Read (nin,*) bl(1:(n+nclin+ncnln))
      Read (nin,*) bu(1:(n+nclin+ncnln))

!     Initialize the solver.

      ifail = 0
      Call e05zkf('Initialize = E05USF',iopts,liopts,opts,lopts,ifail)

!     Solve the problem.

      ifail = -1
      Call e05usf(m,n,nclin,ncnln,a,lda,bl,bu,y,confun,objfun,npts,x,ldx,      &
        start,repeat1,nb,objf,f,fjac,ldfjac,sdfjac,iter,c,ldc,cjac,ldcjac,     &
        sdcjac,clamda,ldclda,istate,listat,iopts,opts,iuser,ruser,info,ifail)

      Select Case (ifail)
      Case (0)
        l = nb
      Case (8)
        l = info(nb)
        Write (nout,99999) iter(nb)
      Case Default
        Go To 100
      End Select

loop: Do i = 1, l
        Write (nout,99998) i
        Write (nout,99997) info(i)
        Write (nout,99996) 'Varbl'
        Do j = 1, n
          Write (nout,99995) 'V', j, istate(j,i), x(j,i), clamda(j,i)
        End Do
        If (nclin>0) Then
          Write (nout,99996) 'L Con'

!         Below is a call to the level 2 BLAS routine DGEMV.
!         This performs the matrix vector multiplication A*X
!         (linear constraint values) and puts the result in
!         the first NCLIN locations of WORK.

          Call dgemv('N',nclin,n,1.0_nag_wp,a,lda,x(1,i),1,0.0_nag_wp,work,1)

          Do k = n + 1, n + nclin
            j = k - n
            Write (nout,99995) 'L', j, istate(k,i), work(j), clamda(k,i)
          End Do
        End If
        If (ncnln>0) Then
          Write (nout,99996) 'N Con'
          Do k = n + nclin + 1, n + nclin + ncnln
            j = k - n - nclin
            Write (nout,99995) 'N', j, istate(k,i), c(j,i), clamda(k,i)
          End Do
        End If
        Write (nout,99994) objf(i)
        Write (nout,99993)
        Write (nout,99992)(clamda(k,i),k=1,n+nclin+ncnln)

        If (l==1) Then
          Exit loop
        End If

        Write (nout,*)

        Write (nout,*)                                                         &
          ' ------------------------------------------------------ '

      End Do loop

100   Continue

99999 Format (1X,I20,'starting points converged')
99998 Format (/,1X,'Solution number',I16)
99997 Format (/,1X,'Local minimization exited with code',I5)
99996 Format (/,1X,A,2X,'Istate',3X,'Value',9X,'Lagr Mult',/)
99995 Format (1X,A,2(1X,I3),4X,F12.4,2X,F12.4)
99994 Format (/,1X,'Final objective value = ',1X,F12.4)
99993 Format (/,1X,'QP multipliers')
99992 Format (1X,F12.4)
    End Program e05usfe