!   D01RJ_A1W_F Example Program Text
!   Mark 30.1 Release. NAG Copyright 2024.

    Module d01rj_a1w_fe_mod

!     .. Use Statements ..
      Use nagad_library, Only: nagad_a1w_w_rtype, sin, sqrt, Operator (==),    &
                               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
          If (x(i)==1.0_nag_wp) Then
!           An undefined result will be generated.
!           Set iflag to force immediate exit and stoe in iuser
            iflag = -1
            iuser(1) = iflag
          Else
            fv(i) = x(i)*sin(ruser(2)*x(i))/sqrt(1.0_nag_wp-x(i)*x(i)/ruser(1) &
              )
          End If
        End Do
        Return
      End Subroutine f
    End Module d01rj_a1w_fe_mod

    Program d01rj_a1w_fe

!     D01RJ_A1W_F Example Main Program

!     .. Use Statements ..
      Use d01rj_a1w_fe_mod, Only: f, nout
      Use, Intrinsic                   :: iso_c_binding, Only: c_ptr
      Use nagad_library, Only: d01rj_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, x01aaf
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Type (nagad_a1w_w_rtype)         :: a, abserr, b, epsabs, epsrel, result
      Type (c_ptr)                     :: ad_handle, cpuser
      Real (Kind=nag_wp)               :: dr1, dr2, pi
      Integer                          :: ifail, liinfo, lrinfo, maxsub
!     .. Local Arrays ..
      Type (nagad_a1w_w_rtype), Allocatable :: rinfo(:)
      Type (nagad_a1w_w_rtype)         :: ruser(2)
      Integer, Allocatable             :: iinfo(:)
      Integer                          :: iuser(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..

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

!     Create AD tape
      Call nagad_a1w_ir_create

      pi = x01aaf(pi)

!     Initialize AD types
      a = 0.0_nag_wp
      b = 2.0_nag_wp*pi
      epsabs = 0.0_nag_wp
      epsrel = 1.0E-4_nag_wp

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

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

      iuser(1) = 0
      ruser(1) = 4.0_nag_wp*pi*pi
      ruser(2) = 30.0_nag_wp

!     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)

!     Evaluate the integral using the AD routine
      ifail = -1
      Call d01rj_a1w_f(ad_handle,f,a,b,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) 'a     ', 'lower limit of integration', a%value
      Write (nout,99998) 'b     ', 'upper limit of integration', b%value
      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,A6,' - ',A30,' = ',F10.4)
99997 Format (1X,A6,' - ',A30,' = ',E10.2)
99996 Format (1X,A6,' - ',A30,' = ',F10.5)
99995 Format (1X,A6,' - ',A30,' = ',I10)

!     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(1))
      dr2 = nagad_a1w_get_derivative(ruser(2))

      Write (nout,*) ' Derivatives:'
      Write (nout,99994) 'd(result)/druser(1) =', dr1
      Write (nout,99994) 'd(result)/druser(2) =', dr2
99994 Format (1X,A15,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 d01rj_a1w_fe