!   E04XAF Example Program Text
!   Mark 24 Release. NAG Copyright 2012.

    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
!     .. Parameters ..
      Integer, Parameter                   :: n = 4, nout = 6
      Integer, Parameter                   :: lhes = n
      Integer, Parameter                   :: 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 nag_library, Only: e04xaf, nag_wp
      Use e04xafe_mod, Only: lhes, lwork, n, nout, objfun
!     .. 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 (80)                       :: 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) = (/3.0E0_nag_wp,-1.0E0_nag_wp,0.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.4)
99997 Format (4(1X,1P,E12.4))
    End Program e04xafe
