Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3 : !!!! !!!!
4 : !!!! ParaMonte: Parallel Monte Carlo and Machine Learning Library. !!!!
5 : !!!! !!!!
6 : !!!! Copyright (C) 2012-present, The Computational Data Science Lab !!!!
7 : !!!! !!!!
8 : !!!! This file is part of the ParaMonte library. !!!!
9 : !!!! !!!!
10 : !!!! LICENSE !!!!
11 : !!!! !!!!
12 : !!!! https://github.com/cdslaborg/paramonte/blob/main/LICENSE.md !!!!
13 : !!!! !!!!
14 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
16 :
17 : !> \brief
18 : !> This include file contains procedure implementation of the generic interfaces of [pm_quadRomb](@ref pm_quadRomb).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Apr 21, 2017, 1:54 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : !%%%%%%%%%%%%%%%%%%
28 : #if getQuadRomb_ENABLED
29 : !%%%%%%%%%%%%%%%%%%
30 :
31 : #if Clos_ENABLED
32 : real(RKC) , parameter :: FACTOR = 0.5_RKC
33 : real(RKC) , parameter :: SHRINKAGE = 0.25_RKC
34 : integer(IK) , parameter :: POINT_INCREMENT_BASE = 2_IK
35 : #elif Open_ENABLED || PWRL_ENABLED || NEXP_ENABLED || PEXP_ENABLED || LBIS_ENABLED || UBIS_ENABLED
36 : real(RKC) , parameter :: FACTOR = 1._RKC / 3._RKC
37 : real(RKC) , parameter :: SHRINKAGE = 1._RKC / 9._RKC
38 : integer(IK) , parameter :: POINT_INCREMENT_BASE = 3_IK
39 : real(RKC) :: resolutionDouble
40 : #if !Open_ENABLED
41 : real(RKC) :: lbTrans, ubTrans
42 : #endif
43 : #if LBIS_ENABLED || UBIS_ENABLED
44 : real(RKC) :: inverseOneMinusExp, exp2InverseOneMinusExp ! 1 / (1 - exponent), exponent / (1 - exponent)
45 : #endif
46 : #else
47 : #error "Unrecognized interface."
48 : #endif
49 : integer(IK) , parameter :: NSTEP = ceiling(log(epsilon(0._RKC)) / log(SHRINKAGE)) ! 31_IK: Beyond this integer, the newer QuadRombAbscissa values become identical, which automatically fail `setInterp()`.
50 : integer(IK) :: refinementStage, km
51 : real(RKC) :: QuadRombAbscissa(NSTEP + 1), QuadRombProxy(NSTEP + 1)
52 : real(RKC) :: integrationRange, nevalNewInverse, resolution, sumFunc, x
53 : integer(IK) :: nevalNew, ieval
54 : #if EM_ENABLED
55 : real(RKC) :: relerr
56 : #endif
57 : #if NP_ENABLED && Clos_ENABLED
58 3 : neval = 2_IK
59 : #elif NP_ENABLED
60 18 : neval = 1_IK
61 : #endif
62 : ! Set the final normalization of the integral, only necessary when there is singularity at a limit.
63 : #if LBIS_ENABLED || UBIS_ENABLED
64 : #define NORMALIZE_QUAD_ROMB quadRomb = quadRomb * inverseOneMinusExp;
65 : #else
66 : #define NORMALIZE_QUAD_ROMB
67 : #endif
68 : ! Set the transformed limits
69 : #if Clos_ENABLED || Open_ENABLED
70 : #define GET_FUNC(x)getFunc(x)
71 : #define lbTrans lb
72 : #define ubTrans ub
73 : #elif PWRL_ENABLED
74 : #define GET_FUNC(x)getFunc(1._RKC / x) / x**2
75 3 : lbTrans = 1._RKC / ub
76 3 : ubTrans = 1._RKC / lb
77 : #elif NEXP_ENABLED
78 : #define GET_FUNC(x)getFunc(-log(x)) / x
79 6 : lbTrans = exp(-ub)
80 6 : ubTrans = exp(-lb)
81 : #elif PEXP_ENABLED
82 : #define GET_FUNC(x)getFunc(+log(x)) / x
83 6 : lbTrans = exp(lb)
84 6 : ubTrans = exp(ub)
85 : #elif LBIS_ENABLED || UBIS_ENABLED
86 : lbTrans = 0._RKC
87 16 : ubTrans = (ub - lb)**(1._RKC + real(Interval%exponent, RKC))
88 16 : inverseOneMinusExp = 1._RKC / (1._RKC + real(Interval%exponent, RKC))
89 16 : exp2InverseOneMinusExp = -real(Interval%exponent, RKC) * inverseOneMinusExp
90 16 : CHECK_ASSERTION(__LINE__, -1._RKC < Interval%exponent .and. Interval%exponent <= 0._RKC, \
91 : SK_"The conditions `-1. < Interval%exponent .and. Interval%exponent <= 0.` must hold: Interval%exponent = "//getStr(Interval%exponent))
92 : #if LBIS_ENABLED
93 : #define GET_FUNC(x)x**exp2InverseOneMinusExp * getFunc(lb + x**inverseOneMinusExp)
94 : #elif UBIS_ENABLED
95 : #define GET_FUNC(x)x**exp2InverseOneMinusExp * getFunc(ub - x**inverseOneMinusExp)
96 : #endif
97 : #endif
98 21 : integrationRange = ubTrans - lbTrans
99 37 : CHECK_ASSERTION(__LINE__, integrationRange >= 0._RKC, SK_"The input lower and upper integration bounds [lb, ub] must satisfy `lb < ub`.")
100 37 : CHECK_ASSERTION(__LINE__, nref > 0_IK, SK_"The input refinement count in the Romberg method must satisfy `nref > 0`. nref = "//getStr(nref))
101 111 : CHECK_ASSERTION(__LINE__, nref <= NSTEP, SK_"The input refinement count in the Romberg method must satisfy `nref > 0`. nref, NSTEP = "//getStr([nref, NSTEP]))
102 37 : km = nref - 1_IK
103 : refinementStage = 1_IK
104 37 : QuadRombAbscissa(refinementStage) = 1._RKC
105 : #if Clos_ENABLED
106 3 : QuadRombProxy(refinementStage) = 0.5_RKC * integrationRange * (GET_FUNC(lbTrans) + GET_FUNC(ubTrans))
107 : #else
108 34 : QuadRombProxy(refinementStage) = integrationRange * GET_FUNC((0.5_RKC * (lbTrans + ubTrans))) ! \warning extra parentheses are important
109 : #endif
110 : ! Define the evaluation instruction.
111 : #define EVAL_QUAD_ROMB \
112 : if (refinementStage >= nref) then; \
113 : call setExtrap(monopol, QuadRombAbscissa(refinementStage - km : refinementStage), QuadRombProxy(refinementStage - km : refinementStage), 0._RKC, quadRomb, relerr); \
114 : relerr = abs(relerr); \
115 : if (relerr <= tol * abs(quadRomb)) then; \
116 : NORMALIZE_QUAD_ROMB \
117 : return; \
118 : end if; \
119 : end if; \
120 : QuadRombProxy(refinementStage + 1_IK) = QuadRombProxy(refinementStage); \
121 : QuadRombAbscissa(refinementStage + 1_IK) = SHRINKAGE * QuadRombAbscissa(refinementStage);
122 : ! Compute the integral.
123 37 : EVAL_QUAD_ROMB
124 200 : do refinementStage = 2_IK, NSTEP
125 3 : nevalNew = POINT_INCREMENT_BASE**(refinementStage - 2_IK)
126 200 : nevalNewInverse = 1._RKC / real(nevalNew, RKC)
127 : #if Clos_ENABLED
128 9 : resolution = integrationRange * nevalNewInverse
129 : #else
130 191 : resolution = integrationRange * nevalNewInverse * FACTOR
131 191 : resolutionDouble = resolution + resolution
132 : #endif
133 200 : x = lbTrans + 0.5_RKC * resolution
134 : sumFunc = 0._RKC
135 120462 : do ieval = 1_IK, nevalNew
136 120262 : sumFunc = sumFunc + GET_FUNC(x)
137 : #if !Clos_ENABLED
138 120241 : x = x + resolutionDouble
139 120241 : sumFunc = sumFunc + GET_FUNC(x)
140 : #endif
141 120462 : x = x + resolution
142 : end do
143 200 : QuadRombProxy(refinementStage) = FACTOR * (QuadRombProxy(refinementStage) + integrationRange * sumFunc * nevalNewInverse)
144 : #if NP_ENABLED && Clos_ENABLED
145 9 : neval = neval + nevalNew
146 : #elif NP_ENABLED
147 87 : neval = neval + nevalNew * 2
148 : #endif
149 200 : EVAL_QUAD_ROMB
150 : end do
151 : #if EM_ENABLED
152 0 : error stop __FILE__//"@getQuadRomb(): The Romberg integration failed to converge."
153 : #elif EP_ENABLED
154 0 : relerr = -huge(relerr)
155 : #else
156 : #error "Unrecognized interface."
157 : #endif
158 :
159 : #else
160 : !%%%%%%%%%%%%%%%%%%%%%%%%
161 : #error "Unrecognized interface."
162 : !%%%%%%%%%%%%%%%%%%%%%%%%
163 : #endif
164 : #undef NORMALIZE_QUAD_ROMB
165 : #undef GET_FUNC
166 : #undef lbTrans
167 : #undef ubTrans
|
