!   E04XAF Example Program Text
!   Mark 28.3 Release. NAG Copyright 2022.

    Module e04xafe_mod

!     E04XAF 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                           :: objfun
!     .. Parameters ..
      Integer, Parameter, Public       :: n = 4, nout = 6
      Integer, Parameter, Public       :: lhes = n
      Integer, Parameter, Public       :: lwork = n*n + n
    Contains
      Subroutine objfun(mode,n,x,objf,objgrd,nstate,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: objf
        Integer, Intent (Inout)        :: mode
        Integer, Intent (In)           :: n, nstate
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: objgrd(n)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: x(n)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: a, b, c, d
!       .. Executable Statements ..
        a = x(1) + 10.0E0_nag_wp*x(2)
        b = x(3) - x(4)
        c = x(2) - 2.0E0_nag_wp*x(3)
        d = x(1) - x(4)
        objf = a**2 + 5.0E0_nag_wp*b**2 + c**4 + 10.0E0_nag_wp*d**4

        If (mode==1) Then
          objgrd(1) = 4.0E1_nag_wp*x(1)**3 + 2.0E0_nag_wp*x(1) -               &
            1.2E2_nag_wp*x(4)*x(1)**2 + 1.2E2_nag_wp*x(1)*x(4)**2 +            &
            2.0E1_nag_wp*x(2) - 4.0E1_nag_wp*x(4)**3
          objgrd(2) = 2.0E2_nag_wp*x(2) + 2.0E1_nag_wp*x(1) +                  &
            4.0E0_nag_wp*x(2)**3 + 4.8E1_nag_wp*x(2)*x(3)**2 -                 &
            2.4E1_nag_wp*x(3)*x(2)**2 - 32.0E0_nag_wp*x(3)**3
          objgrd(3) = 1.0E1_nag_wp*x(3) - 1.0E1_nag_wp*x(4) -                  &
            8.0E0_nag_wp*x(2)**3 + 4.8E1_nag_wp*x(3)*x(2)**2 -                 &
            9.6E1_nag_wp*x(2)*x(3)**2 + 6.4E1_nag_wp*x(3)**3
          objgrd(4) = 1.0E1_nag_wp*x(4) - 1.0E1_nag_wp*x(3) -                  &
            4.0E1_nag_wp*x(1)**3 + 1.2E2_nag_wp*x(4)*x(1)**2 -                 &
            1.2E2_nag_wp*x(1)*x(4)**2 + 4.0E1_nag_wp*x(4)**3
        End If

        Return

      End Subroutine objfun
    End Module e04xafe_mod
    Program e04xafe

!     E04XAF Example Main Program

!     .. Use Statements ..
      Use e04xafe_mod, Only: lhes, lwork, n, nout, objfun
      Use nag_library, Only: e04xaf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: epsrf, objf
      Integer                          :: i, ifail, imode, iwarn, mode, msglvl
!     .. Local Arrays ..
      Real (Kind=nag_wp)               :: hcntrl(n), hesian(lhes,n), hforw(n), &
                                          objgrd(n), user(1), work(lwork),     &
                                          x(n)
      Integer                          :: info(n), iuser(1)
      Character (78)                   :: rc(3)
!     .. Executable Statements ..
      Write (nout,*) 'E04XAF Example Program Results'

      msglvl = 0

!     Set the point at which the derivatives are to be estimated.

      x(1:n) = (/2.0E0_nag_wp,-1.0E0_nag_wp,1.0E0_nag_wp,1.0E0_nag_wp/)

      rc(1) = 'Find gradients and Hessian diagonals given function only'
      rc(2) = 'Find Hessian matrix given function and gradients'
      rc(3) = 'Find gradients and Hessian matrix given function only'

!     Take default value of EPSRF.

      epsrf = -1.0E0_nag_wp

!     Illustrate the different values of MODE.

loop: Do imode = 0, 2
        mode = imode

!       Set HFORW(I) = -1.0 so that E04XAF computes the initial trial
!       interval.

        hforw(1:n) = -1.0E0_nag_wp

        ifail = -1
        Call e04xaf(msglvl,n,epsrf,x,mode,objfun,lhes,hforw,objf,objgrd,       &
          hcntrl,hesian,iwarn,work,iuser,user,info,ifail)

        Select Case (ifail)
        Case (0,2)
          Write (nout,99999) rc(mode+1), mode
          Write (nout,99998) 'Function value is ', objf

          If (mode==1) Then
            Write (nout,*) 'Gradient vector is'
            Write (nout,99997) objgrd(1:n)
          Else
            Write (nout,*) 'Estimated gradient vector is'
            Write (nout,99997) objgrd(1:n)
          End If

          If (mode==0) Then
            Write (nout,*) 'Estimated Hessian matrix diagonal is'
            Write (nout,99997) hesian(1:n,1)
          Else
            Write (nout,*) 'Estimated Hessian matrix (machine dependent) is'
            Write (nout,99997)(hesian(i,1:n),i=1,n)
          End If

        Case Default
          Exit loop
        End Select

      End Do loop

99999 Format (1X,/,1X,A,/,1X,'( i.e. MODE =',I2,' ).')
99998 Format (1X,A,1P,E12.3)
99997 Format (4(1X,1P,E12.3))
    End Program e04xafe