The ParaMonte Documentation Website
Current view: top level - kernel - Integration_mod.f90 (source / functions) Hit Total Coverage
Test: ParaMonte 1.5.1 :: Serial Kernel - Code Coverage Report Lines: 139 144 96.5 %
Date: 2021-01-08 13:03:42 Functions: 9 9 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       2             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       3             : !!!!
       4             : !!!!   MIT License
       5             : !!!!
       6             : !!!!   ParaMonte: plain powerful parallel Monte Carlo library.
       7             : !!!!
       8             : !!!!   Copyright (C) 2012-present, The Computational Data Science Lab
       9             : !!!!
      10             : !!!!   This file is part of the ParaMonte library.
      11             : !!!!
      12             : !!!!   Permission is hereby granted, free of charge, to any person obtaining a
      13             : !!!!   copy of this software and associated documentation files (the "Software"),
      14             : !!!!   to deal in the Software without restriction, including without limitation
      15             : !!!!   the rights to use, copy, modify, merge, publish, distribute, sublicense,
      16             : !!!!   and/or sell copies of the Software, and to permit persons to whom the
      17             : !!!!   Software is furnished to do so, subject to the following conditions:
      18             : !!!!
      19             : !!!!   The above copyright notice and this permission notice shall be
      20             : !!!!   included in all copies or substantial portions of the Software.
      21             : !!!!
      22             : !!!!   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : !!!!   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : !!!!   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : !!!!   IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
      26             : !!!!   DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
      27             : !!!!   OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE
      28             : !!!!   OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : !!!!
      30             : !!!!   ACKNOWLEDGMENT
      31             : !!!!
      32             : !!!!   ParaMonte is an honor-ware and its currency is acknowledgment and citations.
      33             : !!!!   As per the ParaMonte library license agreement terms, if you use any parts of
      34             : !!!!   this library for any purposes, kindly acknowledge the use of ParaMonte in your
      35             : !!!!   work (education/research/industry/development/...) by citing the ParaMonte
      36             : !!!!   library as described on this page:
      37             : !!!!
      38             : !!!!       https://github.com/cdslaborg/paramonte/blob/main/ACKNOWLEDGMENT.md
      39             : !!!!
      40             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      41             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      42             : 
      43             : !>  \brief This module contains classes and procedures to perform numerical integrations.
      44             : !>  \author Amir Shahmoradi
      45             : 
      46             : module Integration_mod
      47             : 
      48             :     use Constants_mod, only: IK, RK
      49             : 
      50             :     implicit none
      51             : 
      52             :     character(*), parameter :: MODULE_NAME = "@Integration_mod"
      53             : 
      54             :     real(RK), parameter :: ONE_THIRD = 1._RK / 3._RK
      55             : 
      56             :     character(len=117)  :: ErrorMessage(3) = &
      57             :     [ "@Integration_mod@doQuadRombClosed(): Too many steps in doQuadRombClosed().                                           " &
      58             :     , "@Integration_mod@doQuadRombClosed()@doPolInterp(): Input lowerLim, upperLim to doQuadRombClosed() are likely equal.  " &
      59             :     , "@Integration_mod@doQuadRombOpen()@doPolInterp(): Input lowerLim, upperLim to doQuadRombOpen() are likely equal.      " &
      60             :     ]
      61             : 
      62             :     abstract interface
      63             : 
      64             :         !> The abstract interface of the integrand.
      65             :         function integrand_proc(x) result(integrand)
      66             :             use Constants_mod, only: RK
      67             :             implicit none
      68             :             real(RK), intent(in)  :: x
      69             :             real(RK)              :: integrand
      70             :         end function integrand_proc
      71             : 
      72             :         !> The abstract interface of the integrator.
      73             :         subroutine integrator_proc(getFunc,lowerLim,upperLim,integral,refinementStage,numFuncEval)
      74             :             use Constants_mod, only: IK, RK
      75             :             import :: integrand_proc
      76             :             implicit none
      77             :             integer(IK) , intent(in)        :: refinementStage
      78             :             real(RK)    , intent(in)        :: lowerLim,upperLim
      79             :             integer(IK) , intent(out)       :: numFuncEval
      80             :             real(RK)    , intent(inout)     :: integral
      81             :             procedure(integrand_proc)       :: getFunc
      82             :         end subroutine integrator_proc
      83             : 
      84             :     end interface
      85             : 
      86             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      87             : 
      88             : contains
      89             : 
      90             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      91             : 
      92             :     !> \brief
      93             :     !> Return the integral of function getFunc in the closed range [lowerLim,upperLim] using Adaptive Romberg extrapolation method.
      94             :     !> \param[in]   getFunc             :   The input function to be integrated. It must have the interface specified by
      95             :     !!                                      [integrand_proc](@ref integrand_proc).
      96             :     !> \param[in]   lowerLim            :   The lower limit of integration.
      97             :     !> \param[in]   upperLim            :   The upper limit of integration.
      98             :     !> \param[in]   maxRelativeError    :   The tolerance that once reach, the integration will stop.
      99             :     !> \param[in]   nRefinement         :   The number of refinements in the Romberg method. A number between 4-6 is typically good.
     100             :     !> \param[out]  integral            :   The result of integration.
     101             :     !> \param[out]  relativeError       :   The final estimated relative error in the result. By definition, this is always
     102             :     !!                                      smaller than `maxRelativeError`.
     103             :     !> \param[out]  numFuncEval         :   The number of function evaluations made.
     104             :     !> \param[out]  ierr                :   Integer flag indicating whether an error has occurred.
     105             :     !!                                      If `ierr == 0`, then no error has occurred.
     106             :     !!                                      Otherwise, the integer value points to the corresponding element of
     107             :     !!                                      [ErrorMessage](@ref errormessage) array.
     108        1317 :     recursive subroutine doQuadRombClosed   ( getFunc &
     109             :                                             , lowerLim &
     110             :                                             , upperLim &
     111             :                                             , maxRelativeError &
     112             :                                             , nRefinement &
     113             :                                             , integral &
     114             :                                             , relativeError &
     115             :                                             , numFuncEval &
     116             :                                             , ierr &
     117             :                                             )
     118             :         !checked by Joshua Osborne on 5/28/2020 at 9:06pm
     119             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     120             :         !DEC$ ATTRIBUTES DLLEXPORT :: doQuadRombClosed
     121             : #endif
     122             :         use, intrinsic :: iso_fortran_env, only: DPI => int64
     123             :         implicit none
     124             :         integer(IK) , intent(in)    :: nRefinement
     125             :         real(RK)    , intent(in)    :: lowerLim,upperLim,maxRelativeError
     126             :         real(RK)    , intent(out)   :: integral, relativeError
     127             :         integer(IK) , intent(out)   :: ierr, numFuncEval
     128             :         integer(IK)                 :: refinementStage,km,numFuncEvalNew
     129             :         integer(IK), parameter      :: NSTEP = 31_IK
     130       84955 :         real(RK)                    :: h(NSTEP+1),s(NSTEP+1)
     131             :         procedure(integrand_proc)   :: getFunc
     132        1307 :         ierr = 0_IK
     133        1307 :         km = nRefinement - 1_IK
     134        1307 :         h(1) = 1._RK
     135        1307 :         numFuncEval = 0_IK
     136        8167 :         do refinementStage = 1, NSTEP
     137        8167 :             call doQuadTrap(getFunc,lowerLim,upperLim,s(refinementStage),refinementStage,numFuncEvalNew)
     138        8167 :             numFuncEval = numFuncEval + numFuncEvalNew
     139        8167 :             if (refinementStage>=nRefinement) then
     140        1471 :                 call doPolInterp(h(refinementStage-km),s(refinementStage-km),nRefinement,0._RK,integral,relativeError,ierr)
     141        1471 :                 if ( abs(relativeError)<=maxRelativeError*abs(integral) .or. ierr/=0_IK ) return
     142             :             end if
     143        6860 :             s(refinementStage+1)=s(refinementStage)
     144       15027 :             h(refinementStage+1)=0.25_RK*h(refinementStage)
     145             :         end do
     146           0 :         ierr = 1_IK    ! too many steps in doQuadRombClosed()
     147             :   end subroutine doQuadRombClosed
     148             : 
     149             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     150             : 
     151             :     !> \brief
     152             :     !> Return the integral of function getFunc in the open range `[lowerLim,upperLim]` using Adaptive Romberg extrapolation method.
     153             :     !> \param[in]   getFunc             :   The input function to be integrated. It must have the interface specified by
     154             :     !!                                      [integrand_proc](@ref integrand_proc).
     155             :     !> \param[in]   integrate           :   The integrator routine. It can be either [midexp](@ref midexp), [midpnt](@ref midpnt),
     156             :     !!                                      or [midinf](@ref midinf).
     157             :     !> \param[in]   lowerLim            :   The lower limit of integration.
     158             :     !> \param[in]   upperLim            :   The upper limit of integration.
     159             :     !> \param[in]   maxRelativeError    :   The tolerance that once reach, the integration will stop.
     160             :     !> \param[in]   nRefinement         :   The number of refinements in the Romberg method. A number between 4-6 is typically good.
     161             :     !> \param[out]  integral            :   The result of integration.
     162             :     !> \param[out]  relativeError       :   The final estimated relative error in the result. By definition, this is always
     163             :     !!                                      smaller than `maxRelativeError`.
     164             :     !> \param[out]  numFuncEval         :   The number of function evaluations made.
     165             :     !> \param[out]  ierr                :   Integer flag indicating whether an error has occurred.
     166             :     !!                                      If `ierr == 0`, then no error has occurred.
     167             :     !!                                      Otherwise, the integer value points to the corresponding element
     168             :     !!                                      of [ErrorMessage](@ref errormessage) array.
     169             :     !> \remark
     170             :     !> Checked by Joshua Osborne on 5/28/2020 at 9:03 pm.
     171          12 :     recursive subroutine doQuadRombOpen ( getFunc &
     172             :                                         , integrate &
     173             :                                         , lowerLim &
     174             :                                         , upperLim &
     175             :                                         , maxRelativeError &
     176             :                                         , nRefinement &
     177             :                                         , integral &
     178             :                                         , relativeError &
     179             :                                         , numFuncEval &
     180             :                                         , ierr &
     181             :                                         )
     182             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     183             :         !DEC$ ATTRIBUTES DLLEXPORT :: doQuadRombOpen
     184             : #endif
     185             :         use, intrinsic :: iso_fortran_env, only: DPI => int64
     186             :         implicit none
     187             :         integer(IK) , intent(in)    :: nRefinement
     188             :         real(RK)    , intent(in)    :: lowerLim,upperLim,maxRelativeError
     189             :         real(RK)    , intent(out)   :: integral, relativeError
     190             :         integer(IK) , intent(out)   :: numFuncEval, ierr
     191             :         integer(IK)                 :: refinementStage,km,numFuncEvalNew
     192             :         integer(IK) , parameter     :: NSTEP = 20_IK
     193             :         real(RK)    , parameter     :: ONE_OVER_NINE = 1._RK / 9._RK
     194         516 :         real(RK)                    :: h(NSTEP+1), s(NSTEP+1)
     195             :         procedure(integrand_proc)   :: getFunc
     196             :         procedure(integrator_proc)  :: integrate
     197          12 :         ierr = 0_IK
     198          12 :         km = nRefinement - 1_IK
     199          12 :         h(1) = 1._RK
     200          12 :         numFuncEval = 0_IK
     201          82 :         do refinementStage = 1, NSTEP
     202          82 :             call integrate(getFunc,lowerLim,upperLim,s(refinementStage),refinementStage,numFuncEvalNew)
     203          82 :             numFuncEval = numFuncEval + numFuncEvalNew
     204          82 :             if (refinementStage>=nRefinement) then
     205          12 :                 call doPolInterp(h(refinementStage-km),s(refinementStage-km),nRefinement,0._RK,integral,relativeError,ierr)
     206          12 :                 if ( abs(relativeError) <= maxRelativeError * abs(integral) .or. ierr/=0_IK ) return
     207             :             end if
     208          70 :             s(refinementStage+1)=s(refinementStage)
     209         152 :             h(refinementStage+1)=h(refinementStage)*ONE_OVER_NINE
     210             :         end do
     211           0 :         ierr = 2_IK    ! too many steps in doQuadRombOpen()
     212             :   end subroutine doQuadRombOpen
     213             : 
     214             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     215             : 
     216        1483 :     recursive subroutine doPolInterp(xa,ya,ndata,x,y,dy,ierr)
     217             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     218             :         !DEC$ ATTRIBUTES DLLEXPORT :: doPolInterp
     219             : #endif
     220             :         implicit none
     221             :         integer , intent(in)    :: ndata
     222             :         real(RK), intent(in)    :: x,xa(ndata),ya(ndata)
     223             :         real(RK), intent(out)   :: dy,y
     224             :         integer , intent(out)   :: ierr
     225       20776 :         real(RK)                :: den,dif,dift,ho,hp,w,c(ndata),d(ndata)
     226             :         integer                 :: i,m,ns
     227        1483 :         ierr = 0
     228        1483 :         ns=1
     229        1483 :         dif=abs(x-xa(1))
     230       10388 :         do i=1,ndata
     231        8905 :             dift=abs(x-xa(i))
     232        8905 :             if (dift<dif) then
     233        7422 :                 ns=i
     234        7422 :                 dif=dift
     235             :             endif
     236        8905 :             c(i)=ya(i)
     237       10388 :             d(i)=ya(i)
     238             :         end do
     239        1483 :         y=ya(ns)
     240        1483 :         ns=ns-1
     241        8905 :         do m=1,ndata-1
     242       30484 :             do i=1,ndata-m
     243       23062 :                 ho=xa(i)-x
     244       23062 :                 hp=xa(i+m)-x
     245       23062 :                 w=c(i+1)-d(i)
     246       23062 :                 den=ho-hp
     247       23062 :                 if(den==0._RK) then
     248           0 :                     ierr = 3 ! failure in doPolInterp()
     249           0 :                     return
     250             :                 end if
     251       23062 :                 den=w/den
     252       23062 :                 d(i)=hp*den
     253       30484 :                 c(i)=ho*den
     254             :             end do
     255        7422 :             if (2*ns<ndata-m)then
     256           0 :                 dy=c(ns+1)
     257             :             else
     258        7422 :                 dy=d(ns)
     259        7422 :                 ns=ns-1
     260             :             endif
     261        8905 :             y=y+dy
     262             :         end do
     263             :     end subroutine doPolInterp
     264             : 
     265             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     266             : 
     267        8167 :     recursive subroutine doQuadTrap(getFunc,lowerLim,upperLim,integral,refinementStage,numFuncEval)
     268             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     269             :         !DEC$ ATTRIBUTES DLLEXPORT :: doQuadTrap
     270             : #endif
     271             :         implicit none
     272             :         integer(IK), intent(in)     :: refinementStage
     273             :         real(RK), intent(in)        :: lowerLim,upperLim
     274             :         real(RK), intent(inout)     :: integral
     275             :         integer(IK), intent(out)    :: numFuncEval
     276             :         integer(IK)                 :: iFuncEval
     277        8167 :         real(RK)                    :: del,sum,tnm,x
     278             :         procedure(integrand_proc)   :: getFunc
     279        8167 :         if (refinementStage==1) then
     280        1307 :             numFuncEval = 2_IK
     281        1307 :             integral = 0.5_RK * (upperLim - lowerLim) * ( getFunc(lowerLim) + getFunc(upperLim) )
     282             :         else
     283        6860 :             numFuncEval = 2**(refinementStage-2)
     284        6860 :             tnm = real(numFuncEval,kind=RK)
     285        6860 :             del = (upperLim-lowerLim) / tnm
     286        6860 :             x = lowerLim + 0.5_RK * del
     287        6860 :             sum = 0._RK
     288       66001 :             do iFuncEval = 1, numFuncEval
     289       59141 :                 sum = sum + getFunc(x)
     290       66001 :                 x = x + del
     291             :             end do
     292        6860 :             integral = 0.5_RK * (integral + (upperLim-lowerLim) * sum / tnm)
     293             :         endif
     294        8167 :     end subroutine doQuadTrap
     295             : 
     296             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     297             : 
     298             :     !> \brief
     299             :     !> Return the refinement of the integration of an exponentially-decaying function on a semi-infinite.
     300             :     !> This function is an integrator driver to be passed to [doQuadRombOpen](@ref doquadrombopen).
     301             :     !>
     302             :     !> \param[in]       getFunc         :   The input function to be integrated. It must have the interface specified by
     303             :     !!                                      [integrand_proc](@ref integrand_proc).
     304             :     !> \param[in]       lowerLim        :   The lower limit of integration.
     305             :     !> \param[in]       upperLim        :   The upper limit of integration (typically set to `huge(1._RK)` to represent \f$+\infty\f$).
     306             :     !> \param[inout]    integral        :   The result of integration.
     307             :     !> \param[in]       refinementStage :   The number of refinements since the first call to the integrator.
     308             :     !> \param[out]      numFuncEval     :   The number of function evaluations made.
     309             :     !>
     310             :     !> \remark
     311             :     !> It is expected that `upperLim > lowerLim > 0.0` must hold for this integrator to function properly.
     312             :     !>
     313             :     !> \remark
     314             :     !> Tested by Joshua Osborne on 5/28/2020 at 8:58 pm.
     315          52 :     recursive subroutine midexp(getFunc,lowerLim,upperLim,integral,refinementStage,numFuncEval)
     316             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     317             :         !DEC$ ATTRIBUTES DLLEXPORT :: midexp
     318             : #endif
     319             :         use Constants_mod, only: LOGHUGE_RK
     320             :         implicit none
     321             :         integer(IK) , intent(in)    :: refinementStage
     322             :         real(RK)    , intent(in)    :: lowerLim,upperLim
     323             :         integer(IK) , intent(out)   :: numFuncEval
     324             :         real(RK)    , intent(inout) :: integral
     325             :         procedure(integrand_proc)   :: getFunc
     326          52 :         real(RK)                    :: ddel,del,summ,x,lowerLimTrans,upperLimTrans
     327          52 :         real(RK)                    :: inverseThreeNumFuncEval
     328             :         integer(IK)                 :: iFuncEval
     329          52 :         upperLimTrans = exp(-lowerLim)
     330          52 :         lowerLimTrans = 0._RK; if (upperLim<LOGHUGE_RK) lowerLimTrans = exp(-upperLim)
     331          52 :         if (refinementStage==1_IK) then
     332           9 :             numFuncEval = 1_IK
     333           9 :             integral = (upperLimTrans-lowerLimTrans)*getTransFunc(0.5_RK*(lowerLimTrans+upperLimTrans))
     334             :         else
     335          43 :             numFuncEval = 3**(refinementStage-2)
     336          43 :             inverseThreeNumFuncEval = ONE_THIRD / numFuncEval
     337          43 :             del = (upperLimTrans-lowerLimTrans) * inverseThreeNumFuncEval
     338          43 :             ddel = del + del
     339          43 :             x = lowerLimTrans + 0.5_RK*del
     340          43 :             summ = 0._RK
     341       10528 :             do iFuncEval = 1,numFuncEval
     342       10485 :                 summ = summ + getTransFunc(x)
     343       10485 :                 x = x + ddel
     344       10485 :                 summ = summ + getTransFunc(x)
     345       10528 :                 x = x + del
     346             :             end do
     347          43 :             integral = ONE_THIRD * integral + (upperLimTrans-lowerLimTrans) * summ * inverseThreeNumFuncEval
     348          43 :             numFuncEval = 2_IK * numFuncEval
     349             :         end if
     350             :     contains
     351       20979 :         function getTransFunc(x)
     352             :             implicit none
     353             :             real(RK), intent(in)    :: x
     354             :             real(RK)                :: getTransFunc
     355       20979 :             getTransFunc = getFunc(-log(x)) / x
     356       20979 :         end function getTransFunc
     357             :     end subroutine midexp
     358             : 
     359             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     360             : 
     361             :     !> \brief
     362             :     !> This routine is an exact replacement for [midpnt](@ref midpnt), i.e., returns as `integral` the nth stage of refinement
     363             :     !> of the integral of funk from `lowerLim` to `upperLim`, except that the function is evaluated at evenly spaced
     364             :     !> points in `1/x` rather than in `x`. This allows the upper limit `upperLim` to be as large and positive
     365             :     !> as the computer allows, or the lower limit `lowerLim` to be as large and negative, but not both.
     366             :     !>
     367             :     !> \param[in]       getFunc         :   The input function to be integrated. It must have the interface specified by
     368             :     !!                                      [integrand_proc](@ref integrand_proc).
     369             :     !> \param[in]       lowerLim        :   The lower limit of integration.
     370             :     !> \param[in]       upperLim        :   The upper limit of integration.
     371             :     !> \param[inout]    integral        :   The result of integration.
     372             :     !> \param[in]       refinementStage :   The number of refinements since the first call to the integrator.
     373             :     !> \param[out]      numFuncEval     :   The number of function evaluations made.
     374             :     !>
     375             :     !> \warning
     376             :     !> `lowerLim * upperLim > 0.0` must hold for this integrator to function properly.
     377             :     !>
     378             :     !> \remark
     379             :     !> Tested by Joshua Osborne on 5/28/2020 at 8:58 pm.
     380          20 :     subroutine midinf(getFunc,lowerLim,upperLim,integral,refinementStage,numFuncEval)
     381             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     382             :         !DEC$ ATTRIBUTES DLLEXPORT :: midinf
     383             : #endif
     384             :         implicit none
     385             :         real(RK)    , intent(in)    :: lowerLim,upperLim
     386             :         integer(IK) , intent(in)    :: refinementStage
     387             :         integer(IK) , intent(out)   :: numFuncEval
     388             :         real(RK)    , intent(inout) :: integral
     389             :         procedure(integrand_proc)   :: getFunc
     390          20 :         real(RK)                    :: lowerLimTrans, upperLimTrans, del, ddel, summ, x
     391          20 :         real(RK)                    :: inverseThreeNumFuncEval
     392             :         integer(IK)                 :: iFuncEval
     393          20 :         upperLimTrans = 1.0_RK / lowerLim
     394          20 :         lowerLimTrans = 1.0_RK / upperLim
     395          40 :         if (refinementStage == 1_IK) then
     396           2 :             numFuncEval = 1_IK
     397           2 :             integral = (upperLimTrans-lowerLimTrans) * getTransFunc(0.5_RK * (lowerLimTrans+upperLimTrans))
     398             :         else
     399          18 :             numFuncEval = 3**(refinementStage-2)
     400          18 :             inverseThreeNumFuncEval = ONE_THIRD / numFuncEval
     401          18 :             del = (upperLimTrans-lowerLimTrans) * inverseThreeNumFuncEval
     402          18 :             ddel = del + del
     403          18 :             x = lowerLimTrans + 0.5_RK * del
     404          18 :             summ = 0._RK
     405       19700 :             do iFuncEval = 1, numFuncEval
     406       19682 :                 summ = summ + getTransFunc(x)
     407       19682 :                 x = x + ddel
     408       19682 :                 summ = summ + getTransFunc(x)
     409       19700 :                 x = x + del
     410             :             end do
     411          18 :             integral = ONE_THIRD * integral + (upperLimTrans-lowerLimTrans) * summ * inverseThreeNumFuncEval
     412          18 :             numFuncEval = 2_IK * numFuncEval
     413             :         end if
     414             :     contains
     415       39366 :         function getTransFunc(x) result(transFunc)
     416             :             implicit none
     417             :             real(RK), intent(in)    :: x
     418             :             real(RK)                :: transFunc
     419       39366 :             transFunc = getFunc(1._RK/x) / x**2
     420       39386 :         end function getTransFunc
     421             :     end subroutine midinf
     422             : 
     423             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     424             : 
     425             :     !> \brief
     426             :     !> This routine computes the nth stage of refinement of an extended midpoint rule.
     427             :     !> When called with `n = 1`, the routine returns as `integral` the crudest estimate of \f$\int_a^b f(x) ~ dx\f$.
     428             :     !> Subsequent calls with `n = 2, 3, ...` (in that sequential order) will improve the accuracy of `integral` by adding
     429             :     !> `(2/3) * 3n-1` additional interior points.
     430             :     !>
     431             :     !> \param[in]       getFunc         :   The input function to be integrated. It must have the interface specified by
     432             :     !!                                      [integrand_proc](@ref integrand_proc).
     433             :     !> \param[in]       lowerLim        :   The lower limit of integration.
     434             :     !> \param[in]       upperLim        :   The upper limit of integration.
     435             :     !> \param[inout]    integral        :   The result of integration.
     436             :     !> \param[in]       refinementStage :   The number of refinements since the first call to the integrator.
     437             :     !> \param[out]      numFuncEval     :   The number of function evaluations made.
     438             :     !>
     439             :     !> \warning
     440             :     !> The argument `integral` should not be modified between sequential calls.
     441             :     !>
     442             :     !> \remark
     443             :     !> Tested by Joshua Osborne on 5/28/2020 at 8:55 pm.
     444          10 :     subroutine midpnt(getFunc,lowerLim,upperLim,integral,refinementStage,numFuncEval)
     445             : #if INTEL_COMPILER_ENABLED && defined DLL_ENABLED && (OS_IS_WINDOWS || defined OS_IS_DARWIN)
     446             :         !DEC$ ATTRIBUTES DLLEXPORT :: midpnt
     447             : #endif
     448             :         implicit none
     449             :         integer(IK) , intent(in)    :: refinementStage
     450             :         real(RK)    , intent(in)    :: lowerLim, upperLim
     451             :         real(RK)    , intent(inout) :: integral
     452             :         integer(IK) , intent(out)   :: numFuncEval
     453             :         procedure(integrand_proc)   :: getFunc
     454             :         integer(IK)                 :: iFuncEval
     455          10 :         real(RK)                    :: ddel,del,summ,x
     456          10 :         real(RK)                    :: inverseThreeNumFuncEval
     457          10 :         if (refinementStage==1) then
     458           1 :             numFuncEval = 1_IK
     459           1 :             integral = (upperLim-lowerLim) * getFunc( 0.5_RK * (lowerLim+upperLim) )
     460             :         else
     461           9 :             numFuncEval = 3_IK**(refinementStage-2)
     462           9 :             inverseThreeNumFuncEval = ONE_THIRD / numFuncEval
     463           9 :             del = (upperLim-lowerLim) * inverseThreeNumFuncEval
     464           9 :             ddel = del+del
     465           9 :             x = lowerLim + 0.5_RK * del
     466           9 :             summ = 0._RK
     467        9850 :             do iFuncEval = 1, numFuncEval
     468        9841 :                 summ = summ + getFunc(x)
     469        9841 :                 x = x + ddel
     470        9841 :                 summ = summ + getFunc(x)
     471        9850 :                 x = x + del
     472             :             end do
     473           9 :             integral = ONE_THIRD * integral + (upperLim-lowerLim) * summ * inverseThreeNumFuncEval
     474           9 :             numFuncEval = 2_IK * numFuncEval
     475             :         end if
     476          10 :     end subroutine midpnt
     477             : 
     478             : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     479             : 
     480             : end module Integration_mod

ParaMonte: Plain Powerful Parallel Monte Carlo Library 
The Computational Data Science Lab
© Copyright 2012 - 2021