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
|
