!   E04XAF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    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)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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