! 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