!   D01ESF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    Module d01esfe_mod

!     D01ESF Example Program Module:
!            User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
      Public                           :: f
    Contains
      Subroutine f(ni,ndim,nx,xtr,nntr,icolzp,irowix,xs,qs,fm,iflag,iuser,     &
        ruser)

!       .. Use Statements ..
        Use nag_library, Only: x06adf
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (In) :: xtr
        Integer, Intent (Inout)        :: iflag
        Integer, Intent (In)           :: ndim, ni, nntr, nx
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: fm(ni,nx)
        Real (Kind=nag_wp), Intent (Inout) :: ruser(*)
        Real (Kind=nag_wp), Intent (In) :: xs(nntr)
        Integer, Intent (In)           :: icolzp(nx+1), irowix(nntr), qs(nntr)
        Integer, Intent (Inout)        :: iuser(*)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: s_ntr, s_tr
        Integer                        :: i, j, logs_hi, logs_lo, s_hi, s_lo,  &
                                          tid
!       .. Intrinsic Procedures ..
        Intrinsic                      :: log, real, sin, sum
!       .. Executable Statements ..

!       For each evaluation point x_i, i = 1, ..., nx, return in fm the
!       computed values of the ni integrands f_j, j = 1, ..., ni defined by

!         fm(j,i) = f_j(x_i)
!                                                         ndim
!                 = sin(j + S(i))*log(S(i)), where S(i) =  Sum  k*x_i(k).
!                                                          k=1

!       Split the S expression into two components, one involving only the
!       'trivial' value xtr:

!                ndim             ndim
!         S(i) =  Sum  (k*xtr)  +  Sum  (k*(x_i(k)-xtr))
!                 k=1              k=1

!                      ndim*(ndim+1)   ndim
!              = xtr * ------------- +  Sum  (k*(x_i(k)-xtr))
!                           2           k=1

!             := s_tr                + s_ntr(i)

!       By definition the summands in the s_ntr(i) term on the right-hand side
!       are zero for those k outside the range of indices defined in irowix.

!       As a demonstration of safely operating with the user arrays iuser and
!       ruser when running in parallel, 'partition' these based on the current
!       thread number. Store some of the s_tr and s_ntr computations in these
!       array sections.

!       The thread number, converted to 1-based numbering.
        tid = x06adf() + 1

        s_lo = iuser(tid)
        s_hi = s_lo + nx - 1
        logs_lo = s_hi + 1
        logs_hi = logs_lo + nx - 1

        If (iflag==0) Then
!         First call: nx=1, no non-trivial dimensions.
!         The constant s_tr can be reused by all subsequent calculations.
          s_tr = 0.5E0_nag_wp*xtr*real(ndim*(ndim+1),kind=nag_wp)
          ruser(1) = s_tr
          ruser(s_lo) = s_tr
          ruser(logs_lo) = log(s_tr)
        Else
!         Calculate S(i) = s_tr + s_ntr(i).
          s_tr = ruser(1)
          Do i = 1, nx
            s_ntr = sum(real(irowix(icolzp(i):icolzp(i+1)-                     &
              1),kind=nag_wp)*(xs(icolzp(i):icolzp(i+1)-1)-xtr))
            ruser(s_lo+i-1) = s_ntr + s_tr
            ruser(logs_lo+i-1) = log(s_ntr+s_tr)
          End Do
        End If

!       Finally we obtain fm(j,:) = sin(j+S(:))*log(S(:))
        Do j = 1, ni
          fm(j,:) = sin(real(j,kind=nag_wp)+ruser(s_lo:s_hi))*                 &
            ruser(logs_lo:logs_hi)
        End Do

        Return
      End Subroutine f
    End Module d01esfe_mod
    Program d01esfe

!     .. Use Statements ..
      Use d01esfe_mod, Only: f
      Use nag_library, Only: d01esf, d01zkf, d01zlf, nag_wp, x06acf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: rvalue
      Integer                          :: ifail, j, liuser, lruser, maxnx,     &
                                          ndim, ni, optype, smpthd
      Character (16)                   :: cvalue
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: dinest(:), errest(:), opts(:),       &
                                          ruser(:)
      Integer, Allocatable             :: iopts(:), iuser(:), ivalid(:),       &
                                          maxdlv(:)
!     .. Executable Statements ..
      Write (nout,*) 'D01ESF Example Program Results'
      Write (nout,*)
      ni = 10
      ndim = 4

      Allocate (iopts(100),opts(100),ivalid(ni),dinest(ni),errest(ni),         &
        maxdlv(ndim))

!     Initialize option arrays.
      ifail = 0
      Call d01zkf('Initialize = D01ESF',iopts,100,opts,100,ifail)

!     Set any required options.
      Call d01zkf('Absolute Tolerance = 0.0',iopts,100,opts,100,ifail)
      Call d01zkf('Relative Tolerance = 1.0e-3',iopts,100,opts,100,ifail)
      Call d01zkf('Maximum Level = 6',iopts,100,opts,100,ifail)
      Call d01zkf('Index Level = 5',iopts,100,opts,100,ifail)

!     Set any required maximum dimension levels.
      maxdlv(:) = 0

!     As a demonstration of safely operating with the user arrays iuser and
!     ruser when running in parallel, we will 'partition' these in the user-
!     supplied function f based on the current thread number.
!     The size of these arrays is a function of Maximum Nx and the maximum
!     allowed number of OpenMP threads.

      ifail = 0
      Call d01zlf('Maximum Nx',maxnx,rvalue,cvalue,optype,iopts,opts,ifail)

      smpthd = x06acf()

      lruser = 1 + 2*maxnx*smpthd
      liuser = smpthd
      Allocate (iuser(liuser),ruser(lruser))

!     iuser stores the partition indices for ruser:
      iuser(1) = 2
      Do j = 2, smpthd
        iuser(j) = iuser(j-1) + 2*maxnx
      End Do

!     Approximate the integrals.
      ifail = -1
      Call d01esf(ni,ndim,f,maxdlv,dinest,errest,ivalid,iopts,opts,iuser,      &
        ruser,ifail)
      Select Case (ifail)
      Case (0,1,2,-1)
!       0: The result returned satisfies the requested accuracy requirements.
!       1, 2: The result returned is inaccurate for at least one integral.
!       -1: Exit was requested by setting iflag negative in f.
!           A result will be returned if at least one call to f was
!           successful.
        Write (nout,99999)
        Do j = 1, ni
          Write (nout,99998) j, dinest(j), errest(j), ivalid(j)
        End Do
      Case Default
!       If internal memory allocation failed consider reducing the options
!       'Maximum Nx' and 'Index Level', or run with fewer threads.
        Write (nout,99997) ifail
      End Select

99999 Format (1X,'Integral # | Estimated value | Error estimate | ',           &
        'Final state of integral')
99998 Format (1X,I11,'|',Es17.5,'|',Es16.5,'|',I8)
99997 Format (1X,'D01ESF exited with IFAIL = ',I8)
    End Program d01esfe