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_quadTest](@ref pm_quadTest).
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 test_isFailedQuad_ENABLED
29 : !%%%%%%%%%%%%%%%%%%%%%%%%
30 :
31 : integer(IK) :: numFuncEval
32 : real(RKC) :: integral, abserr, lb, ub, truth
33 : real(RKC) , allocatable :: break(:)
34 :
35 24 : truth = real(integrand%integral, RKC)
36 24 : lb = real(integrand%lb, RKC)
37 24 : ub = real(integrand%ub, RKC)
38 :
39 24 : call disp%show("integrand%desc")
40 24 : call disp%show( integrand%desc , deliml = SK_"""" )
41 :
42 24 : if (present(abstol)) then
43 0 : call disp%show("abstol")
44 0 : call disp%show( abstol )
45 : end if
46 :
47 24 : if (present(reltol)) then
48 0 : call disp%show("reltol")
49 0 : call disp%show( reltol )
50 : end if
51 :
52 24 : call disp%skip()
53 24 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
54 24 : call disp%show("! "//integrand%desc)
55 24 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
56 24 : call disp%skip()
57 :
58 24 : call disp%skip()
59 24 : call disp%show("! Regular adaptive global quadrature.")
60 24 : call disp%skip()
61 :
62 24 : call disp%show("if (isFailedQuad(getFunc, lb, ub, integral, abserr, neval = numFuncEval)) call disp%show(' ******** Integration did NOT converge. ********', tmsize = 1_IK, bmsize = 1_IK)")
63 24 : if (isFailedQuad(getFunc, lb, ub, integral, abserr, neval = numFuncEval)) call disp%show(' ******** Integration did NOT converge. ********', tmsize = 1_IK, bmsize = 1_IK)
64 24 : call disp%show("getStr([truth, integral, abserr, abs(integral - truth)])//SK_' (unbiased)? '//getStr(abs(integral - truth) <= abserr)")
65 120 : call disp%show( getStr([truth, integral, abserr, abs(integral - truth)])//SK_' (unbiased)? '//getStr(abs(integral - truth) <= abserr) )
66 24 : call disp%show("numFuncEval")
67 24 : call disp%show( numFuncEval )
68 24 : call disp%skip()
69 :
70 24 : call disp%skip()
71 24 : call disp%show("! Assisted adaptive global quadrature by the Wynn Epsilon extrapolation.")
72 24 : call disp%skip()
73 :
74 24 : call disp%show("if (isFailedQuad(getFunc, lb, ub, weps, integral, abserr, neval = numFuncEval)) call disp%show(' ******** Integration did NOT converge. ********', tmsize = 1_IK, bmsize = 1_IK)")
75 24 : if (isFailedQuad(getFunc, lb, ub, weps, integral, abserr, neval = numFuncEval)) call disp%show(' ******** Integration did NOT converge. ********', tmsize = 1_IK, bmsize = 1_IK)
76 24 : call disp%show("getStr([truth, integral, abserr, abs(integral - truth)])//SK_' (unbiased)? '//getStr(abs(integral - truth) <= abserr)")
77 120 : call disp%show( getStr([truth, integral, abserr, abs(integral - truth)])//SK_' (unbiased)? '//getStr(abs(integral - truth) <= abserr) )
78 24 : call disp%show("numFuncEval")
79 24 : call disp%show( numFuncEval )
80 24 : call disp%skip()
81 :
82 24 : if (allocated(integrand%break)) then
83 :
84 6 : call disp%skip()
85 6 : call disp%show("! Assisted adaptive global quadrature by specifying points of difficulties of the integrand.")
86 6 : call disp%skip()
87 :
88 6 : call disp%show("break = integrand%break")
89 36 : break = integrand%break
90 6 : call disp%show("break")
91 6 : call disp%show( break )
92 6 : call disp%show("if (isFailedQuad(getFunc, lb, ub, break, integral, abserr, neval = numFuncEval)) call disp%show(' ******** Integration did NOT converge. ********', tmsize = 1_IK, bmsize = 1_IK)")
93 6 : if (isFailedQuad(getFunc, lb, ub, break, integral, abserr, neval = numFuncEval)) call disp%show(' ******** Integration did NOT converge. ********', tmsize = 1_IK, bmsize = 1_IK)
94 6 : call disp%show("getStr([truth, integral, abserr, abs(integral - truth)])//SK_' (unbiased)? '//getStr(abs(integral - truth) <= abserr)")
95 30 : call disp%show( getStr([truth, integral, abserr, abs(integral - truth)])//SK_' (unbiased)? '//getStr(abs(integral - truth) <= abserr) )
96 6 : call disp%show("numFuncEval")
97 6 : call disp%show( numFuncEval )
98 6 : call disp%skip()
99 :
100 18 : elseif (allocated(integrand%wcauchy)) then
101 :
102 2 : call disp%skip()
103 2 : call disp%show("! Assisted adaptive global quadrature by specifying Cauchy weight of the integrand.")
104 2 : call disp%skip()
105 :
106 2 : call disp%show("integrand%wcauchy%cs")
107 2 : call disp%show( integrand%wcauchy%cs )
108 2 : call disp%show("if (isFailedQuad(getFuncUnweighted, lb, ub, integrand%wcauchy, integral, abserr, neval = numFuncEval)) call disp%show(' ******** Integration did NOT converge. ********', tmsize = 1_IK, bmsize = 1_IK)")
109 2 : if (isFailedQuad(getFuncUnweighted, lb, ub, integrand%wcauchy, integral, abserr, neval = numFuncEval)) call disp%show(' ******** Integration did NOT converge. ********', tmsize = 1_IK, bmsize = 1_IK)
110 2 : call disp%show("getStr([truth, integral, abserr, abs(integral - truth)])//SK_' (unbiased)? '//getStr(abs(integral - truth) <= abserr)")
111 10 : call disp%show( getStr([truth, integral, abserr, abs(integral - truth)])//SK_' (unbiased)? '//getStr(abs(integral - truth) <= abserr) )
112 2 : call disp%show("numFuncEval")
113 2 : call disp%show( numFuncEval )
114 2 : call disp%skip()
115 :
116 : end if
117 :
118 : !%%%%%%%%%%%%%%%%%%%%%%
119 : #elif test_getQuadErr_ENABLED
120 : !%%%%%%%%%%%%%%%%%%%%%%
121 :
122 : integer(IK) :: nintmax_def
123 : integer(IK) :: err, numFuncEval, numInterval
124 : real(RKC) :: integral, abserr, abstol, reltol, lb, ub, truth
125 : real(RKC) , allocatable :: nodeK(:), weightK(:), weightG(:), sinfo(:,:), break(:)
126 : integer(IK) , allocatable :: sindex(:)
127 :
128 21 : abstol = getOption(0._RKC, atol)
129 21 : reltol = getOption(epsilon(0._RKC)**0.66, rtol)
130 21 : nintmax_def = getOption(2000_IK, nintmax)
131 21 : call setResized(sindex, nintmax_def)
132 63 : call setResized(sinfo, [4_IK, nintmax_def])
133 21 : truth = real(integrand%integral, RKC)
134 21 : lb = real(integrand%lb, RKC)
135 21 : ub = real(integrand%ub, RKC)
136 :
137 21 : call disp%skip()
138 21 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
139 21 : call disp%show("! "//integrand%desc)
140 21 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
141 21 : call disp%skip()
142 :
143 : #define DESC_INTEGRAND \
144 : call disp%show("integrand%desc"); \
145 : call disp%show( integrand%desc , deliml = SK_"""" ); \
146 : call disp%show("[abstol, reltol]"); \
147 : call disp%show( [abstol, reltol] ); \
148 : call disp%show("[lb, ub]"); \
149 : call disp%show( [lb, ub] ); \
150 : call disp%skip();
151 :
152 : #define DISP_INTEGRAL \
153 : call disp%show("if (err /= 0_IK) call disp%show(getStrUpper(SK_' **** integration failed: err = '//getStr(err)//SK_' ****'), tmsize = 1_IK, bmsize = 1_IK)"); \
154 : if (err /= 0_IK) call disp%show(getStrUpper(SK_' **** integration failed: err = '//getStr(err)//SK_' ****'), tmsize = 1_IK, bmsize = 1_IK); \
155 : call disp%show("getStr([truth, integral, abserr, abs(integral - truth)])//SK_' >= 1 (unbiased)? '//getStr(abs(integral - truth) <= abserr)"); \
156 : call disp%show( getStr([truth, integral, abserr, abs(integral - truth)])//SK_' >= 1 (unbiased)? '//getStr(abs(integral - truth) <= abserr) ); \
157 : call disp%show("[numFuncEval, numInterval]"); \
158 : call disp%show( [numFuncEval, numInterval] ); \
159 : call disp%skip();
160 :
161 : #define DISP_ARB_WEIGHT \
162 : call disp%show("call setNodeWeightGK(order = 35_IK, nodeK = nodeK, weightK = weightK, weightG = weightG) ! Compute the 35-71 Gauss-Kronrod nodes and weights."); \
163 : call setNodeWeightGK(order = 35_IK, nodeK = nodeK, weightK = weightK, weightG = weightG);
164 :
165 : ! Define the integration with arbitrary Gauss-Kronrod rule for a given integration range.
166 : #define TRIPLET_ARB(LBV,LBS,UBV,UBS) \
167 : DESC_INTEGRAND \
168 : DISP_ARB_WEIGHT \
169 : call disp%show("err = getQuadErr(getFunc, "//LBS//", "//UBS//", abstol, reltol, nodeK, weightK, weightG, integral, abserr, sinfo, sindex, numFuncEval, numInterval)"); \
170 : err = getQuadErr(getFunc, LBV, UBV , abstol, reltol, nodeK, weightK, weightG, integral, abserr, sinfo, sindex, numFuncEval, numInterval); \
171 : DISP_INTEGRAL \
172 : DESC_INTEGRAND \
173 : DISP_ARB_WEIGHT \
174 : call disp%show("err = getQuadErr(getFunc, "//LBS//", "//UBS//", abstol, reltol, nodeK, weightK, weightG, weps, integral, abserr, sinfo, sindex, numFuncEval, numInterval)"); \
175 : err = getQuadErr(getFunc, LBV, UBV , abstol, reltol, nodeK, weightK, weightG, weps, integral, abserr, sinfo, sindex, numFuncEval, numInterval); \
176 : DISP_INTEGRAL \
177 : if (allocated(integrand%break)) then; \
178 : DESC_INTEGRAND \
179 : call disp%show("break = integrand%break"); \
180 : break = integrand%break; \
181 : call disp%show("break"); \
182 : call disp%show( break ); \
183 : DISP_ARB_WEIGHT \
184 : call disp%show("err = getQuadErr(getFunc, "//LBS//", "//UBS//", abstol, reltol, nodeK, weightK, weightG, break, integral, abserr, sinfo, sindex, numFuncEval, numInterval)"); \
185 : err = getQuadErr(getFunc, LBV, UBV , abstol, reltol, nodeK, weightK, weightG, break, integral, abserr, sinfo, sindex, numFuncEval, numInterval); \
186 : DISP_INTEGRAL \
187 : elseif (allocated(integrand%wcauchy)) then; \
188 : DESC_INTEGRAND \
189 : call disp%show("integrand%wcauchy%cs"); \
190 : call disp%show( integrand%wcauchy%cs ); \
191 : DISP_ARB_WEIGHT \
192 : call disp%show("err = getQuadErr(getFuncUnweighted, "//LBS//", "//UBS//", abstol, reltol, nodeK, weightK, weightG, integrand%wcauchy, integral, abserr, sinfo, sindex, numFuncEval, numInterval)"); \
193 : err = getQuadErr(getFuncUnweighted, LBV, UBV , abstol, reltol, nodeK, weightK, weightG, integrand%wcauchy, integral, abserr, sinfo, sindex, numFuncEval, numInterval); \
194 : DISP_INTEGRAL \
195 : end if;
196 :
197 : ! Define the integration with predefined Gauss-Kronrod rule for a given integration range.
198 : #define TRIPLET_GKX(LBV,LBS,UBV,UBS,GKV,GKS) \
199 : DESC_INTEGRAND \
200 : DISP_ARB_WEIGHT \
201 : call disp%show("err = getQuadErr(getFunc, "//LBS//", "//UBS//", abstol, reltol, "//GKS//", integral, abserr, sinfo, sindex, numFuncEval, numInterval)"); \
202 : err = getQuadErr(getFunc, LBV, UBV , abstol, reltol, GKV , integral, abserr, sinfo, sindex, numFuncEval, numInterval); \
203 : DISP_INTEGRAL \
204 : DESC_INTEGRAND \
205 : DISP_ARB_WEIGHT \
206 : call disp%show("err = getQuadErr(getFunc, "//LBS//", "//UBS//", abstol, reltol, "//GKS//", weps, integral, abserr, sinfo, sindex, numFuncEval, numInterval)"); \
207 : err = getQuadErr(getFunc, LBV, UBV , abstol, reltol, GKV , weps, integral, abserr, sinfo, sindex, numFuncEval, numInterval); \
208 : DISP_INTEGRAL \
209 : if (allocated(integrand%break)) then; \
210 : DESC_INTEGRAND \
211 : call disp%show("break = integrand%break"); \
212 : break = integrand%break; \
213 : call disp%show("break"); \
214 : call disp%show( break ); \
215 : DISP_ARB_WEIGHT \
216 : call disp%show("err = getQuadErr(getFunc, "//LBS//", "//UBS//", abstol, reltol, "//GKS//", break, integral, abserr, sinfo, sindex, numFuncEval, numInterval)"); \
217 : err = getQuadErr(getFunc, LBV, UBV , abstol, reltol, GKV , break, integral, abserr, sinfo, sindex, numFuncEval, numInterval); \
218 : DISP_INTEGRAL \
219 : elseif (allocated(integrand%wcauchy)) then; \
220 : DESC_INTEGRAND \
221 : call disp%show("integrand%wcauchy%cs"); \
222 : call disp%show( integrand%wcauchy%cs ); \
223 : DISP_ARB_WEIGHT \
224 : call disp%show("err = getQuadErr(getFuncUnweighted, "//LBS//", "//UBS//", abstol, reltol, "//GKS//", integrand%wcauchy, integral, abserr, sinfo, sindex, numFuncEval, numInterval)"); \
225 : err = getQuadErr(getFuncUnweighted, LBV, UBV , abstol, reltol, GKV , integrand%wcauchy, integral, abserr, sinfo, sindex, numFuncEval, numInterval); \
226 : DISP_INTEGRAL \
227 : end if;
228 :
229 21 : if (getInfNeg(0._RKC) < lb .and. ub < getInfPos(0._RKC)) then
230 :
231 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
232 11 : call disp%show("! Using arbitrary Gauss-Kronrod nodes and weights (here: 35-71).")
233 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
234 11 : call disp%skip()
235 :
236 281 : TRIPLET_ARB(lb,"lb",ub,"ub");
237 :
238 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
239 11 : call disp%show("! Using predefined Gauss-Kronrod 30-61 nodes and weights.")
240 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
241 11 : call disp%skip()
242 :
243 279 : TRIPLET_GKX(lb,"lb",ub,"ub",GK61,"GK61")
244 :
245 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
246 11 : call disp%show("! Using predefined Gauss-Kronrod 25-51 nodes and weights.")
247 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
248 11 : call disp%skip()
249 :
250 279 : TRIPLET_GKX(lb,"lb",ub,"ub",GK51,"GK51")
251 :
252 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
253 11 : call disp%show("! Using predefined Gauss-Kronrod 20-41 nodes and weights.")
254 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
255 11 : call disp%skip()
256 :
257 279 : TRIPLET_GKX(lb,"lb",ub,"ub",GK41,"GK41")
258 :
259 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
260 11 : call disp%show("! Using predefined Gauss-Kronrod 15-31 nodes and weights.")
261 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
262 11 : call disp%skip()
263 :
264 279 : TRIPLET_GKX(lb,"lb",ub,"ub",GK31,"GK31")
265 :
266 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
267 11 : call disp%show("! Using predefined Gauss-Kronrod 10-21 nodes and weights.")
268 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
269 11 : call disp%skip()
270 :
271 279 : TRIPLET_GKX(lb,"lb",ub,"ub",GK21,"GK21")
272 :
273 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
274 11 : call disp%show("! Using predefined Gauss-Kronrod 7-15 nodes and weights.")
275 11 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
276 11 : call disp%skip()
277 :
278 279 : TRIPLET_GKX(lb,"lb",ub,"ub",GK15,"GK15")
279 :
280 10 : elseif (getInfNeg(0._RKC) < lb .and. huge(0._RKC) <= ub) then
281 :
282 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
283 4 : call disp%show("! Using arbitrary Gauss-Kronrod nodes and weights (here: 35-71).")
284 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
285 4 : call disp%skip()
286 :
287 98 : TRIPLET_ARB(lb,"lb",pinf,"pinf");
288 :
289 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
290 4 : call disp%show("! Using predefined Gauss-Kronrod 30-61 nodes and weights.")
291 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
292 4 : call disp%skip()
293 :
294 97 : TRIPLET_GKX(lb,"lb",pinf,"pinf",GK61,"GK61")
295 :
296 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
297 4 : call disp%show("! Using predefined Gauss-Kronrod 25-51 nodes and weights.")
298 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
299 4 : call disp%skip()
300 :
301 97 : TRIPLET_GKX(lb,"lb",pinf,"pinf",GK51,"GK51")
302 :
303 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
304 4 : call disp%show("! Using predefined Gauss-Kronrod 20-41 nodes and weights.")
305 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
306 4 : call disp%skip()
307 :
308 97 : TRIPLET_GKX(lb,"lb",pinf,"pinf",GK41,"GK41")
309 :
310 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
311 4 : call disp%show("! Using predefined Gauss-Kronrod 15-31 nodes and weights.")
312 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
313 4 : call disp%skip()
314 :
315 97 : TRIPLET_GKX(lb,"lb",pinf,"pinf",GK31,"GK31")
316 :
317 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
318 4 : call disp%show("! Using predefined Gauss-Kronrod 10-21 nodes and weights.")
319 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
320 4 : call disp%skip()
321 :
322 97 : TRIPLET_GKX(lb,"lb",pinf,"pinf",GK21,"GK21")
323 :
324 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
325 4 : call disp%show("! Using predefined Gauss-Kronrod 7-15 nodes and weights.")
326 4 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
327 4 : call disp%skip()
328 :
329 97 : TRIPLET_GKX(lb,"lb",pinf,"pinf",GK15,"GK15")
330 :
331 6 : elseif (lb <= -huge(0._RKC) .and. ub < getInfPos(0._RKC)) then
332 :
333 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
334 3 : call disp%show("! Using arbitrary Gauss-Kronrod nodes and weights (here: 35-71).")
335 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
336 3 : call disp%skip()
337 :
338 76 : TRIPLET_ARB(ninf,"ninf",ub,"ub");
339 :
340 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
341 3 : call disp%show("! Using predefined Gauss-Kronrod 30-61 nodes and weights.")
342 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
343 3 : call disp%skip()
344 :
345 76 : TRIPLET_GKX(ninf,"ninf",ub,"ub",GK61,"GK61")
346 :
347 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
348 3 : call disp%show("! Using predefined Gauss-Kronrod 25-51 nodes and weights.")
349 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
350 3 : call disp%skip()
351 :
352 76 : TRIPLET_GKX(ninf,"ninf",ub,"ub",GK51,"GK51")
353 :
354 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
355 3 : call disp%show("! Using predefined Gauss-Kronrod 20-41 nodes and weights.")
356 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
357 3 : call disp%skip()
358 :
359 76 : TRIPLET_GKX(ninf,"ninf",ub,"ub",GK41,"GK41")
360 :
361 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
362 3 : call disp%show("! Using predefined Gauss-Kronrod 15-31 nodes and weights.")
363 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
364 3 : call disp%skip()
365 :
366 76 : TRIPLET_GKX(ninf,"ninf",ub,"ub",GK31,"GK31")
367 :
368 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
369 3 : call disp%show("! Using predefined Gauss-Kronrod 10-21 nodes and weights.")
370 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
371 3 : call disp%skip()
372 :
373 76 : TRIPLET_GKX(ninf,"ninf",ub,"ub",GK21,"GK21")
374 :
375 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
376 3 : call disp%show("! Using predefined Gauss-Kronrod 7-15 nodes and weights.")
377 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
378 3 : call disp%skip()
379 :
380 76 : TRIPLET_GKX(ninf,"ninf",ub,"ub",GK15,"GK15")
381 :
382 3 : elseif (lb <= -huge(0._RKC) .and. huge(0._RKC) <= ub) then
383 :
384 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
385 3 : call disp%show("! Using arbitrary Gauss-Kronrod nodes and weights (here: 35-71).")
386 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
387 3 : call disp%skip()
388 :
389 96 : TRIPLET_ARB(ninf,"ninf",pinf,"pinf");
390 :
391 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
392 3 : call disp%show("! Using predefined Gauss-Kronrod 30-61 nodes and weights.")
393 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
394 3 : call disp%skip()
395 :
396 95 : TRIPLET_GKX(ninf,"ninf",pinf,"pinf",GK61,"GK61")
397 :
398 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
399 3 : call disp%show("! Using predefined Gauss-Kronrod 25-51 nodes and weights.")
400 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
401 3 : call disp%skip()
402 :
403 95 : TRIPLET_GKX(ninf,"ninf",pinf,"pinf",GK51,"GK51")
404 :
405 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
406 3 : call disp%show("! Using predefined Gauss-Kronrod 20-41 nodes and weights.")
407 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
408 3 : call disp%skip()
409 :
410 95 : TRIPLET_GKX(ninf,"ninf",pinf,"pinf",GK41,"GK41")
411 :
412 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
413 3 : call disp%show("! Using predefined Gauss-Kronrod 15-31 nodes and weights.")
414 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
415 3 : call disp%skip()
416 :
417 95 : TRIPLET_GKX(ninf,"ninf",pinf,"pinf",GK31,"GK31")
418 :
419 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
420 3 : call disp%show("! Using predefined Gauss-Kronrod 10-21 nodes and weights.")
421 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
422 3 : call disp%skip()
423 :
424 95 : TRIPLET_GKX(ninf,"ninf",pinf,"pinf",GK21,"GK21")
425 :
426 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
427 3 : call disp%show("! Using predefined Gauss-Kronrod 7-15 nodes and weights.")
428 3 : call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
429 3 : call disp%skip()
430 :
431 95 : TRIPLET_GKX(ninf,"ninf",pinf,"pinf",GK15,"GK15")
432 :
433 : end if
434 :
435 : #undef COMPUTE_INTEGRAL
436 : #undef DISP_ARB_WEIGHT
437 : #undef DESC_INTEGRAND
438 : #undef DISP_INTEGRAL
439 : #undef TRIPLET_GKX
440 : #undef TRIPLET_ARB
441 :
442 : #else
443 : !%%%%%%%%%%%%%%%%%%%%%%%%
444 : #error "Unrecognized interface."
445 : !%%%%%%%%%%%%%%%%%%%%%%%%
446 : #endif
447 :
448 : contains
449 :
450 4295490 : function getFunc(x) result(func)
451 : use pm_kind, only: RKH
452 : real(RKC) , intent(in) :: x
453 : real(RKC) :: func
454 4295490 : func = real(integrand%get(real(x, RKH)), RKC)
455 4295490 : if (allocated(integrand%wcauchy)) func = func / (x - real(integrand%wcauchy%cs, RKC))
456 4295490 : end function
457 :
458 6184 : function getFuncUnweighted(x) result(func)
459 : use pm_kind, only: RKH
460 : real(RKC) , intent(in) :: x
461 : real(RKC) :: func
462 6184 : func = real(integrand%get(real(x, RKH)), RKC)
463 6184 : end function
|
