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

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

    Module e04usae_mod

!     E04USA 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
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter, Public :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter, Public :: zero = 0.0_nag_wp
      Integer, Parameter, Public       :: inc1 = 1, lcwsav = 1, liwsav = 610,  &
                                          llwsav = 120, lrwsav = 475, nin = 5, &
                                          nout = 6
    Contains
      Subroutine objfun(mode,m,n,ldfj,needfi,x,f,fjac,nstate,iuser,ruser)
!       Routine to evaluate the subfunctions and their 1st derivatives.

!       .. Parameters ..
        Real (Kind=nag_wp), Parameter  :: a(44) = (/8.0_nag_wp,8.0_nag_wp,     &
                                          10.0_nag_wp,10.0_nag_wp,10.0_nag_wp, &
                                          10.0_nag_wp,12.0_nag_wp,12.0_nag_wp, &
                                          12.0_nag_wp,12.0_nag_wp,14.0_nag_wp, &
                                          14.0_nag_wp,14.0_nag_wp,16.0_nag_wp, &
                                          16.0_nag_wp,16.0_nag_wp,18.0_nag_wp, &
                                          18.0_nag_wp,20.0_nag_wp,20.0_nag_wp, &
                                          20.0_nag_wp,22.0_nag_wp,22.0_nag_wp, &
                                          22.0_nag_wp,24.0_nag_wp,24.0_nag_wp, &
                                          24.0_nag_wp,26.0_nag_wp,26.0_nag_wp, &
                                          26.0_nag_wp,28.0_nag_wp,28.0_nag_wp, &
                                          30.0_nag_wp,30.0_nag_wp,30.0_nag_wp, &
                                          32.0_nag_wp,32.0_nag_wp,34.0_nag_wp, &
                                          36.0_nag_wp,36.0_nag_wp,38.0_nag_wp, &
                                          38.0_nag_wp,40.0_nag_wp,             &
                                          42.0_nag_wp/)
!       .. Scalar Arguments ..
        Integer, Intent (In)           :: ldfj, m, n, needfi, nstate
        Integer, Intent (Inout)        :: mode
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: f(m)
        Real (Kind=nag_wp), Intent (Inout) :: fjac(ldfj,n), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: ai, temp, x1, x2
        Integer                        :: i
        Logical                        :: mode02, mode12
!       .. Intrinsic Procedures ..
        Intrinsic                      :: exp
!       .. Executable Statements ..
        x1 = x(1)
        x2 = x(2)

        mode02 = (mode==0 .Or. mode==2)
        mode12 = (mode==1 .Or. mode==2)

loop:   Do i = 1, m

          If (needfi==i) Then
            f(i) = x1 + (0.49_nag_wp-x1)*exp(-x2*(a(i)-8.0_nag_wp))
            Exit loop
          End If

          ai = a(i)
          temp = exp(-x2*(ai-8.0_nag_wp))

          If (mode02) Then
            f(i) = x1 + (0.49_nag_wp-x1)*temp
          End If

          If (mode12) Then
            fjac(i,1) = one - temp
            fjac(i,2) = -(0.49_nag_wp-x1)*(ai-8.0_nag_wp)*temp
          End If

        End Do loop

        Return

      End Subroutine objfun
      Subroutine confun(mode,ncnln,n,ldcj,needc,x,c,cjac,nstate,iuser,ruser)
!       Routine to evaluate the nonlinear constraint and its 1st
!       derivatives.

!       .. Scalar Arguments ..
        Integer, Intent (In)           :: ldcj, n, ncnln, nstate
        Integer, Intent (Inout)        :: mode
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: c(ncnln)
        Real (Kind=nag_wp), Intent (Inout) :: cjac(ldcj,n), ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
        Integer, Intent (In)           :: needc(ncnln)
!       .. Executable Statements ..
        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.2).

          cjac(1:ncnln,1:n) = zero
        End If

        If (needc(1)>0) Then

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

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

        End If

        Return

      End Subroutine confun
    End Module e04usae_mod
    Program e04usae

!     E04USA Example Main Program

