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

    Module d01alfe_mod

!     D01ALF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: f
!     .. Parameters ..
      Integer, Parameter, Public       :: lw = 800, nout = 6, npts = 1
      Integer, Parameter, Public       :: liw = lw/2
    Contains
      Function f(x)

!       .. Function Return Value ..
        Real (Kind=nag_wp)             :: f
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: x
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: a
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs
!       .. Executable Statements ..
        a = abs(x-1.0E0_nag_wp/7.0E0_nag_wp)

        If (a/=0.0E0_nag_wp) Then
          f = a**(-0.5E0_nag_wp)
        Else
          f = 0.0E0_nag_wp
        End If

        Return

      End Function f
    End Module d01alfe_mod
    Program d01alfe

!     D01ALF Example Main Program

!     .. Use Statements ..
      Use d01alfe_mod, Only: f, liw, lw, nout, npts
      Use nag_library, Only: d01alf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: a, abserr, b, epsabs, epsrel, result
      Integer                          :: ifail
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: points(:), w(:)
      Integer, Allocatable             :: iw(:)
!     .. Executable Statements ..
      Write (nout,*) 'D01ALF Example Program Results'

      Allocate (points(npts),w(lw),iw(liw))

      epsabs = 0.0E0_nag_wp
      epsrel = 1.0E-03_nag_wp
      a = 0.0E0_nag_wp
      b = 1.0E0_nag_wp
      points(1) = 1.0E0_nag_wp/7.0E0_nag_wp

      ifail = -1
      Call d01alf(f,a,b,npts,points,epsabs,epsrel,result,abserr,w,lw,iw,liw,   &
        ifail)

      If (ifail>=0) Then
        Write (nout,*)
        Write (nout,99999) 'A     ', 'lower limit of integration', a
        Write (nout,99999) 'B     ', 'upper limit of integration', b
        Write (nout,99998) 'EPSABS', 'absolute accuracy requested', epsabs
        Write (nout,99998) 'EPSREL', 'relative accuracy requested', epsrel
        Write (nout,99995) 'POINTS(1)', 'given break-point', points(1)
      End If

      If (ifail>=0 .And. ifail<=5) Then
        Write (nout,*)
        Write (nout,99997) 'RESULT', 'approximation to the integral', result
        Write (nout,99998) 'ABSERR', 'estimate of the absolute error', abserr
        Write (nout,99996) 'IW(1) ', 'number of subintervals used', iw(1)
      End If

99999 Format (1X,A6,'    - ',A32,' = ',F10.4)
99998 Format (1X,A6,'    - ',A32,' = ',E9.2)
99997 Format (1X,A6,'    - ',A32,' = ',F9.5)
99996 Format (1X,A6,'    - ',A32,' = ',I4)
99995 Format (1X,A9,' - ',A32,' = ',F10.4)

    End Program d01alfe