! D01RAF Example Program Text
! Mark 26.2 Release. NAG Copyright 2017.
Module d01rafe_mod
! .. Use Statements ..
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Accessibility Statements ..
Private
Public :: display_integration_details, &
display_option
! .. Parameters ..
Integer, Parameter, Public :: nout = 6
Logical, Parameter, Public :: disp_integration_info = .True.
Contains
Subroutine display_integration_details(ni,iopts,opts,icom,licom,com, &
lcom)
! .. Use Statements ..
Use nag_library, Only: d01zlf
! .. Implicit None Statement ..
Implicit None
! .. Scalar Arguments ..
Integer, Intent (In) :: lcom, licom, ni
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: com(lcom), opts(100)
Integer, Intent (In) :: icom(licom), iopts(100)
! .. Local Scalars ..
Real (Kind=nag_wp) :: rvalue
Integer :: dsp, esp, evalsp, fcp, ifail, index, &
ldi, ldr, nsdiv, nseg, optype, &
sinfoip, sinforp
Character (16) :: cvalue
! .. Executable Statements ..
! Request communication array indices.
ifail = 0
Call d01zlf('Index nseg',index,rvalue,cvalue,optype,iopts,opts,ifail)
nseg = icom(index)
Call d01zlf('Index nsdiv',index,rvalue,cvalue,optype,iopts,opts,ifail)
nsdiv = icom(index)
Call d01zlf('Index ldi',index,rvalue,cvalue,optype,iopts,opts,ifail)
ldi = icom(index)
Call d01zlf('Index ldr',index,rvalue,cvalue,optype,iopts,opts,ifail)
ldr = icom(index)
Call d01zlf('Index fcp',index,rvalue,cvalue,optype,iopts,opts,ifail)
fcp = icom(index)
Call d01zlf('Index evalsp',index,rvalue,cvalue,optype,iopts,opts, &
ifail)
evalsp = icom(index)
Call d01zlf('Index sinfoip',index,rvalue,cvalue,optype,iopts,opts, &
ifail)
sinfoip = icom(index)
Call d01zlf('Index dsp',index,rvalue,cvalue,optype,iopts,opts,ifail)
dsp = icom(index)
Call d01zlf('Index esp',index,rvalue,cvalue,optype,iopts,opts,ifail)
esp = icom(index)
Call d01zlf('Index sinforp',index,rvalue,cvalue,optype,iopts,opts, &
ifail)
sinforp = icom(index)
Write (nout,*)
Write (nout,99999)
Write (nout,99998) ni, nseg, nsdiv
Call display_subdivision_strategy(ni,nseg,icom(fcp),icom(sinfoip), &
icom(evalsp),ldi,com(sinforp),com(dsp),com(esp),ldr)
Return
99999 Format (' Information on integration: ')
99998 Format (' NI = ',I3,', nseg = ',I3,', nsdiv = ',I3,'.')
End Subroutine display_integration_details
Subroutine display_subdivision_strategy(ni,nseg,fcount,sinfoi,evals,ldi, &
sinfor,ds,es,ldr)
! .. Scalar Arguments ..
Integer, Intent (In) :: ldi, ldr, ni, nseg
! .. Array Arguments ..
Real (Kind=nag_wp), Intent (In) :: ds(ldr,*), es(ldr,*), sinfor(ldr,*)
Integer, Intent (In) :: evals(ldi,*), fcount(ni), &
sinfoi(ldi,*)
! .. Local Scalars ..
Real (Kind=nag_wp) :: lbnd, ubnd
Integer :: child1, child2, j, k, level, parent, &
sid
! .. Executable Statements ..
! Display information on individual segments.
Do j = 1, ni
Write (nout,99991) j, fcount(j)
End Do
Write (nout,*)
Write (nout,99992)
Do k = 1, nseg
Write (nout,*)
sid = sinfoi(1,k)
parent = sinfoi(2,k)
child1 = sinfoi(3,k)
child2 = sinfoi(4,k)
level = sinfoi(5,k)
lbnd = sinfor(1,k)
ubnd = sinfor(2,k)
Write (nout,99999) k
Write (nout,99998) sid, parent, level
If (child1>0) Then
Write (nout,99997) child1, child2
End If
Write (nout,99996) lbnd, ubnd
Do j = 1, ni
If (evals(j,k)/=0) Then
Write (nout,99995) j, ds(j,k)
Write (nout,99994) j, es(j,k)
If (evals(j,k)==3) Then
Write (nout,99993) j
End If
End If
End Do
End Do
Return
99999 Format (' Segment ',I3,'.')
99998 Format (' Sid = ',I3,', Parent = ',I3,', Level = ',I3,'.')
99997 Format (' Children = (',I3,',',I3,').')
99996 Format (' Bounds (',Es11.4,',',Es11.4,').')
99995 Format (' Integral ',I2,' approximation :',1X,Es11.4,'.')
99994 Format (' Integral ',I2,' error estimate:',1X,Es11.4,'.')
99993 Format (' Integral ',I2, &
' evaluation has been superseded by descendants.')
99992 Format (' Information on subdivision and evaluations over segments.')
99991 Format (' Integral ',I2,' total approximations: ',I3,'.')
End Subroutine display_subdivision_strategy
Subroutine display_option(optstr,optype,ivalue,rvalue,cvalue)
! Subroutine to query optype and print the appropriate option values
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: rvalue
Integer, Intent (In) :: ivalue, optype
Character (*), Intent (In) :: cvalue, optstr
! .. Executable Statements ..
Select Case (optype)
Case (1)
Write (nout,99999) optstr, ivalue
Case (2)
Write (nout,99998) optstr, rvalue
Case (3)
Write (nout,99997) optstr, cvalue
Case (4)
Write (nout,99996) optstr, ivalue, cvalue
Case (5)
Write (nout,99995) optstr, rvalue, cvalue
Case Default
End Select
Flush (nout)
Return
99999 Format (3X,A30,' : ',I16)
99998 Format (3X,A30,' : ',Es16.4)
99997 Format (3X,A30,' : ',12X,A16)
99996 Format (3X,A30,' : ',I16,3X,A16)
99995 Format (3X,A30,' : ',Es16.4,3X,A16)
End Subroutine display_option
End Module d01rafe_mod
Program d01rafe
! .. Use Statements ..
Use d01rafe_mod, Only: display_integration_details, display_option, &
disp_integration_info, nout
Use nag_library, Only: d01raf, d01rcf, d01zkf, d01zlf, nag_wp, x01aaf
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: a, b, pi, rvalue
Integer :: ifail, irevcm, ivalue, j, lcmax, &
lcmin, lcom, ldfm, ldfmrq, lenx, &
lenxrq, licmax, licmin, licom, ni, &
nx, optype, sdfm, sdfmrq, sid
Character (16) :: cvalue
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: com(:), dinest(:), errest(:), &
fm(:,:), opts(:), x(:)
Integer, Allocatable :: icom(:), iopts(:), needi(:)
! .. Intrinsic Procedures ..
Intrinsic :: cos, sin
! .. Executable Statements ..
Continue
Write (nout,*) 'D01RAF Example Program Results'
Write (nout,*)
pi = x01aaf(pi)
! Setup phase.
! Set problem parameters
ni = 2
! Lower (a) and upper (b) bounds
a = 0.0E0_nag_wp
b = pi
Allocate (opts(100),iopts(100))
! Initialize option arrays
ifail = 0
Call d01zkf('Initialize = d01raf',iopts,100,opts,100,ifail)
! Set any non-default options required
Call d01zkf('Quadrature Rule = gk41',iopts,100,opts,100,ifail)
Call d01zkf('Absolute Tolerance = 1.0e-7',iopts,100,opts,100,ifail)
Call d01zkf('Relative Tolerance = 1.0e-7',iopts,100,opts,100,ifail)
! Determine maximum required array lengths
ifail = -1
Call d01rcf(ni,lenxrq,ldfmrq,sdfmrq,licmin,licmax,lcmin,lcmax,iopts, &
opts,ifail)
ldfm = ldfmrq
sdfm = sdfmrq
lenx = lenxrq
licom = licmax
lcom = lcmax
! Allocate remaining arrays
Allocate (icom(licom),needi(ni),com(lcom),fm(ldfm,sdfm),dinest(ni), &
errest(ni),x(lenx))
! Solve phase.
! Use D01RAF to evaluate the definite integrals of:
! f_1 = (x*sin(2*x))*cos(15*x)
! f_2 = (x*sin(2*x))*(x*cos(50*x))
! Set initial irevcm
irevcm = 1
ifail = -1
Do While (irevcm/=0)
Call d01raf(irevcm,ni,a,b,sid,needi,x,lenx,nx,fm,ldfm,dinest,errest, &
iopts,opts,icom,licom,com,lcom,ifail)
Select Case (irevcm)
Case (11)
! Initial returns.
! These will occur during the non-adaptive phase.
! All values must be supplied.
! DINEST and ERREST do not contain approximations
! over the complete interval at this stage.
! Calculate x*sin(2*x), storing the result in fm(2,1:nx) for re-use.
fm(2,1:nx) = x(1:nx)*sin(2.0E0_nag_wp*x(1:nx))
! Calculate f1
fm(1,1:nx) = fm(2,1:nx)*cos(15.0E0_nag_wp*x(1:nx))
! Complete f2 calculation.
fm(2,1:nx) = fm(2,1:nx)*x(1:nx)*cos(50.0E0_nag_wp*x(1:nx))
Case (12)
! Intermediate returns.
! These will occur during the adaptive phase.
! All requested values must be supplied.
! DINEST and ERREST do not contain approximations
! over the complete interval at this stage.
! Calculate x*sin(2*x).
fm(2,1:nx) = x(1:nx)*sin(2.0E0_nag_wp*x(1:nx))
! Calculate f1 if required.
If (needi(1)==1) Then
fm(1,1:nx) = fm(2,1:nx)*cos(15.0E0_nag_wp*x(1:nx))
End If
! Complete f2 calculation if required.
If (needi(2)==1) Then
fm(2,1:nx) = fm(2,1:nx)*x(1:nx)*cos(50.0E0_nag_wp*x(1:nx))
End If
Case (0)
! Final return. Test IFAIL.
Select Case (ifail)
Case (0:3)
! Useful information has been returned.
Case Default
! An unrecoverable error has been detected.
Go To 100
End Select
End Select
End Do
! Query some currently set options and statistics.
ifail = 0
Call d01zlf('Quadrature rule',ivalue,rvalue,cvalue,optype,iopts,opts, &
ifail)
Call display_option('Quadrature rule',optype,ivalue,rvalue,cvalue)
Call d01zlf('Maximum Subdivisions',ivalue,rvalue,cvalue,optype,iopts, &
opts,ifail)
Call display_option('Maximum Subdivisions',optype,ivalue,rvalue,cvalue)
Call d01zlf('Extrapolation',ivalue,rvalue,cvalue,optype,iopts,opts, &
ifail)
Call display_option('Extrapolation',optype,ivalue,rvalue,cvalue)
Call d01zlf('Extrapolation Safeguard',ivalue,rvalue,cvalue,optype,iopts, &
opts,ifail)
Call display_option('Extrapolation safeguard',optype,ivalue,rvalue, &
cvalue)
! Print solution
Write (nout,*)
Write (nout,99999)
Do j = 1, ni
Write (nout,99998) j, needi(j), dinest(j), errest(j)
End Do
! Investigate integration strategy
If (disp_integration_info) Then
Call display_integration_details(ni,iopts,opts,icom,licom,com,lcom)
End If
100 Continue
99999 Format (' Integral | NEEDI | DINEST | ERREST ')
99998 Format (2(1X,I9),2(1X,Es12.4))
End Program d01rafe