!     .. Use Statements ..
      Use e04usae_mod, Only: confun, inc1, lcwsav, liwsav, llwsav, lrwsav,     &
                             nin, nout, objfun, one, zero
      Use nag_library, Only: dgemv, e04usa, e04wbf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: objf
      Integer                          :: i, ifail, iter, j, lda, ldcj, ldfj,  &
                                          ldr, liwork, lwork, m, n, nclin,     &
                                          ncnln, sda, sdcjac
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), bl(:), bu(:), c(:),          &
                                          cjac(:,:), clamda(:), f(:),          &
                                          fjac(:,:), r(:,:), work(:), x(:),    &
                                          y(:)
      Real (Kind=nag_wp)               :: rwsav(lrwsav), user(1)
      Integer, Allocatable             :: istate(:), iwork(:)
      Integer                          :: iuser(1), iwsav(liwsav)
      Logical                          :: lwsav(llwsav)
      Character (80)                   :: cwsav(lcwsav)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..
      Write (nout,*) 'E04USA Example Program Results'

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

      Read (nin,*) m, n
      Read (nin,*) nclin, ncnln
      liwork = 3*n + nclin + 2*ncnln
      lda = max(1,nclin)

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

      ldcj = max(1,ncnln)

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

      ldfj = m
      ldr = n

      If (ncnln==0 .And. nclin>0) Then
        lwork = 2*n**2 + 20*n + 11*nclin + m*(n+3)
      Else If (ncnln>0 .And. nclin>=0) Then
        lwork = 2*n**2 + n*nclin + 2*n*ncnln + 20*n + 11*nclin + 21*ncnln +    &
          m*(n+3)
      Else
        lwork = 20*n + m*(n+3)
      End If

      Allocate (istate(n+nclin+ncnln),iwork(liwork),a(lda,sda),                &
        bl(n+nclin+ncnln),bu(n+nclin+ncnln),y(m),c(max(1,                      &
        ncnln)),cjac(ldcj,sdcjac),f(m),fjac(ldfj,n),clamda(n+nclin+ncnln),     &
        r(ldr,n),x(n),work(lwork))

      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))
      Read (nin,*) x(1:n)

!     Initialise E04USA

      ifail = 0
      Call e04wbf('E04USA',cwsav,lcwsav,lwsav,llwsav,iwsav,liwsav,rwsav,       &
        lrwsav,ifail)

!     Solve the problem

      ifail = -1
      Call e04usa(m,n,nclin,ncnln,lda,ldcj,ldfj,ldr,a,bl,bu,y,confun,objfun,   &
        iter,istate,c,cjac,f,fjac,clamda,objf,r,x,iwork,liwork,work,lwork,     &
        iuser,user,lwsav,iwsav,rwsav,ifail)

      Select Case (ifail)
      Case (0:8)
        Write (nout,*)
        Write (nout,99999)
        Write (nout,*)

        Do i = 1, n
          Write (nout,99998) i, istate(i), x(i), clamda(i)
        End Do

        If (nclin>0) Then

!         A*x --> work
!         The NAG name equivalent of dgemv is f06paf
          Call dgemv('N',nclin,n,one,a,lda,x,inc1,zero,work,inc1)

          Write (nout,*)
          Write (nout,*)
          Write (nout,99997)
          Write (nout,*)

          Do i = n + 1, n + nclin
            j = i - n
            Write (nout,99996) j, istate(i), work(j), clamda(i)
          End Do

        End If

        If (ncnln>0) Then
          Write (nout,*)
          Write (nout,*)
          Write (nout,99995)
          Write (nout,*)

          Do i = n + nclin + 1, n + nclin + ncnln
            j = i - n - nclin
            Write (nout,99994) j, istate(i), c(j), clamda(i)
          End Do

        End If

        Write (nout,*)
        Write (nout,*)
        Write (nout,99993) objf
      End Select

99999 Format (1X,'Varbl',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99998 Format (1X,'V',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99997 Format (1X,'L Con',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99996 Format (1X,'L',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99995 Format (1X,'N Con',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99994 Format (1X,'N',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99993 Format (1X,'Final objective value = ',G15.7)
    End Program e04usae