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

    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
!      .. Parameters ..
       INTEGER, PARAMETER                  :: lw = 800, nout = 6, npts = 1
       INTEGER, PARAMETER                  :: liw = lw/2
    CONTAINS
       FUNCTION f(x)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. 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 Functions ..
          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 nag_library, ONLY : d01alf, nag_wp
       USE d01alfe_mod, ONLY : f, liw, lw, nout, npts
!      .. 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