!   D01RL_A1W_F Example Program Text
!   Mark 30.0 Release. NAG Copyright 2024.

    Module d01rl_a1w_fe_mod

!     .. Use Statements ..
      Use nagad_library, Only: nagad_a1w_w_rtype, sqrt, Assignment (=),        &
                               Operator (==), 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 ..
        iflag = 0
        Do i = 1, nx
          fv(i) = x(i) - 1.0_nag_wp/7.0_nag_wp
          If (fv(i)==0.0_nag_wp) Then
!           singular point has been hit
            iflag = iflag + 1
            ruser(iflag) = x(i)
          Else If (fv(i)<0.0_nag_wp) Then
            fv(i) = -fv(i)
          End If
        End Do
        iuser(1) = iflag
!       signal abort by setting iflag <0
        If (iflag==0) Then
          Do i = 1, nx
            fv(i) = 1.0_nag_wp/sqrt(fv(i))
          End Do
        Else
          iflag = -iflag
        End If
        Return
      End Subroutine f
    End Module d01rl_a1w_fe_mod

    Program d01rl_a1w_fe

!     D01RL_A1W_F Example Main Program

!     .. Use Statements ..
      Use d01rl_a1w_fe_mod, Only: f, nout
      Use, Intrinsic                   :: iso_c_binding, Only: c_null_ptr,     &
                                          c_ptr
      Use nagad_library, Only: d01rl_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)         :: a, abserr, b, epsabs, epsrel, result
      Type (c_ptr)                     :: ad_handle, cpuser
      Real (Kind=nag_wp)               :: da, db
      Integer                          :: ifail, liinfo, lrinfo, maxsub, npts
!     .. Local Arrays ..
      Type (nagad_a1w_w_rtype), Allocatable :: points(:), rinfo(:)
      Type (nagad_a1w_w_rtype)         :: ruser(20)
      Integer, Allocatable             :: iinfo(:)
      Integer                          :: iuser(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..

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

!     Create AD tape
      Call nagad_a1w_ir_create

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

      npts = 1
      maxsub = 20
      liinfo = 2*max(maxsub,npts) + npts + 4
      lrinfo = 4*max(maxsub,npts) + npts + 6

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

      points(1) = 1.0_nag_wp/7.0_nag_wp
      iuser(1) = 0
      ruser(1:20) = 0.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(a)
      Call nagad_a1w_ir_register_variable(b)

!     Evaluate the integral using the AD routine
      ifail = -1
      Call d01rl_a1w_f(ad_handle,f,a,b,npts,points,epsabs,epsrel,maxsub,       &
        result,abserr,rinfo,iinfo,iuser,ruser,cpuser,ifail)

      If (ifail<0) Then
        If (ifail==-1) Then
          Write (nout,99999) 'A user requested exit was issued from F '
        Else
          Write (nout,99999) 'The routine has failed with ifail = ', ifail
        End If
        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,99998) 'points(1)', 'given break-point', points(1)%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,A9,' - ',A32,' = ',F9.4)
99997 Format (1X,A9,' - ',A32,' = ',E9.2)
99996 Format (1X,A9,' - ',A32,' = ',F9.5)
99995 Format (1X,A9,' - ',A32,' = ',I4)

!     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
      da = nagad_a1w_get_derivative(a)
      db = nagad_a1w_get_derivative(b)

      Write (nout,*) ' Derivatives:'
      Write (nout,99994) 'd/da(x) =', da
      Write (nout,99994) 'd/db(x) =', db
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 d01rl_a1w_fe