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 3951 : 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 254865 : real(RK) :: h(NSTEP+1),s(NSTEP+1)
131 : procedure(integrand_proc) :: getFunc
132 3921 : ierr = 0_IK
133 3921 : km = nRefinement - 1_IK
134 3921 : h(1) = 1._RK
135 3921 : numFuncEval = 0_IK
136 24501 : do refinementStage = 1, NSTEP
137 24501 : call doQuadTrap(getFunc,lowerLim,upperLim,s(refinementStage),refinementStage,numFuncEvalNew)
138 24501 : numFuncEval = numFuncEval + numFuncEvalNew
139 24501 : if (refinementStage>=nRefinement) then
140 4413 : call doPolInterp(h(refinementStage-km),s(refinementStage-km),nRefinement,0._RK,integral,relativeError,ierr)
141 4413 : if ( abs(relativeError)<=maxRelativeError*abs(integral) .or. ierr/=0_IK ) return
142 : end if
143 20580 : s(refinementStage+1)=s(refinementStage)
144 45081 : 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 36 : 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 1548 : real(RK) :: h(NSTEP+1), s(NSTEP+1)
195 : procedure(integrand_proc) :: getFunc
196 : procedure(integrator_proc) :: integrate
197 36 : ierr = 0_IK
198 36 : km = nRefinement - 1_IK
199 36 : h(1) = 1._RK
200 36 : numFuncEval = 0_IK
201 246 : do refinementStage = 1, NSTEP
202 246 : call integrate(getFunc,lowerLim,upperLim,s(refinementStage),refinementStage,numFuncEvalNew)
203 246 : numFuncEval = numFuncEval + numFuncEvalNew
204 246 : if (refinementStage>=nRefinement) then
205 36 : call doPolInterp(h(refinementStage-km),s(refinementStage-km),nRefinement,0._RK,integral,relativeError,ierr)
206 36 : if ( abs(relativeError) <= maxRelativeError * abs(integral) .or. ierr/=0_IK ) return
207 : end if
208 210 : s(refinementStage+1)=s(refinementStage)
209 456 : 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 4449 : 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 62328 : real(RK) :: den,dif,dift,ho,hp,w,c(ndata),d(ndata)
226 : integer :: i,m,ns
227 4449 : ierr = 0
228 4449 : ns=1
229 4449 : dif=abs(x-xa(1))
230 31164 : do i=1,ndata
231 26715 : dift=abs(x-xa(i))
232 26715 : if (dift<dif) then
233 22266 : ns=i
234 22266 : dif=dift
235 : endif
236 26715 : c(i)=ya(i)
237 31164 : d(i)=ya(i)
238 : end do
239 4449 : y=ya(ns)
240 4449 : ns=ns-1
241 26715 : do m=1,ndata-1
242 91452 : do i=1,ndata-m
243 69186 : ho=xa(i)-x
244 69186 : hp=xa(i+m)-x
245 69186 : w=c(i+1)-d(i)
246 69186 : den=ho-hp
247 69186 : if(den==0._RK) then
248 0 : ierr = 3 ! failure in doPolInterp()
249 0 : return
250 : end if
251 69186 : den=w/den
252 69186 : d(i)=hp*den
253 91452 : c(i)=ho*den
254 : end do
255 22266 : if (2*ns<ndata-m)then
256 0 : dy=c(ns+1)
257 : else
258 22266 : dy=d(ns)
259 22266 : ns=ns-1
260 : endif
261 26715 : y=y+dy
262 : end do
263 : end subroutine doPolInterp
264 :
265 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
266 :
267 24501 : 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 24501 : real(RK) :: del,sum,tnm,x
278 : procedure(integrand_proc) :: getFunc
279 24501 : if (refinementStage==1) then
280 3921 : numFuncEval = 2_IK
281 3921 : integral = 0.5_RK * (upperLim - lowerLim) * ( getFunc(lowerLim) + getFunc(upperLim) )
282 : else
283 20580 : numFuncEval = 2**(refinementStage-2)
284 20580 : tnm = real(numFuncEval,kind=RK)
285 20580 : del = (upperLim-lowerLim) / tnm
286 20580 : x = lowerLim + 0.5_RK * del
287 20580 : sum = 0._RK
288 198003 : do iFuncEval = 1, numFuncEval
289 177423 : sum = sum + getFunc(x)
290 198003 : x = x + del
291 : end do
292 20580 : integral = 0.5_RK * (integral + (upperLim-lowerLim) * sum / tnm)
293 : endif
294 24501 : 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 156 : 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 156 : real(RK) :: ddel,del,summ,x,lowerLimTrans,upperLimTrans
327 156 : real(RK) :: inverseThreeNumFuncEval
328 : integer(IK) :: iFuncEval
329 156 : upperLimTrans = exp(-lowerLim)
330 156 : lowerLimTrans = 0._RK; if (upperLim<LOGHUGE_RK) lowerLimTrans = exp(-upperLim)
331 156 : if (refinementStage==1_IK) then
332 27 : numFuncEval = 1_IK
333 27 : integral = (upperLimTrans-lowerLimTrans)*getTransFunc(0.5_RK*(lowerLimTrans+upperLimTrans))
334 : else
335 129 : numFuncEval = 3**(refinementStage-2)
336 129 : inverseThreeNumFuncEval = ONE_THIRD / numFuncEval
337 129 : del = (upperLimTrans-lowerLimTrans) * inverseThreeNumFuncEval
338 129 : ddel = del + del
339 129 : x = lowerLimTrans + 0.5_RK*del
340 129 : summ = 0._RK
341 31584 : do iFuncEval = 1,numFuncEval
342 31455 : summ = summ + getTransFunc(x)
343 31455 : x = x + ddel
344 31455 : summ = summ + getTransFunc(x)
345 31584 : x = x + del
346 : end do
347 129 : integral = ONE_THIRD * integral + (upperLimTrans-lowerLimTrans) * summ * inverseThreeNumFuncEval
348 129 : numFuncEval = 2_IK * numFuncEval
349 : end if
350 : contains
351 62937 : function getTransFunc(x)
352 : implicit none
353 : real(RK), intent(in) :: x
354 : real(RK) :: getTransFunc
355 62937 : getTransFunc = getFunc(-log(x)) / x
356 62937 : 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 60 : 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 60 : real(RK) :: lowerLimTrans, upperLimTrans, del, ddel, summ, x
391 60 : real(RK) :: inverseThreeNumFuncEval
392 : integer(IK) :: iFuncEval
393 60 : upperLimTrans = 1.0_RK / lowerLim
394 60 : lowerLimTrans = 1.0_RK / upperLim
395 120 : if (refinementStage == 1_IK) then
396 6 : numFuncEval = 1_IK
397 6 : integral = (upperLimTrans-lowerLimTrans) * getTransFunc(0.5_RK * (lowerLimTrans+upperLimTrans))
398 : else
399 54 : numFuncEval = 3**(refinementStage-2)
400 54 : inverseThreeNumFuncEval = ONE_THIRD / numFuncEval
401 54 : del = (upperLimTrans-lowerLimTrans) * inverseThreeNumFuncEval
402 54 : ddel = del + del
403 54 : x = lowerLimTrans + 0.5_RK * del
404 54 : summ = 0._RK
405 59100 : do iFuncEval = 1, numFuncEval
406 59046 : summ = summ + getTransFunc(x)
407 59046 : x = x + ddel
408 59046 : summ = summ + getTransFunc(x)
409 59100 : x = x + del
410 : end do
411 54 : integral = ONE_THIRD * integral + (upperLimTrans-lowerLimTrans) * summ * inverseThreeNumFuncEval
412 54 : numFuncEval = 2_IK * numFuncEval
413 : end if
414 : contains
415 118098 : function getTransFunc(x) result(transFunc)
416 : implicit none
417 : real(RK), intent(in) :: x
418 : real(RK) :: transFunc
419 118098 : transFunc = getFunc(1._RK/x) / x**2
420 118158 : 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 30 : 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 30 : real(RK) :: ddel,del,summ,x
456 30 : real(RK) :: inverseThreeNumFuncEval
457 30 : if (refinementStage==1) then
458 3 : numFuncEval = 1_IK
459 3 : integral = (upperLim-lowerLim) * getFunc( 0.5_RK * (lowerLim+upperLim) )
460 : else
461 27 : numFuncEval = 3_IK**(refinementStage-2)
462 27 : inverseThreeNumFuncEval = ONE_THIRD / numFuncEval
463 27 : del = (upperLim-lowerLim) * inverseThreeNumFuncEval
464 27 : ddel = del+del
465 27 : x = lowerLim + 0.5_RK * del
466 27 : summ = 0._RK
467 29550 : do iFuncEval = 1, numFuncEval
468 29523 : summ = summ + getFunc(x)
469 29523 : x = x + ddel
470 29523 : summ = summ + getFunc(x)
471 29550 : x = x + del
472 : end do
473 27 : integral = ONE_THIRD * integral + (upperLim-lowerLim) * summ * inverseThreeNumFuncEval
474 27 : numFuncEval = 2_IK * numFuncEval
475 : end if
476 30 : end subroutine midpnt
477 :
478 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
479 :
480 : end module Integration_mod
|
