!   D01RM_A1W_F Example Program Text
!   Mark 29.1 Release. NAG Copyright 2023.

    Module d01rm_a1w_fe_mod

!     .. Use Statements ..
      Use nagad_library, Only: nagad_a1w_w_rtype, sqrt, Assignment (=),        &
                               Operator (/), Operator (+), Operator (*)
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: f
!     .. Parameters ..
      Integer, Parameter, Public       :: nout = 6
    Contains
      Subroutine f(ad_handle,x,nx,fv,iflag,iuser,ruser,cpuser)
!       .. Use Statements ..
        Use, Intrinsic                 :: iso_c_binding, Only: c_ptr
!       .. Scalar Arguments ..
        Type (c_ptr), Intent (Inout)   :: ad_handle
        Type (c_ptr), Intent (In)      :: cpuser
        Integer, Intent (Inout)        :: iflag
        Integer, Intent (In)           :: nx
!       .. Array Arguments ..
        Type (nagad_a1w_w_rtype), Intent (Out) :: fv(nx)
        Type (nagad_a1w_w_rtype), Intent (Inout) :: ruser(*)
        Type (nagad_a1w_w_rtype), Intent (In) :: x(nx)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Integer                        :: i
!       .. Executable Statements ..
        Do i = 1, nx
          fv(i) = 1.0_nag_wp/((x(i)+ruser(21))*sqrt(ruser(22)*x(i)))
        End Do
        Return
      End Subroutine f
    End Module d01rm_a1w_fe_mod

    Program d01rm_a1w_fe

!     D01RM_A1W_F Example Main Program

!     .. Use Statements ..
      Use d01rm_a1w_fe_mod, Only: f, nout
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nagad_library, Only: d01rm_a1w_f, nagad_a1w_get_derivative,          &
                               nagad_a1w_inc_derivative,                       &
                               nagad_a1w_ir_create => x10za_a1w_f,             &
                               nagad_a1w_ir_interpret_adjoint_sparse,          &
                               nagad_a1w_ir_register_variable,                 &
                               nagad_a1w_ir_remove, nagad_a1w_w_rtype,         &
                               x10aa_a1w_f, x10ab_a1w_f, Assignment (=)
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (nagad_a1w_w_rtype)         :: abserr, bound, epsabs, epsrel,       &
                                          result
      Type (c_ptr)                     :: ad_handle, cpuser
      Real (Kind=nag_wp)               :: dr1, dr2
      Integer                          :: ifail, inf, liinfo, lrinfo, maxsub
!     .. Local Arrays ..
      Type (nagad_a1w_w_rtype), Allocatable :: rinfo(:)
      Type (nagad_a1w_w_rtype)         :: ruser(22)
      Integer, Allocatable             :: iinfo(:)
      Integer                          :: iuser(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..

      Write (nout,*) 'D01RM_A1W_F Example Program Results'

!     Create AD tape
      Call nagad_a1w_ir_create

!     Initialize AD types
      bound = 0.0_nag_wp
      epsabs = 0.0_nag_wp
      epsrel = 1.0E-4_nag_wp

      inf = 1
      maxsub = 20
      lrinfo = 4*maxsub
      liinfo = max(maxsub,4)

      Allocate (rinfo(lrinfo),iinfo(liinfo))

      iuser(1) = 0
      ruser(1:20) = 0.0_nag_wp
      ruser(21) = 1.0_nag_wp
      ruser(22) = 1.0_nag_wp
      cpuser = c_null_ptr

!     Create AD configuration data object
      ifail = 0
      Call x10aa_a1w_f(ad_handle,ifail)

!     Register variables to differentiate w.r.t.
      Call nagad_a1w_ir_register_variable(ruser(21))
      Call nagad_a1w_ir_register_variable(ruser(22))

!     Evaluate the integral using the AD routine
      ifail = -1
      Call d01rm_a1w_f(ad_handle,f,bound,inf,epsabs,epsrel,maxsub,result,      &
        abserr,rinfo,iinfo,iuser,ruser,cpuser,ifail)

      If (ifail<0) Then
        Write (nout,99999) 'The routine has failed with ifail = ', ifail
        Go To 100
99999   Format (1X,A,I0)
      End If
!     Print inputs and primal outputs
      Write (nout,*)
      Write (nout,99998) 'bound    ', 'lower limit of integration',            &
        bound%value
      Write (nout,99994) 'Inf      ', 'upper limit of integration', 'Infinity'
      Write (nout,99997) 'epsabs   ', 'absolute accuracy requested',           &
        epsabs%value
      Write (nout,99997) 'epsrel   ', 'relative accuracy requested',           &
        epsrel%value
      Write (nout,99995) 'maxsub   ', 'max number of subintervals', maxsub
      Write (nout,*)
      If (ifail>=0) Then
        Write (nout,99996) 'result   ', 'approximation to the integral',       &
          result%value
        Write (nout,99997) 'abserr   ', 'estimate of the absolute error',      &
          abserr%value
        Write (nout,99995) 'iinfo(1) ', 'number of subintervals used',         &
          iinfo(1)
      End If
99998 Format (1X,A9,' - ',A32,' = ',F9.4)
99997 Format (1X,A9,' - ',A32,' = ',E9.2)
99996 Format (1X,A9,' - ',A32,' = ',F9.5)
99995 Format (1X,A9,' - ',A32,' = ',I4)
99994 Format (1X,A9,' - ',A32,' = ',A8)

!     Setup evaluation of derivatives via adjoints
      Call nagad_a1w_inc_derivative(result,1.0E0_nag_wp)
      Call nagad_a1w_ir_interpret_adjoint_sparse(ifail)

      Write (nout,*)
      Write (nout,*) ' Derivatives calculated: First order adjoints'
      Write (nout,*) ' Computational mode    : algorithmic'

!     Get derivatives
      dr1 = nagad_a1w_get_derivative(ruser(21))
      dr2 = nagad_a1w_get_derivative(ruser(22))

      Write (nout,*) ' Derivatives:'
      Write (nout,99993) 'd(result)/druser(21) =', dr1
      Write (nout,99993) 'd(result)/druser(22) =', dr2
99993 Format (1X,A22,1X,F10.5)

100   Continue
!     Remove computational data object and tape
      Call x10ab_a1w_f(ad_handle,ifail)
      Call nagad_a1w_ir_remove

    End Program d01rm_a1w_fe