!   D01RAF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.
    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