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 the implementations of procedures of [pm_quadPack](@ref pm_quadPack).
19 : !>
20 : !> \finmain
21 : !>
22 : !> \author
23 : !> \AmirShahmoradi, Oct 16, 2009, 11:14 AM, Michigan
24 :
25 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26 :
27 : ! Define the help argument.
28 : #if QAGD_ENABLED
29 : #define HELP_ARG
30 : #elif QAGP_ENABLED && (IF_ENABLED || FI_ENABLED || II_ENABLED)
31 : #define HELP_ARG, breakTrans
32 : #elif QAGS_ENABLED || QAGP_ENABLED || QAWC_ENABLED || QAWFS_ENABLED || QAWFC_ENABLED
33 : #define HELP_ARG, help
34 : #elif !(constructWcauchy_ENABLED || constructWsin_ENABLED || constructWcos_ENABLED || getQuadGK_ENABLED || setNodeWeightGK_ENABLED || setErrSorted_ENABLED || setSeqLimEps_ENABLED || setChebExpan_ENABLED)
35 : #error "Unrecognized interface."
36 : #endif
37 : ! Define the GK quadrature arguments.
38 : #if GKXX_ENABLED
39 : #define QRULE_ARG nodeK, weightK, weightG
40 : #elif GK15_ENABLED || GK21_ENABLED || GK31_ENABLED || GK41_ENABLED || GK51_ENABLED || GK61_ENABLED
41 : #define QRULE_ARG qrule
42 : #elif !(constructWcauchy_ENABLED || constructWsin_ENABLED || constructWcos_ENABLED || setNodeWeightGK_ENABLED || setErrSorted_ENABLED || setSeqLimEps_ENABLED || setChebExpan_ENABLED || isFailedQuad_ENABLED)
43 : #error "Unrecognized interface."
44 : #endif
45 : !%%%%%%%%%%%%%%%%%%%%%%%
46 : #if constructWcauchy_ENABLED
47 : !%%%%%%%%%%%%%%%%%%%%%%%
48 :
49 4 : wcauchy%cs = real(cs, RKC)
50 :
51 : !%%%%%%%%%%%%%%%%%%%%
52 : #elif constructWsin_ENABLED
53 : !%%%%%%%%%%%%%%%%%%%%
54 :
55 0 : wsin%omega = real(omega, RKC)
56 :
57 : !%%%%%%%%%%%%%%%%%%%%
58 : #elif constructWcos_ENABLED
59 : !%%%%%%%%%%%%%%%%%%%%
60 :
61 0 : wcos%omega = real(omega, RKC)
62 :
63 : !%%%%%%%%%%%%%%%%%%%%%%
64 : #elif setNodeWeightGK_ENABLED
65 : !%%%%%%%%%%%%%%%%%%%%%%
66 :
67 : logical(LK) :: orderIsEven
68 : integer(IK) :: i, k, kk, l, ll, halfOrder, halfOrderPlusOne
69 : real(RKC) :: ak, an, bb, c, coefNode, coefWeight, d, x1, xx, s, y
70 : real(RKC) , parameter :: TOL = epsilon(0._RKC) * 10._RKC
71 : #if Fixed_ENABLED
72 252 : real(RKC) :: B((size(nodeK, kind = IK) / 2_IK) + 1_IK), Tau(size(nodeK, kind = IK) / 2_IK)
73 : integer(IK) :: order
74 126 : order = size(nodeK, kind = IK) - 1_IK
75 : #elif Alloc_ENABLED
76 700 : real(RKC) :: B(((order + 1_IK) / 2_IK) + 1_IK), Tau((order + 1_IK) / 2_IK)
77 3150 : allocate(nodeK(order + 1_IK), weightK(order + 1_IK), weightG((order + 1_IK) / 2_IK))
78 : #else
79 : #error "Unrecognized intefrace."
80 : #endif
81 476 : halfOrder = size(nodeK, kind = IK) / 2_IK
82 476 : halfOrderPlusOne = halfOrder + 1_IK
83 476 : orderIsEven = (2_IK * halfOrder == order)
84 :
85 476 : CHECK_ASSERTION(__LINE__, 0_IK <= order, SK_"@setNodeWeightGK(): The condition `0 <= order` must hold. order = "//getStr(order))
86 1428 : CHECK_ASSERTION(__LINE__, size(weightK, kind = IK) == order + 1_IK, SK_"@setNodeWeightGK(): The condition `size(weightK) == order + 1` must hold. size(weightK), order = "//getStr([size(weightK, kind = IK), order]))
87 1428 : CHECK_ASSERTION(__LINE__, size(weightG, kind = IK) == halfOrder, SK_"@setNodeWeightGK(): The condition `size(weightG) == halfOrder` must hold. size(weightG), halfOrder = "//getStr([size(weightG, kind = IK), halfOrder]))
88 :
89 476 : d = 2._RKC
90 476 : an = 0._RKC
91 14693 : do k = 1_IK, order
92 0 : an = an + 1._RKC
93 14693 : d = d * an / (an + 0.5_RKC)
94 : end do
95 :
96 : ! Compute the Chebyshev coefficients of the orthogonal polynomial.
97 :
98 476 : Tau(1) = (an + 2._RKC) / (an + an + 3._RKC)
99 0 : B(halfOrder) = Tau(1) - 1._RKC
100 476 : ak = an
101 7315 : do l = 1_IK, halfOrder - 1_IK
102 6839 : ak = ak + 2._RKC
103 6839 : Tau(l + 1_IK) = ((ak - 1._RKC) * ak - an * (an + 1._RKC)) * (ak + 2._RKC) * Tau(l) / (ak * ((ak + 3._RKC) * (ak + 2._RKC) - an * (an + 1._RKC)))
104 6839 : B(halfOrder - l) = Tau(l + 1_IK)
105 65617 : do ll = 1_IK, l
106 65141 : B(halfOrder - l) = B(halfOrder - l) + Tau(ll) * B(halfOrder - l + ll)
107 : end do
108 : end do
109 476 : B(halfOrderPlusOne) = 1._RKC
110 :
111 : ! Calculation of approximate values for the abscissas.
112 :
113 476 : bb = sin(1.570796_RKC / (an + an + 1._RKC))
114 476 : x1 = sqrt(1._RKC - bb * bb)
115 476 : s = 2._RKC * bb * x1
116 476 : c = sqrt(1._RKC - s * s)
117 476 : coefNode = 1._RKC - (1._RKC - 1._RKC / an) / (8._RKC * an * an)
118 476 : xx = coefNode * x1
119 :
120 : ! Coefficient needed for weights.
121 :
122 : ! COEF2 = 2^(2*order+1) * order! * order! / (2n+1)!
123 :
124 476 : coefWeight = 2._RKC / real(2_IK * order + 1_IK, RKC)
125 14693 : do i = 1_IK, order
126 14693 : coefWeight = coefWeight * 4._RKC * real(i, RKC) / real(order + i, RKC)
127 : end do
128 :
129 : ! Calculation of the K-th abscissa (a Kronrod abscissa) and the corresponding weight.
130 :
131 : !do k = 1_IK, order, 2_IK
132 7791 : do k = 1_IK, (order + 1_IK) / 2_IK
133 :
134 7315 : kk = 2_IK * k - 1_IK
135 7315 : call setNodeWeightK(xx, weightK(kk))
136 : !weightG(k) = 0._RKC
137 7315 : nodeK(kk) = xx
138 7315 : y = x1
139 7315 : x1 = y * c - bb * s
140 7315 : bb = y * s + bb * c
141 7315 : if (kk == order) then
142 413 : xx = 0._RKC
143 : else
144 6902 : xx = coefNode * x1
145 : end if
146 :
147 : ! Calculation of the K+1 abscissa (a Gaussian abscissa) and the corresponding weights.
148 :
149 7315 : call setNodeWeightGK(xx, weightK(kk + 1_IK), weightG(k))
150 :
151 7315 : nodeK(kk + 1_IK) = xx
152 7315 : y = x1
153 7315 : x1 = y * c - bb * s
154 7315 : bb = y * s + bb * c
155 7791 : xx = coefNode * x1
156 :
157 : end do
158 :
159 : ! Compute the Kronrod abscissa at the origin if order is even.
160 :
161 476 : if (orderIsEven) then
162 63 : xx = 0._RKC
163 63 : call setNodeWeightK(xx, weightK(order + 1_IK))
164 : !weightG(order + 1_IK) = 0._RKC
165 63 : nodeK(order + 1_IK) = xx
166 : end if
167 :
168 : contains
169 :
170 7378 : pure subroutine setNodeWeightK(node, weight)
171 :
172 : real(RKC) , intent(inout) :: node
173 : real(RKC) , intent(out) :: weight
174 : real(RKC) :: ai, b0, ub1, ub2, d0, d1, d2, delta, dif, getFunc, fd, yy
175 : logical(LK) :: convergenceOccurred
176 : integer(IK) :: ii, iter, kk
177 :
178 7378 : convergenceOccurred = node == 0._RKC
179 :
180 : ! Iteratively compute a Kronrod node.
181 :
182 35328 : do iter = 1_IK, 50_IK
183 :
184 35328 : ub1 = 0._RKC
185 35328 : ub2 = B(halfOrderPlusOne)
186 35328 : yy = 4._RKC * node * node - 2._RKC
187 35328 : d1 = 0._RKC
188 :
189 35328 : if (orderIsEven) then
190 2602 : ai = halfOrder + halfOrderPlusOne
191 2602 : d2 = ai * B(halfOrderPlusOne)
192 2602 : dif = 2._RKC
193 : else
194 32726 : ai = real(halfOrderPlusOne, RKC)
195 32726 : d2 = 0._RKC
196 32726 : dif = 1._RKC
197 : end if
198 :
199 630205 : do kk = 1_IK, halfOrder
200 594877 : ai = ai - dif
201 594877 : ii = halfOrderPlusOne - kk
202 594877 : b0 = ub1
203 594877 : ub1 = ub2
204 594877 : d0 = d1
205 594877 : d1 = d2
206 594877 : ub2 = yy * ub1 - b0 + B(ii)
207 594877 : if (.not. orderIsEven) ii = ii + 1_IK
208 630205 : d2 = yy * d1 - d0 + ai * B(ii)
209 : end do
210 :
211 35328 : if (orderIsEven) then
212 2602 : getFunc = node * (ub2 - ub1)
213 2602 : fd = d2 + d1
214 : else
215 32726 : getFunc = 0.5_RKC * (ub2 - b0)
216 32726 : fd = 4._RKC * node * d2
217 : end if
218 :
219 : ! Newton correction.
220 :
221 35328 : delta = getFunc / fd
222 35328 : node = node - delta
223 35328 : if (convergenceOccurred) then
224 7378 : d0 = 1._RKC
225 7378 : d1 = node
226 7378 : ai = 0._RKC
227 242053 : do kk = 2_IK, order
228 234675 : ai = ai + 1._RKC
229 234675 : d2 = ((ai + ai + 1._RKC) * node * d1 - ai * d0) / (ai + 1._RKC)
230 234675 : d0 = d1
231 242053 : d1 = d2
232 : end do
233 7378 : weight = coefWeight / (fd * d2)
234 7378 : return
235 : end if
236 27950 : convergenceOccurred = logical(abs(delta) <= TOL, LK)
237 :
238 : end do
239 :
240 0 : error stop "Iterative computation of the Kronrod node and weight failed."
241 :
242 : end subroutine
243 :
244 7315 : pure subroutine setNodeWeightGK(node, weightKronrod, weightGauss)
245 :
246 : real(RKC) , intent(inout) :: node
247 : real(RKC) , intent(out) :: weightKronrod, weightGauss
248 : real(RKC) :: ai, delta, p0, p1, p2, pd0, pd1, pd2, yy
249 : logical(LK) :: convergenceOccurred
250 : integer(IK) :: ii, iter, kk
251 7315 : convergenceOccurred = logical(node == 0._RKC, LK)
252 :
253 : ! Iteratively compute a Kronrod node.
254 :
255 35382 : do iter = 1_IK, 50_IK
256 :
257 35382 : p1 = node
258 35382 : p0 = 1._RKC
259 35382 : pd0 = 0._RKC
260 35382 : pd1 = 1._RKC
261 :
262 35382 : if (order <= 1) then ! Initialize P2 and PD2 to avoid problems with `delta`.
263 4 : if (epsilon(node) < abs(node)) then
264 0 : p2 = 0.5_RKC * (3._RKC * node * node - 1._RKC)
265 0 : pd2 = 3._RKC * node
266 : else
267 4 : p2 = 3._RKC * node
268 4 : pd2 = 3._RKC
269 : end if
270 : end if
271 :
272 35382 : ai = 0._RKC
273 1166419 : do kk = 2_IK, order
274 1131037 : ai = ai + 1._RKC
275 1131037 : p2 = ((ai + ai + 1._RKC) * node * p1 - ai * p0) / (ai + 1._RKC)
276 1131037 : pd2 = ((ai + ai + 1._RKC) * (p1 + node * pd1) - ai * pd0) / (ai + 1._RKC)
277 1131037 : p0 = p1
278 1131037 : p1 = p2
279 1131037 : pd0 = pd1
280 1166419 : pd1 = pd2
281 : end do
282 :
283 : ! Newton correction.
284 :
285 35382 : delta = p2 / pd2
286 35382 : node = node - delta
287 35382 : if (convergenceOccurred) then
288 7315 : weightGauss = 2._RKC / (order * pd2 * p0)
289 7315 : yy = 4._RKC * node * node - 2._RKC
290 7315 : p2 = B(halfOrderPlusOne)
291 7315 : p1 = 0._RKC
292 131234 : do kk = 1_IK, halfOrder
293 123919 : ii = halfOrderPlusOne - kk
294 123919 : p0 = p1
295 123919 : p1 = p2
296 131234 : p2 = yy * p1 - p0 + B(ii)
297 : end do
298 7315 : if (orderIsEven) then
299 510 : weightKronrod = weightGauss + coefWeight / (pd2 * node * (p2 - p1))
300 : else
301 6805 : weightKronrod = weightGauss + 2._RKC * coefWeight / (pd2 * (p2 - p0))
302 : end if
303 7315 : return
304 : end if
305 28067 : convergenceOccurred = logical(abs(delta) <= TOL, LK)
306 :
307 : end do
308 :
309 0 : error stop "Iterative computation of Gauss-Kronrod nodes and weights failed."
310 :
311 : end subroutine
312 :
313 : !%%%%%%%%%%%%%%%%%%%
314 : #elif setErrSorted_ENABLED
315 : !%%%%%%%%%%%%%%%%%%%
316 :
317 : ! set maxErrLoc and maxErrVal.
318 : #define SET_ERROR_VALUE maxErrLoc = sindex(nrmax); maxErrVal = sinfoErr(maxErrLoc);
319 :
320 : real(RKC) :: errmin
321 : integer(IK) :: i, ibeg, ido, isucc, j, jbnd, jupbn, k
322 : integer(IK) :: nint
323 :
324 171138 : nint = size(sindex, kind = IK)
325 513414 : CHECK_ASSERTION(__LINE__, nintmax >= nint, SK_"@setErrSorted(): The condition `nintmax >= nint` must hold. nintmax, nint = "//getStr([nintmax, nint]))
326 513414 : CHECK_ASSERTION(__LINE__, size(sinfoErr, kind = IK) == nint, SK_"@setErrSorted(): The condition `size(sinfoErr, kind = IK) == nint` must hold. size(sinfoErr, kind = IK), nint = "//getStr([size(sinfoErr, kind = IK), nint]))
327 :
328 : ! Check if the list contains more than two error estimates.
329 :
330 171138 : if (nint > 2_IK) then
331 :
332 : ! \remark
333 : ! This part of the routine is only executed if, due to a difficult integrand, subdivision increased the error estimate.
334 : ! The insert procedure should start after the nrmax-th largest error estimate, in the normal case.
335 :
336 131160 : maxErrVal = sinfoErr(maxErrLoc)
337 131160 : if (nrmax /= 1_IK) then
338 : ido = nrmax - 1_IK
339 20209 : do i = 1_IK, ido
340 20145 : isucc = sindex(nrmax - 1_IK)
341 20145 : if (maxErrVal <= sinfoErr(isucc)) exit
342 151 : sindex(nrmax) = isucc
343 20209 : nrmax = nrmax - 1_IK
344 : end do
345 : end if
346 :
347 : ! Compute the number of elements in the list to be maintained in descending order.
348 : ! This number depends on the number of subdivisions still allowed.
349 :
350 : jupbn = nint
351 131160 : if (nint > (nintmax / 2_IK + 2_IK)) jupbn = nintmax + 3_IK - nint
352 131160 : errmin = sinfoErr(nint)
353 :
354 : ! Insert maxErrVal by traversing the list top-down, starting comparison from the element sinfoErr(sindex(nrmax + 1_IK)).
355 :
356 131160 : jbnd = jupbn - 1_IK
357 131160 : ibeg = nrmax + 1_IK
358 131160 : if (ibeg <= jbnd) then
359 4235389 : do i = ibeg, jbnd
360 4228430 : isucc = sindex(i)
361 4228430 : if (maxErrVal >= sinfoErr(isucc)) then
362 : ! insert errmin by traversing the list bottom-up.
363 124166 : sindex(i - 1) = maxErrLoc
364 : k = jbnd
365 2571330 : do j = i, jbnd
366 2537103 : isucc = sindex(k)
367 2537103 : if (errmin < sinfoErr(isucc)) then
368 89939 : sindex(k + 1) = nint
369 89939 : SET_ERROR_VALUE ! fpp
370 89939 : return
371 : end if
372 2447164 : sindex(k + 1) = isucc
373 2481391 : k = k - 1_IK
374 : end do
375 34227 : sindex(i) = nint
376 34227 : SET_ERROR_VALUE ! fpp
377 34227 : return
378 : end if
379 4111223 : sindex(i - 1) = isucc
380 : end do
381 : end if
382 6994 : sindex(jbnd) = maxErrLoc
383 6994 : sindex(jupbn) = nint
384 : else
385 39978 : sindex(1) = 1
386 39978 : sindex(2) = 2
387 : end if
388 46972 : SET_ERROR_VALUE ! fpp
389 : #undef SET_ERROR_VALUE
390 :
391 : !%%%%%%%%%%%%%%%%%%%
392 : #elif setSeqLimEps_ENABLED
393 : !%%%%%%%%%%%%%%%%%%%
394 :
395 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
396 : real(RKC) , parameter :: EPS_RKC = epsilon(0._RKC)
397 : real(RKC) :: delta1, delta2, delta3, epsinf, err1, err2, err3, e0, e1, e1abs, e2, e3, res, ss, tol1, tol2, tol3
398 : real(RKC) :: error ! `error = abs(e1-e0)+abs(e2-e1)+abs(new-e2)`
399 : integer(IK) :: newelm ! number of elements to be computed in the new diagonal
400 : integer(IK) :: i, ib, ib2, indx, k1, k2, k3, num
401 : ! `seqlim` is the element in the new diagonal with least value of error.
402 24048 : ical = ical + 1_IK
403 24048 : abserr = HUGE_RKC
404 24048 : seqlim = EpsTable(inew)
405 24048 : blockNewCondition: if (inew >= 3_IK) then
406 24048 : EpsTable(inew + 2) = EpsTable(inew)
407 24048 : newelm = (inew - 1_IK) / 2_IK
408 24048 : EpsTable(inew) = HUGE_RKC
409 : num = inew
410 : k1 = inew
411 110064 : loopNewElement: do i = 1_IK, newelm
412 86620 : k2 = k1 - 1_IK
413 86620 : k3 = k1 - 2_IK
414 86620 : res = EpsTable(k1 + 2_IK)
415 86620 : e0 = EpsTable(k3)
416 86620 : e1 = EpsTable(k2)
417 70224 : e2 = res
418 86620 : e1abs = abs(e1)
419 86620 : delta2 = e2 - e1
420 86620 : err2 = abs(delta2)
421 86620 : tol2 = max(abs(e2), e1abs) * EPS_RKC
422 86620 : delta3 = e1 - e0
423 86620 : err3 = abs(delta3)
424 86620 : tol3 = max(e1abs, abs(e0)) * EPS_RKC
425 110064 : if (err2 > tol2 .or. err3 > tol3) then
426 86528 : e3 = EpsTable(k1)
427 86528 : EpsTable(k1) = e1
428 86528 : delta1 = e1 - e3
429 86528 : err1 = abs(delta1)
430 86528 : tol1 = max(e1abs, abs(e3)) * EPS_RKC
431 : ! If two elements are very close to each other, omit a part of the table by adjusting the value of `inew`.
432 86528 : if (err1 > tol1 .and. err2 > tol2 .and. err3 > tol3) then
433 86016 : ss = 1._RKC / delta1 + 1._RKC / delta2 - 1._RKC / delta3
434 86016 : epsinf = abs(ss * e1)
435 : ! Test to detect irregular behavior in the table and eventually omit a part of the table adjusting the value of `inew`.
436 86016 : if (epsinf > 1.e-4_RKC) then
437 86016 : res = e1 + 1._RKC / ss
438 86016 : EpsTable(k1) = res
439 : k1 = k1 - 2_IK
440 86016 : error = err2 + abs(res - e2) + err3
441 86016 : if (error <= abserr) then
442 34327 : abserr = error
443 34327 : seqlim = res
444 : end if
445 : cycle loopNewElement
446 : end if
447 : end if
448 512 : inew = i + i - 1_IK
449 512 : exit loopNewElement
450 : else ! e0, e1 and e2 are equal to within machine accuracy. Convergence is assumed.
451 : ! abserr = abs(e1 - e0) + abs(e2 - e1)
452 : ! seqlim = e2
453 92 : seqlim = res
454 92 : abserr = err2 + err3
455 92 : exit blockNewCondition
456 : end if
457 : end do loopNewElement
458 : ! shift the table.
459 : !if (inew > MAXLEN_EPSTAB) error stop "inew > MAXLEN_EPSTAB cannot happen."
460 23956 : if (inew == MAXLEN_EPSTAB) inew = 2_IK * (MAXLEN_EPSTAB / 2_IK) - 1_IK
461 23956 : ib = 1_IK; if ((num / 2_IK) * 2_IK == num) ib = 2_IK
462 134747 : do i = 1_IK, newelm + 1_IK
463 110791 : ib2 = ib + 2_IK
464 110791 : EpsTable(ib) = EpsTable(ib2)
465 23956 : ib = ib2
466 : end do
467 23956 : if (num /= inew) then
468 1321 : indx = num - inew + 1_IK
469 42628 : do i = 1_IK, inew
470 41307 : EpsTable(i) = EpsTable(indx)
471 42628 : indx = indx + 1_IK
472 : end do
473 : end if
474 23956 : if (ical < 4_IK) then
475 12961 : seqlims(ical) = seqlim
476 12961 : abserr = HUGE_RKC
477 : else
478 43980 : abserr = sum(abs(seqlim - seqlims)) ! abs(seqlim - seqlims(3)) + abs(seqlim - seqlims(2)) + abs(seqlim - seqlims(1))
479 10995 : seqlims(1) = seqlims(2)
480 10995 : seqlims(2) = seqlims(3)
481 10995 : seqlims(3) = seqlim
482 : end if
483 : end if blockNewCondition
484 24048 : abserr = max(abserr, 5._RKC * EPS_RKC * abs(seqlim)) ! \warning This is set incorrectly in QuadPack of John Burkardt.
485 :
486 : !%%%%%%%%%%%%%%%%%%%
487 : #elif isFailedQuad_ENABLED
488 : !%%%%%%%%%%%%%%%%%%%
489 :
490 : real(RKC), parameter :: vinf = huge(0._RKC) / 1000._RKC ! virtual infinity.
491 : integer(IK) , parameter :: MAX_NUM_INTERVAL = 2000_IK
492 : integer(IK) :: err, nint_def, neval_def, sindex(MAX_NUM_INTERVAL)
493 : real(RKC) :: abserr_def, abstol_def, reltol_def
494 : real(RKC) :: sinfo(4, MAX_NUM_INTERVAL)
495 5843 : abstol_def = 0._RKC; if (present(abstol)) abstol_def = abstol
496 5843 : reltol_def = epsilon(0._RKC)**(2./3.); if (present(reltol)) reltol_def = reltol
497 5843 : if (-vinf < lb .and. ub < vinf) then ! finite
498 5733 : err = getQuadErr(getFunc, lb, ub, abstol_def, reltol_def, GK21 HELP_ARG, integral, abserr_def, sinfo, sindex, neval_def, nint_def)
499 110 : elseif (lb <= -vinf .and. vinf <= ub) then ! both infinite
500 12 : err = getQuadErr(getFunc, ninf, pinf, abstol_def, reltol_def, GK21 HELP_ARG, integral, abserr_def, sinfo, sindex, neval_def, nint_def)
501 98 : elseif (vinf <= ub) then ! positive semi-infinite
502 91 : err = getQuadErr(getFunc, lb, pinf, abstol_def, reltol_def, GK21 HELP_ARG, integral, abserr_def, sinfo, sindex, neval_def, nint_def)
503 7 : elseif (lb <= -vinf) then ! negative semi-infinite
504 7 : err = getQuadErr(getFunc, ninf, ub, abstol_def, reltol_def, GK21 HELP_ARG, integral, abserr_def, sinfo, sindex, neval_def, nint_def)
505 : else ! possibly NAN
506 : failed = .true._LK
507 0 : if (present(msg)) msg = SK_"The specified integration range is invalid. lb, ub = "//getStr([lb, ub])
508 0 : return
509 : end if
510 5843 : failed = logical(err /= 0_IK, kind = LK)
511 5843 : if (present(abserr)) abserr = abserr_def
512 5843 : if (present(neval)) neval = neval_def
513 5843 : if (present(nint)) nint = nint_def
514 5843 : if (failed) then
515 10 : if (present(msg)) then
516 : if (err == 1_IK) then
517 0 : msg = SK_"The maximum allowed number of subdivisions was reached without convergence."
518 : elseif (err == 2_IK) then
519 0 : msg = SK_"Roundoff error was detected, which prevents the requested tolerance from being achieved."
520 : elseif (err == 3_IK) then
521 0 : msg = SK_"Roundoff error was detected, which prevents the requested tolerance from being achieved."
522 : elseif (err == 4_IK) then
523 0 : msg = SK_"An extremely bad integrand behavior occurs at some points of the integration interval."
524 : elseif (err == 5_IK) then
525 0 : msg = SK_"The algorithm did not converge because the integral is likely divergent or slowly convergent."
526 : else
527 0 : error stop SK_"Internal error occurred. Invalid error code detected: "//getStr(err)
528 : end if
529 : end if
530 : end if
531 :
532 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
533 : #elif getQuadErr_ENABLED && (QAGD_ENABLED || QAGS_ENABLED || QAGP_ENABLED || QAWC_ENABLED) && (FI_ENABLED || IF_ENABLED || II_ENABLED)
534 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
535 :
536 : ! Setup the transformation coefficient.
537 : #if QAWC_ENABLED
538 : real(RKC) :: cstrans
539 : #if II_ENABLED
540 : #define TRANS_COEF invx * (-cstrans)
541 0 : cstrans = 1._RKC / (1._RKC + real(help%cs, RKC))
542 : #elif IF_ENABLED
543 : #define TRANS_COEF invx * cstrans
544 0 : cstrans = 1._RKC / (1._RKC - real(help%cs, RKC) + ub)
545 : #elif FI_ENABLED
546 : #define TRANS_COEF invx * (-cstrans)
547 0 : cstrans = 1._RKC / (1._RKC + real(help%cs, RKC) - lb)
548 : #endif
549 : #elif QAGD_ENABLED || QAGS_ENABLED || QAGP_ENABLED
550 : #define TRANS_COEF invx**2
551 : #else
552 : #error "Unrecognized interface."
553 : #endif
554 : ! Transform QAGP break points.
555 : #if QAGP_ENABLED && II_ENABLED
556 16 : real(RKC) :: breakTrans(size(help, kind = IK))
557 : ! Double infinite range requires special care to remove duplicate breaks.
558 : integer(IK) :: locLastNegBreak, lenBreak
559 96 : CHECK_ASSERTION(__LINE__, all(help /= -1._RKC), SK_"@getQuadErr(): The condition `all(help /= -1.)` must hold. help = "//getStr(help))
560 : lenBreak = size(help, kind = IK)
561 128 : locLastNegBreak = maxloc(help, mask = help < -1._RKC, dim = 1)
562 96 : call setMerged(breakTrans, 1._RKC / (1._RKC - help(1:locLastNegBreak)), 1._RKC / (1._RKC + help(lenBreak : locLastNegBreak + 1 : -1)))
563 : !print *, help(1:locLastNegBreak); print *, help(locLastNegBreak+1:); print *, breakTrans
564 : #elif QAGP_ENABLED && IF_ENABLED
565 16 : real(RKC) :: breakTrans(size(help, kind = IK))
566 24 : breakTrans = 1._RKC / (1._RKC - help + ub)
567 : #elif QAGP_ENABLED && FI_ENABLED
568 16 : real(RKC) :: breakTrans(size(help, kind = IK))
569 24 : breakTrans = 1._RKC / (1._RKC + help(size(help) : 1 : -1) - lb)
570 : #elif !(QAGD_ENABLED || QAGS_ENABLED || QAWC_ENABLED || QAWFS_ENABLED || QAWFC_ENABLED)
571 : #error "Unrecognized interface."
572 : #endif
573 1786 : err = getQuadErr(getFuncTrans, 0._RKC, 1._RKC, abstol, reltol, QRULE_ARG HELP_ARG, integral, abserr, sinfo, sindex, neval, nint)
574 : contains
575 2216674 : function getFuncTrans(x) result(func)
576 : implicit none
577 : real(RKC), intent(in) :: x
578 : real(RKC) :: func
579 : real(RKC) :: invx
580 : real(RKC) :: xMinuxOne2x
581 2216674 : invx = 1._RKC / x
582 2216674 : xMinuxOne2x = (x - 1._RKC) * invx
583 : #if II_ENABLED
584 963036 : func = (getFunc(xMinuxOne2x) + getFunc(-xMinuxOne2x)) * TRANS_COEF
585 : #elif IF_ENABLED
586 312262 : func = getFunc(ub + xMinuxOne2x) * TRANS_COEF
587 : #elif FI_ENABLED
588 941376 : func = getFunc(lb - xMinuxOne2x) * TRANS_COEF
589 : #else
590 : #error "Unrecognized interface."
591 : #endif
592 2216674 : end function
593 :
594 : !%%%%%%%%%%%%%%%%
595 : #elif getQuadGK_ENABLED
596 : !%%%%%%%%%%%%%%%%
597 :
598 : #if GKXX_ENABLED
599 : #define SET_PARAMETER(x) x
600 : #elif GK15_ENABLED || GK21_ENABLED || GK31_ENABLED || GK41_ENABLED || GK51_ENABLED || GK61_ENABLED
601 : #define SET_PARAMETER(x) parameter(x)
602 : #else
603 : #error "Unrecognized interface."
604 : #endif
605 : ! Define the nodes for different rules.
606 : #if GK15_ENABLED
607 : !real(RKC) , parameter :: nodeG(*) = real(nodeG7(size(nodeG7):1:-1), RKC)
608 : real(RKC) , parameter :: nodeK(*) = real(nodeK15(size(nodeK15):1:-1), RKC)
609 : real(RKC) , parameter :: weightG(*) = real(weightG7(size(weightG7):1:-1), RKC)
610 : real(RKC) , parameter :: weightK(*) = real(weightK15(size(weightK15):1:-1), RKC)
611 : #elif GK21_ENABLED
612 : !real(RKC) , parameter :: nodeG(*) = real(nodeG10(size(nodeG10):1:-1), RKC)
613 : real(RKC) , parameter :: nodeK(*) = real(nodeK21(size(nodeK21):1:-1), RKC)
614 : real(RKC) , parameter :: weightG(*) = real(weightG10(size(weightG10):1:-1), RKC)
615 : real(RKC) , parameter :: weightK(*) = real(weightK21(size(weightK21):1:-1), RKC)
616 : #elif GK31_ENABLED
617 : !real(RKC) , parameter :: nodeG(*) = real(nodeG15(size(nodeG15):1:-1), RKC)
618 : real(RKC) , parameter :: nodeK(*) = real(nodeK31(size(nodeK31):1:-1), RKC)
619 : real(RKC) , parameter :: weightG(*) = real(weightG15(size(weightG15):1:-1), RKC)
620 : real(RKC) , parameter :: weightK(*) = real(weightK31(size(weightK31):1:-1), RKC)
621 : #elif GK41_ENABLED
622 : !real(RKC) , parameter :: nodeG(*) = real(nodeG20(size(nodeG20):1:-1), RKC)
623 : real(RKC) , parameter :: nodeK(*) = real(nodeK41(size(nodeK41):1:-1), RKC)
624 : real(RKC) , parameter :: weightG(*) = real(weightG20(size(weightG20):1:-1), RKC)
625 : real(RKC) , parameter :: weightK(*) = real(weightK41(size(weightK41):1:-1), RKC)
626 : #elif GK51_ENABLED
627 : !real(RKC) , parameter :: nodeG(*) = real(nodeG25(size(nodeG25):1:-1), RKC)
628 : real(RKC) , parameter :: nodeK(*) = real(nodeK51(size(nodeK51):1:-1), RKC)
629 : real(RKC) , parameter :: weightG(*) = real(weightG25(size(weightG25):1:-1), RKC)
630 : real(RKC) , parameter :: weightK(*) = real(weightK51(size(weightK51):1:-1), RKC)
631 : #elif GK61_ENABLED
632 : !real(RKC) , parameter :: nodeG(*) = real(nodeG30(size(nodeG30):1:-1), RKC)
633 : real(RKC) , parameter :: nodeK(*) = real(nodeK61(size(nodeK61):1:-1), RKC)
634 : real(RKC) , parameter :: weightG(*) = real(weightG30(size(weightG30):1:-1), RKC)
635 : real(RKC) , parameter :: weightK(*) = real(weightK61(size(weightK61):1:-1), RKC)
636 : #elif !GKXX_ENABLED
637 : #error "Unrecognized interface."
638 : #endif
639 : real(RKC) , parameter :: TINY_RKC = tiny(0._RKC)
640 : real(RKC) , parameter :: EPS_RKC = epsilon(0._RKC)
641 : integer(IK) :: i, ishared, sizeNodeG
642 : real(RKC) :: sumFunc
643 : real(RKC) :: midFunc ! function value at the center of the interval.
644 : real(RKC) :: avgFunc ! approximation to the mean value of `getFunc` over `(lb, ub)`, i.e. to `quadGK / (ub - lb)`
645 : real(RKC) :: quadG ! integral via the 7-point gauss formula
646 : real(RKC) :: quadK ! integral via the 15-point kronrod formula
647 : real(RKC) :: abscissa
648 : real(RKC) :: midInterval
649 : real(RKC) :: halfInterval
650 16428 : real(RKC) :: absHalfInterval, NegFunc(size(nodeK, kind = IK) - 1_IK), PosFunc(size(nodeK, kind = IK) - 1_IK)
651 : integer(IK) :: offset
652 :
653 : ! integration limits.
654 : #if FF_ENABLED
655 : #define LBT lb
656 : #define UBT ub
657 : #else
658 : real(RKC) :: invx, xMinuxOne2x
659 : #define LBT 0._RKC
660 : #define UBT 1._RKC
661 : #endif
662 : ! integration function.
663 : #if FF_ENABLED
664 : #define EVAL_EXPR(y, x) y = getFunc(x)
665 : #elif FI_ENABLED
666 : #define EVAL_EXPR(y, x) invx = 1._RKC / (x); xMinuxOne2x = (x - 1._RKC) * invx; y = getFunc(lb - xMinuxOne2x) * invx**2
667 : #elif IF_ENABLED
668 : #define EVAL_EXPR(y, x) invx = 1._RKC / (x); xMinuxOne2x = (x - 1._RKC) * invx; y = getFunc(ub + xMinuxOne2x) * invx**2
669 : #elif II_ENABLED
670 : #define EVAL_EXPR(y, x) invx = 1._RKC / (x); xMinuxOne2x = (x - 1._RKC) * invx; y = (getFunc(xMinuxOne2x) + getFunc(-xMinuxOne2x)) * invx**2
671 : #else
672 : #error "Unrecognized interface."
673 : #endif
674 8214 : SET_PARAMETER(sizeNodeG = int(0.5 * real(size(nodeK, kind = IK)), IK))
675 8214 : SET_PARAMETER(offset = merge(1_IK, 0_IK, mod(size(nodeK, kind = IK) - 1_IK, 2_IK) == 1_IK)) ! is 1 when the Gauss abscissas includes the origin 0.
676 : #if GKXX_ENABLED
677 24642 : CHECK_ASSERTION(__LINE__, size(weightK) == size(nodeK), SK_"@getQuadGK(): The condition `size(weightK) == size(nodeK)` must hold. size(weightK), size(nodeK) = "//getStr([size(weightK), size(nodeK)]))
678 24642 : CHECK_ASSERTION(__LINE__, size(weightG) == size(nodeK) / 2, SK_"@getQuadGK(): The condition `size(weightG) == size(nodeK)/2` must hold. size(weightG), size(nodeK)/2 = "//getStr([size(weightG), size(nodeK)/2]))
679 : #else
680 425776 : CHECK_ASSERTION(__LINE__, precision(nodeK) <= 100, SK_"@getQuadGK(): The condition `precision(nodeK) <= 100` must hold. precision(nodeK) = "//getStr(precision(nodeK)))
681 : #endif
682 433990 : midInterval = 0.5_RKC * (LBT + UBT)
683 433990 : halfInterval = 0.5_RKC * (UBT - LBT)
684 433990 : absHalfInterval = abs(halfInterval)
685 :
686 : ! Initialize the summation depending on whether zero is in the abscissas of the Gauss summation.
687 :
688 433990 : EVAL_EXPR(midFunc, midInterval) ! fpp
689 8214 : if (offset == 1_IK) then
690 8154 : quadG = midFunc * weightG(size(weightG))
691 : else
692 60 : quadG = 0._RKC
693 : end if
694 433990 : quadK = midFunc * weightK(size(weightK))
695 425776 : intAbsFunc = abs(quadK)
696 :
697 : ! Gauss summation.
698 :
699 2871370 : do i = 1_IK, sizeNodeG - offset
700 2437380 : ishared = 2_IK * i
701 2437380 : abscissa = halfInterval * nodeK(ishared)
702 2437380 : EVAL_EXPR(NegFunc(ishared), midInterval - abscissa) ! fpp
703 2437380 : EVAL_EXPR(PosFunc(ishared), midInterval + abscissa) ! fpp
704 2437380 : sumFunc = NegFunc(ishared) + PosFunc(ishared)
705 2437380 : quadG = quadG + sumFunc * weightG(i)
706 2437380 : quadK = quadK + sumFunc * weightK(ishared)
707 2871370 : intAbsFunc = intAbsFunc + weightK(ishared) * (abs(NegFunc(ishared)) + abs(PosFunc(ishared)))
708 : end do
709 :
710 : ! Kronrod summation.
711 :
712 8214 : do i = 1_IK, 2_IK * sizeNodeG - 1_IK, 2_IK
713 2496593 : abscissa = halfInterval * nodeK(i)
714 2496593 : EVAL_EXPR(NegFunc(i), midInterval - abscissa) ! fpp
715 2496593 : EVAL_EXPR(PosFunc(i), midInterval + abscissa) ! fpp
716 2496593 : sumFunc = NegFunc(i) + PosFunc(i)
717 2496593 : quadK = quadK + weightK(i) * sumFunc
718 2496593 : intAbsFunc = intAbsFunc + weightK(i) * (abs(NegFunc(i)) + abs(PosFunc(i)))
719 : end do
720 :
721 : ! Estimate the error.
722 :
723 433990 : avgFunc = quadK * 0.5_RKC
724 433990 : smoothness = weightK(1) * abs(midFunc - avgFunc)
725 5367963 : do i = 1_IK, size(nodeK, kind = IK) - 1_IK
726 5367963 : smoothness = smoothness + weightK(i) * (abs(NegFunc(i) - avgFunc) + abs(PosFunc(i) - avgFunc))
727 : end do
728 433990 : quadGK = quadK * halfInterval
729 433990 : intAbsFunc = intAbsFunc*absHalfInterval
730 433990 : smoothness = smoothness*absHalfInterval
731 433990 : abserr = abs((quadK - quadG) * halfInterval)
732 : !write(*,*) "======================================================================================"
733 : !print *, "nodeK", nodeK
734 : !print *, "weightK", weightK
735 : !print *, "weightG", weightG
736 : !print *, "size(nodeK), size(weightK), size(weightG)", size(nodeK), size(weightK), size(weightG)
737 : !print *, "abserr, quadK, quadG", abserr, quadK, quadG, halfInterval
738 : !write(*,*) "======================================================================================"
739 : !write(*,*)
740 433990 : if (smoothness /= 0._RKC .and. abserr /= 0._RKC) then
741 390454 : abserr = 200._RKC * abserr / smoothness ! `abserr` represents relative error here.
742 390454 : if (abserr < 1._RKC) then
743 336694 : abserr = smoothness * sqrt(abserr)**3
744 : else
745 53760 : abserr = smoothness
746 : end if
747 : !abserr = smoothness * min(1._RKC, sqrt(200._RKC * abserr / smoothness)**3) ! the sqrt and exponentiation can be avoided here.
748 : end if
749 433990 : if (intAbsFunc > TINY_RKC / (50._RKC * EPS_RKC)) abserr = max((EPS_RKC * 50._RKC) * intAbsFunc, abserr)
750 : #undef SET_PARAMETER
751 : #undef EVAL_EXPR
752 : #undef LBT
753 : #undef UBT
754 :
755 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
756 : #elif getQuadErr_ENABLED && (QAGD_ENABLED || QAGS_ENABLED || QAGP_ENABLED || QAWC_ENABLED) && FF_ENABLED
757 : !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
758 :
759 : #if QAWC_ENABLED
760 : real(RKC) :: cs
761 : #endif
762 : real(RKC) , parameter :: HUGE_RKC = huge(0._RKC)
763 : real(RKC) , parameter :: TINY_RKC = tiny(0._RKC)
764 : real(RKC) , parameter :: EPS_RKC = epsilon(0._RKC)
765 : integer(IK) :: iroff1, iroff2
766 : integer(IK) :: nintmax, nrmax, maxErrLoc ! pointer to the interval with largest error estimate.
767 : real(RKC) :: result1, lb1, ub1, abserr1 ! variables for the left subinterval.
768 : real(RKC) :: result2, lb2, ub2, abserr2 ! variables for the right subinterval.
769 : real(RKC) :: result ! sum of the integrals over the subintervals
770 : real(RKC) :: result12 ! `result1 + result2`
771 : real(RKC) :: abserr12 ! `abserr1 + abserr2`
772 : real(RKC) :: maxErrVal ! `elist(maxErrLoc)`.
773 : real(RKC) :: errsum ! sum of the errors over the subintervals.
774 : real(RKC) :: errbnd ! requested accuracy `max(abstol, reltol * abs(integral))`.
775 : real(RKC) :: intAbsFunc1, intAbsFunc2
776 : real(RKC) :: smoothness1, smoothness2
777 : #if QAGD_ENABLED || QAWC_ENABLED
778 : #define FINALIZE_INTEGRATION integral = sum(sinfo(3,1:nint)); abserr = errsum; return
779 : real(RKC) , parameter :: FACTOR = 1._RKC
780 : #elif QAGS_ENABLED || QAGP_ENABLED
781 : #define FINALIZE_INTEGRATION integral = sum(sinfo(3,1:nint)); abserr = errsum; if (err > 2_IK) err = err - 1_IK; return
782 : real(RKC) , parameter :: FACTOR = 2._RKC
783 : integer(IK) :: ksgn
784 : real(RKC) :: EpsTable(MAXLEN_EPSTAB + 2)
785 : real(RKC) :: absErrEps, seqlim, seqlims(3)
786 : ! error on the interval currently subdivided;
787 : ! sum of the errors over the intervals larger than the smallest interval considered up to now.
788 : real(RKC) :: erlast, erlarg, correc, ertest
789 : integer(IK) :: lenEpsTable, nextrap, iroffnew, k, ktmin, jupbnd, ierro ! roundoff error flag
790 : logical(LK) :: extrapolating, extrapAllowed
791 : #if QAGS_ENABLED
792 : #define CYCLE_NEEDED abs(sinfo(2,maxErrLoc) - sinfo(1,maxErrLoc)) > small
793 : #define EPS_TABLE_START_OFFSET 1_IK
794 : real(RKC) :: small
795 : #elif QAGP_ENABLED
796 : #define CYCLE_NEEDED level(maxErrLoc) < maximumLevel
797 : #define EPS_TABLE_START_OFFSET 0_IK
798 : integer(IK) :: i, j, jlow, ind1, ind2, lenBreak, lenBreakPlusOne
799 96 : integer(IK) :: currentLevel, maximumLevel, level(size(sindex, kind = IK))
800 96 : logical(LK) :: refinementNeeded(size(help, kind = IK) + 1_IK)
801 : #else
802 : #error "Unrecognized interface."
803 : #endif
804 : #else
805 : #error "Unrecognized interface."
806 : #endif
807 : #if GK15_ENABLED
808 : integer(IK) , parameter :: NEVAL0 = 15_IK
809 : #elif GK21_ENABLED
810 : integer(IK) , parameter :: NEVAL0 = 21_IK
811 : #elif GK31_ENABLED
812 : integer(IK) , parameter :: NEVAL0 = 31_IK
813 : #elif GK41_ENABLED
814 : integer(IK) , parameter :: NEVAL0 = 41_IK
815 : #elif GK51_ENABLED
816 : integer(IK) , parameter :: NEVAL0 = 51_IK
817 : #elif GK61_ENABLED
818 : integer(IK) , parameter :: NEVAL0 = 61_IK
819 : #elif GKXX_ENABLED
820 : integer(IK) :: NEVAL0
821 50 : NEVAL0 = 2_IK * size(nodeK, kind = IK) - 1_IK
822 : #else
823 : #error "Unrecognized interface."
824 : #endif
825 91494 : nintmax = size(sindex, kind = IK)
826 : #if QAWC_ENABLED
827 16 : cs = real(help%cs, RKC)
828 48 : CHECK_ASSERTION(__LINE__, lb < cs, SK_"@getQuadErr(): The condition `lb < cs` must hold. lb, cs = "//getStr([lb, cs]))
829 48 : CHECK_ASSERTION(__LINE__, cs < ub, SK_"@getQuadErr(): The condition `cs < ub` must hold. cs, ub = "//getStr([cs, ub]))
830 : #endif
831 : ! \bug Intel ifort 2021.6
832 : ! The intel ifort yields an internal compiler error which appears to be related to the following macro placements.
833 : ! As such they are commented out and implemented explicitly as all-in-one checking block.
834 : ! Update Dec 2022: This bug is likely related to the Intel bug with the maximum number
835 : ! of `use` statements in subdmoules in [pm_araryRebill](@ref pm_araryRebill).
836 274482 : CHECK_ASSERTION(__LINE__, lb < ub, SK_"@getQuadErr(): The condition `lb < ub` must hold. lb, ub = "//getStr([lb, ub]))
837 91494 : CHECK_ASSERTION(__LINE__, abstol >= 0._RKC, SK_"@getQuadErr(): The condition `abstol >= 0.` must hold. abstol = "//getStr(abstol))
838 91494 : CHECK_ASSERTION(__LINE__, reltol >= 0._RKC, SK_"@getQuadErr(): The condition `reltol >= 0.` must hold. reltol = "//getStr(reltol))
839 731952 : CHECK_ASSERTION(__LINE__, all(shape(sinfo, kind = IK) == [4_IK, nintmax]), SK_"@getQuadErr(): The condition `all(shape(sinfo) == [4, size(sindex)])` must hold. shape(sinfo), size(sindex) = "//getStr([shape(sinfo, kind = IK), size(sindex, kind = IK)]))
840 274482 : CHECK_ASSERTION(__LINE__, size(sinfo, 2, IK) == nintmax, SK_"@getQuadErr(): The condition `size(sinfo,2) == size(sindex)` must hold. size(sinfo,2), size(sindex) = "//getStr([size(sinfo, 2, IK), nintmax]))
841 91494 : CHECK_ASSERTION(__LINE__, nintmax >= 1_IK, SK_"@getQuadErr(): The condition `size(sindex, kind = IK) >= 1` must hold. size(sinfo,2) = "//getStr(size(sindex, kind = IK)))
842 : #if QAGP_ENABLED
843 48 : neval = 0_IK
844 : lenBreak = size(help, kind = IK)
845 : lenBreakPlusOne = lenBreak + 1_IK
846 48 : CHECK_ASSERTION(__LINE__, size(help, kind = IK) > 0_IK, SK_"@getQuadErr(): The condition `size(help) > 0` must hold. size(help) = "//getStr(size(help, kind = IK)))
847 144 : CHECK_ASSERTION(__LINE__, size(help, kind = IK) < nintmax, SK_"@getQuadErr(): The condition `size(help) < size(sindex)` must hold. size(help), size(sindex) = "//getStr([size(help, kind = IK), nintmax]))
848 576 : CHECK_ASSERTION(__LINE__, all(help > lb), SK_"@getQuadErr(): The condition `all(lb < help)` must hold for the input arguments `help` and `lb`. lb, help = "//getStr([lb, help]))
849 576 : CHECK_ASSERTION(__LINE__, all(help < ub), SK_"@getQuadErr(): The condition `all(help < ub)` must hold for the input arguments `help` and `ub`. help, ub = "//getStr([help, ub]))
850 144 : CHECK_ASSERTION(__LINE__, size(sinfo, 2, IK) == nintmax, SK_"@getQuadErr(): The condition `size(sinfo,2) == size(sindex)` must hold. size(sinfo,2), size(sindex) = "//getStr([size(sinfo, 2, IK), nintmax]))
851 48 : CHECK_ASSERTION(__LINE__, isAscending(help), SK_"@getQuadErr(): The input argument `help` must be in ascending order. help = "//getStr(help))
852 : ! Define the first set of interval edges.
853 48 : sinfo(1,1) = lb
854 208 : do i = 1_IK, lenBreak
855 160 : sinfo(2, i) = help(i)
856 208 : sinfo(1, i + 1) = sinfo(2, i)
857 : end do
858 48 : sinfo(2,i) = ub
859 : ! Compute the first approximation to the integral.
860 48 : abserr = 0._RKC
861 0 : integral = 0._RKC
862 48 : intAbsFunc1 = 0._RKC
863 256 : do i = 1_IK, lenBreakPlusOne ! loop over the first set of intervals.
864 208 : result1 = getQuadGK(getFunc, sinfo(1,i), sinfo(2,i), QRULE_ARG, abserr1, intAbsFunc2, smoothness1) ! Compute the first approximation to the integral.
865 208 : neval = neval + NEVAL0
866 208 : refinementNeeded(i) = abserr1 /= 0._RKC .and. abserr1 == smoothness1
867 208 : intAbsFunc1 = intAbsFunc1 + intAbsFunc2
868 208 : integral = integral + result1
869 208 : abserr = abserr + abserr1
870 208 : sinfo(3,i) = result1
871 208 : sinfo(4,i) = abserr1
872 208 : level(i) = 0_IK
873 256 : sindex(i) = i
874 : end do
875 : ! Set the overall error estimate.
876 48 : errsum = 0._RKC
877 256 : do i = 1_IK, lenBreakPlusOne
878 208 : if (refinementNeeded(i)) sinfo(4,i) = abserr
879 256 : errsum = errsum + sinfo(4,i)
880 : end do
881 48 : nint = lenBreakPlusOne
882 : ! Sort the error estimate indices in descending order of error estimates.
883 208 : do i = 1_IK, lenBreak
884 160 : jlow = i + 1_IK
885 160 : ind1 = sindex(i)
886 552 : do j = jlow, lenBreakPlusOne
887 392 : ind2 = sindex(j)
888 552 : if (sinfo(4,ind1) <= sinfo(4,ind2)) then
889 : ind1 = ind2
890 : k = j
891 : end if
892 : end do
893 208 : if (ind1 /= sindex(i)) then
894 99 : sindex(k) = sindex(i)
895 99 : sindex(i) = ind1
896 : end if
897 : end do
898 : #elif QAGD_ENABLED || QAGS_ENABLED || QAWC_ENABLED
899 : #if QAWC_ENABLED
900 16 : neval = 0_IK
901 16 : call setQuadQAWC(integral, lb, ub, abserr, intAbsFunc1, smoothness1, neval)
902 : #elif QAGD_ENABLED || QAGS_ENABLED
903 91430 : integral = getQuadGK(getFunc, lb, ub, QRULE_ARG, abserr, intAbsFunc1, smoothness1)
904 91430 : neval = NEVAL0
905 : #endif
906 91446 : sinfo(1,1) = lb
907 91446 : sinfo(2,1) = ub
908 91446 : sinfo(3,1) = integral
909 91446 : sinfo(4,1) = abserr
910 91446 : sindex(1) = 1_IK
911 91446 : nint = 1_IK
912 : #else
913 : #error "Unrecognized interface."
914 : #endif
915 : ! Test the accuracy of the zeroth order integration.
916 : err = 0_IK
917 91494 : errbnd = max(abstol, reltol * abs(integral))
918 91494 : if (nint == nintmax) then; err = 1_IK; return; end if
919 91494 : if (errbnd < abserr .and. abserr <= FACTOR * 50._RKC * EPS_RKC * intAbsFunc1) then; err = 2_IK; return; end if
920 : #if QAGP_ENABLED
921 48 : if (abserr <= errbnd) return
922 : #elif QAGD_ENABLED || QAGS_ENABLED || QAWC_ENABLED
923 91446 : if (abserr == 0._RKC) return
924 91446 : if (abserr <= errbnd .and. abserr /= smoothness1) return
925 : #else
926 : #error "Unrecognized interface."
927 : #endif
928 : !#if QAWC_ENABLED
929 : ! if (abserr <= min(0.01_RKC * abs(integral), errbnd)) return
930 : !#endif
931 : ! Refine the integration.
932 : #if QAGP_ENABLED
933 48 : erlarg = errsum
934 48 : ertest = errbnd
935 : maximumLevel = 1_IK
936 : #elif QAGD_ENABLED || QAGS_ENABLED || QAWC_ENABLED
937 6915 : errsum = abserr
938 : #else
939 : #error "Unrecognized interface."
940 : #endif
941 6963 : result = integral
942 40026 : maxErrLoc = sindex(1)
943 40026 : maxErrVal = sinfo(4,maxErrLoc)
944 : iroff1 = 0_IK
945 : iroff2 = 0_IK
946 40026 : nrmax = 1_IK
947 : #if QAGS_ENABLED || QAGP_ENABLED
948 : ierro = 0_IK
949 : ktmin = 0_IK
950 7353 : nextrap = 0_IK
951 : iroffnew = 0_IK
952 7353 : abserr = HUGE_RKC
953 7353 : EpsTable(1) = integral
954 7353 : lenEpsTable = 1_IK + EPS_TABLE_START_OFFSET ! remark: The second element of `EpsTable` for QAGS will be set later on.
955 : extrapAllowed = .true._LK
956 : extrapolating = .false._LK
957 : ! \warning \bug The quadpack.F90 version of John Burkardt is different from the original quadpack, where 50 is replaced with 0.5 in QAGP.
958 7353 : ksgn = -1_IK; if (abs(integral) >= (1._RKC - 50._RKC * EPS_RKC) * intAbsFunc1) ksgn = 1_IK
959 : #elif !(QAGD_ENABLED || QAWC_ENABLED)
960 : #error "Unrecognized interface."
961 : #endif
962 : !print *, "abserr, errsum", abserr, errsum
963 : ! Adaptively integrate the integrand.
964 171138 : loopAdaptation: do nint = nint + 1_IK, nintmax
965 : ! Bisect the subinterval that has the largest error estimate.
966 171138 : lb1 = sinfo(1,maxErrLoc)
967 171138 : ub2 = sinfo(2,maxErrLoc)
968 171138 : ub1 = 0.5_RKC * (lb1 + ub2)
969 : #if QAWC_ENABLED
970 112 : if (lb1 < cs .and. cs <= ub1) then
971 40 : ub1 = 0.5_RKC * (cs + ub2)
972 72 : elseif (ub1 < cs .and. cs < ub2) then
973 24 : ub1 = 0.5_RKC * (lb1 + cs)
974 : end if
975 : #endif
976 171138 : lb2 = ub1
977 : #if QAGS_ENABLED || QAGP_ENABLED
978 63947 : erlast = maxErrVal
979 : #if QAGP_ENABLED
980 3930 : currentLevel = level(maxErrLoc) + 1_IK
981 : #endif
982 : #elif !(QAGD_ENABLED || QAWC_ENABLED)
983 : #error "Unrecognized interface."
984 : #endif
985 : #if QAWC_ENABLED
986 112 : call setQuadQAWC(result1, lb1, ub1, abserr1, intAbsFunc2, smoothness1, neval)
987 112 : call setQuadQAWC(result2, lb2, ub2, abserr2, intAbsFunc2, smoothness2, neval)
988 : !print *, abserr1, abserr2, lb1 < cs .and. cs <= ub1
989 : #elif QAGD_ENABLED || QAGS_ENABLED || QAGP_ENABLED
990 171026 : result1 = getQuadGK(getFunc, lb1, ub1, QRULE_ARG, abserr1, intAbsFunc2, smoothness1)
991 171026 : result2 = getQuadGK(getFunc, lb2, ub2, QRULE_ARG, abserr2, intAbsFunc2, smoothness2)
992 171026 : neval = neval + 2_IK * NEVAL0
993 : #else
994 : #error "Unrecognized interface."
995 : #endif
996 : ! Improve the previous approximations to integral and error and verify the accuracy.
997 171138 : result12 = result1 + result2
998 171138 : abserr12 = abserr1 + abserr2
999 171138 : errsum = errsum + abserr12 - maxErrVal
1000 171138 : result = result + result12 - sinfo(3,maxErrLoc)
1001 171138 : if (smoothness1 /= abserr1 .and. smoothness2 /= abserr2) then
1002 124621 : if (nint > 10_IK .and. abserr12 > maxErrVal) iroff2 = iroff2 + 1_IK
1003 124621 : if (abs(sinfo(3,maxErrLoc) - result12) <= 1.e-05_RKC * abs(result12) .and. abserr12 >= 0.99_RKC * maxErrVal) then
1004 : #if QAGD_ENABLED || QAWC_ENABLED
1005 37 : iroff1 = iroff1 + 1_IK
1006 : #elif QAGS_ENABLED || QAGP_ENABLED
1007 30 : if (extrapolating) then
1008 16 : iroffnew = iroffnew + 1_IK
1009 : else
1010 14 : iroff1 = iroff1 + 1_IK
1011 : end if
1012 : #else
1013 : #error "Unrecognized interface."
1014 : #endif
1015 : end if
1016 : end if
1017 : #if QAGP_ENABLED
1018 3930 : level(maxErrLoc) = currentLevel
1019 3930 : level(nint) = currentLevel
1020 : #endif
1021 : ! Append the newly-created intervals to the list.
1022 171138 : CHECK_ASSERTION(__LINE__, nint /= maxErrLoc, SK_"@getQuadErr(): The condition `nint /= maxErrLoc` happened. This in an internal error. Please report this issue to the developers.")
1023 171138 : if (abserr2 > abserr1) then ! \todo is this order necessary?
1024 24625 : sinfo(1,maxErrLoc) = lb2
1025 24625 : sinfo(1,nint) = lb1
1026 24625 : sinfo(2,nint) = ub1
1027 24625 : sinfo(3,maxErrLoc) = result2
1028 24625 : sinfo(3,nint) = result1
1029 24625 : sinfo(4,maxErrLoc) = abserr2
1030 24625 : sinfo(4,nint) = abserr1
1031 : else
1032 146513 : sinfo(1,nint) = lb2
1033 146513 : sinfo(2,maxErrLoc) = ub1
1034 146513 : sinfo(2,nint) = ub2
1035 146513 : sinfo(3,maxErrLoc) = result1
1036 146513 : sinfo(3,nint) = result2
1037 146513 : sinfo(4,maxErrLoc) = abserr1
1038 146513 : sinfo(4,nint) = abserr2
1039 : end if
1040 : ! Ensure the descending ordering of the list of error estimates and select the subinterval with the largest error estimate (to be bisected next).
1041 171138 : call setErrSorted(nintmax, sinfo(4,1:nint), sindex(1:nint), nrmax = nrmax, maxErrLoc = maxErrLoc, maxErrVal = maxErrVal)
1042 : ! Check for convergence/divergence.
1043 171138 : errbnd = max(abstol, reltol * abs(result))
1044 171138 : if (errsum <= errbnd) then ! convergence occurred.
1045 208304 : FINALIZE_INTEGRATION
1046 : end if
1047 : #if QAGD_ENABLED || QAWC_ENABLED
1048 74577 : if (iroff1 >= 6_IK .or. iroff2 >= 20_IK) err = 2_IK ! test for roundoff error and eventually set error flag.
1049 74577 : if (nint == nintmax) err = 1_IK ! set error flag in the case that the number of subintervals equals nintmax.
1050 : ! \warning \bug QuadPack version of John Burkardt, replaces 1000 in the original code with 10000 in the condition below.
1051 74577 : if (max(abs(lb1), abs(ub2)) <= (1._RKC + NEVAL0 * 100._RKC * EPS_RKC) * (abs(lb2) + 1000._RKC * TINY_RKC)) err = 3_IK !> \warning \todo The impact of NEVAL0 in this condition may need further investigation for high-order > 61 rules.
1052 : #elif QAGS_ENABLED || QAGP_ENABLED
1053 57743 : if (iroff1 + iroffnew >= 10_IK .or. iroff2 >= 20) err = 2_IK ! roundoff error.
1054 57743 : if (iroffnew >= 5_IK) ierro = 3_IK
1055 : !if (iroffnew >= 5_IK) error stop __LINE__
1056 57743 : if (nint == nintmax) err = 1_IK
1057 : ! \warning \bug QuadPack version of John Burkardt, replaces 100 in the original code with 1000 in the condition below.
1058 57743 : if (max(abs(lb1), abs(ub2)) <= (1._RKC + 100._RKC * EPS_RKC) * (abs(lb2) + 1000._RKC * TINY_RKC)) err = 4_IK ! bad integrand behavior.
1059 : #else
1060 : #error "Unrecognized interface."
1061 : #endif
1062 132376 : if (err /= 0_IK) exit loopAdaptation
1063 : #if QAGS_ENABLED
1064 : ! Extrapolate, if allowed.
1065 53810 : if (nint == 2_IK) then
1066 6965 : small = abs(ub - lb) * 0.375_RKC
1067 3256 : erlarg = errsum
1068 3256 : ertest = errbnd
1069 6965 : EpsTable(2) = result
1070 6965 : cycle loopAdaptation
1071 : end if
1072 : #endif
1073 : #if QAGS_ENABLED || QAGP_ENABLED
1074 51924 : if (extrapAllowed) then
1075 50775 : erlarg = erlarg - erlast
1076 : if ( & ! LCOV_EXCL_LINE
1077 : #if QAGS_ENABLED
1078 : small < abs(ub1 - lb1) & ! LCOV_EXCL_LINE
1079 : #elif QAGP_ENABLED
1080 : currentLevel < maximumLevel & ! LCOV_EXCL_LINE
1081 : #else
1082 : #error "Unrecognized intefrace."
1083 : #endif
1084 362 : ) erlarg = erlarg + abserr12
1085 50775 : if (.not. extrapolating) then ! test whether the interval to be bisected next is the smallest interval.
1086 30713 : if (CYCLE_NEEDED) cycle loopAdaptation
1087 : extrapolating = .true._LK
1088 24098 : nrmax = 2_IK
1089 : end if
1090 44160 : if (ierro /= 3_IK .and. erlarg > ertest) then
1091 : ! The smallest interval has the largest error. Before bisecting decrease the
1092 : ! sum of the errors over the larger intervals (erlarg) and perform extrapolation.
1093 : jupbnd = nint
1094 20064 : if (nint > 2_IK + nintmax / 2_IK) jupbnd = nintmax + 3_IK - nint
1095 28907 : do k = nrmax, jupbnd
1096 28907 : maxErrLoc = sindex(nrmax)
1097 28907 : maxErrVal = sinfo(4,maxErrLoc)
1098 28907 : if (CYCLE_NEEDED) cycle loopAdaptation
1099 8843 : nrmax = nrmax + 1_IK
1100 : end do
1101 : end if
1102 : ! perform extrapolation.
1103 24096 : lenEpsTable = lenEpsTable + 1_IK
1104 24096 : EpsTable(lenEpsTable) = result
1105 : #if QAGP_ENABLED
1106 634 : if (lenEpsTable > 2_IK) then
1107 : #elif !QAGS_ENABLED
1108 : #error "Unrecognized intefrace."
1109 : #endif
1110 : !print *, "lenEpsTable, nextrap = ", lenEpsTable, nextrap
1111 24048 : call setSeqLimEps(lenEpsTable, nextrap, EpsTable, seqlims, seqlim, absErrEps)
1112 24048 : ktmin = ktmin + 1_IK
1113 : !print *, "end0", "ktmin", ktmin
1114 : !print *, "end0", "errsum", errsum
1115 : !print *, "end0", "abserr", abserr
1116 24048 : if (ktmin > 5_IK .and. abserr < 1.e-03_RKC * errsum) err = 5_IK
1117 : !print *, "end1", MAXLEN_EPSTAB, lenEpsTable
1118 24048 : if (absErrEps < abserr) then
1119 : ktmin = 0_IK
1120 10368 : abserr = absErrEps
1121 10368 : integral = seqlim
1122 7607 : correc = erlarg
1123 10368 : ertest = max(abstol, reltol * abs(seqlim))
1124 10368 : if (abserr <= ertest) exit loopAdaptation
1125 : end if
1126 : ! prepare bisection of the smallest interval.
1127 22908 : if (lenEpsTable == 1_IK) extrapAllowed = .false._LK
1128 22908 : if (err == 5_IK) exit loopAdaptation
1129 : #if QAGP_ENABLED
1130 : end if
1131 586 : maximumLevel = maximumLevel + 1_IK
1132 : #elif QAGS_ENABLED
1133 22364 : small = small * 0.5_RKC
1134 : #else
1135 : #error "Unrecognized intefrace."
1136 : #endif
1137 22950 : maxErrLoc = sindex(1)
1138 22950 : maxErrVal = sinfo(4,maxErrLoc)
1139 : extrapolating = .false._LK
1140 13846 : erlarg = errsum
1141 22950 : nrmax = 1_IK
1142 : end if
1143 : #endif
1144 :
1145 : end do loopAdaptation
1146 3624 : CHECK_ASSERTION(__LINE__, nint <= nintmax, \
1147 : SK_"@getQuadErr(): Internal library error: The condition `nint <= nintmax` must hold. nint, nintmax = ["//\
1148 : getStr([nint, nintmax])//SK_"]. Please report this error to the ParaMonte library developers.")
1149 :
1150 : #if QAGS_ENABLED || QAGP_ENABLED
1151 1149 : if (abserr /= HUGE_RKC) then
1152 1149 : if (err + ierro /= 0_IK) then
1153 9 : if (ierro == 3_IK) abserr = abserr + correc
1154 9 : if (err == 0_IK) err = 3_IK
1155 : !if (err == 3) error stop __LINE__
1156 9 : if (integral == 0._RKC .or. result == 0._RKC) then
1157 0 : if (abserr > errsum) then
1158 0 : FINALIZE_INTEGRATION
1159 : end if
1160 0 : if (result == 0._RKC) then
1161 0 : if (err > 2_IK) err = err - 1_IK
1162 0 : return
1163 : end if
1164 9 : elseif (abserr / abs(integral) > errsum / abs(result)) then
1165 2001 : FINALIZE_INTEGRATION
1166 : end if
1167 : end if
1168 1148 : if (ksgn /= -1_IK .or. max(abs(integral), abs(result)) > intAbsFunc1 * 0.01_RKC) then ! test on divergence.
1169 1140 : if (0.01_RKC > integral / result .or. integral / result > 100._RKC .or. errsum > abs(result)) err = 6_IK
1170 : end if
1171 1148 : if (err > 2_IK) err = err - 1_IK
1172 1148 : return
1173 : end if
1174 : #elif !(QAGD_ENABLED || QAWC_ENABLED)
1175 : #error "Unrecognized interface."
1176 : #endif
1177 14168 : FINALIZE_INTEGRATION
1178 :
1179 : contains
1180 :
1181 : #if QAWC_ENABLED
1182 : ! Approximate the integral of the input function using,
1183 : ! <ul>
1184 : ! <li> a generalized clenshaw-curtis method if `cs` lies within ten percent of the integration limits, otherwise,
1185 : ! <li> a user-specified Gauss-Kronrod quadrature rule.
1186 : ! </ul>
1187 : ! This routine uses `getFunc`, `getFuncWC`, `NEVAL0`, and `cs` objects from the parent routine.
1188 : ! However, it does not have any side effects (does not change any global variables).
1189 240 : subroutine setQuadQAWC(quadQAWC, lb, ub, abserr, intAbsFunc, smoothness, neval)
1190 : real(RKC) , intent(in) :: lb, ub ! the integration lower/upper bound.
1191 : real(RKC) , intent(out) :: abserr ! the abslue value of the integral error estimate, which equals or exceeds `abs(i - quadQAWC)`.
1192 : real(RKC) , intent(out) :: intAbsFunc ! integral of the absolute of the function.
1193 : real(RKC) , intent(out) :: smoothness ! An integrand smoothness measure.
1194 : integer(IK) , intent(inout) :: neval ! number of function evaluations.
1195 : real(RKC) , intent(out) :: quadQAWC ! approximation to the integral.
1196 : integer(IK) :: i, j
1197 : real(RKC) , parameter :: NODE(11) = [(cos(i * acos(-1._RKC) / 24._RKC), i = 1_IK, 11_IK)]
1198 : real(RKC) :: ak22, amom0, amom1, amom2, u, cstrans
1199 : real(RKC) :: cheb12(13) ! Chebyshev series expansion of degree 12.
1200 : real(RKC) :: cheb24(25) ! Chebyshev series expansion of degree 24.
1201 : real(RKC) :: func(25) ! Function values at `cos(k*pi/24)`, `k = 0, ..., 24`.
1202 : real(RKC) :: result12, result24 ! approximations to the integral Via Clenshaw-Curtis of orders 13 and 25.
1203 : real(RKC) :: midInterval, halfInterval
1204 240 : cstrans = (2._RKC * cs - ub - lb) / (ub - lb) ! Compute the position of cs.
1205 : !print *, cstrans
1206 240 : if (abs(cstrans) < 1.1_RKC) then ! Use the generalized Clenshaw-Curtis method.
1207 80 : smoothness = huge(smoothness)
1208 80 : halfInterval = 0.5_RKC * (ub - lb)
1209 80 : midInterval = 0.5_RKC * (ub + lb)
1210 80 : func(1) = 0.5_RKC * getFunc(halfInterval + midInterval)
1211 80 : func(13) = getFunc(midInterval)
1212 80 : func(25) = 0.5_RKC * getFunc(midInterval - halfInterval)
1213 960 : do i = 2_IK, 12_IK
1214 880 : u = halfInterval * NODE(i - 1_IK)
1215 880 : func(i) = getFunc(u + midInterval)
1216 960 : func(26_IK - i) = getFunc(midInterval - u)
1217 : end do
1218 80 : call setChebExpan(func, cheb12, cheb24) ! compute the Chebyshev series expansion.
1219 : ! the modified Chebyshev moments are computed by forward recursion, using amom0 and amom1 as starting values.
1220 80 : amom0 = log(abs((1._RKC - cstrans) / (1._RKC + cstrans)))
1221 80 : amom1 = 2._RKC + cstrans * amom0
1222 80 : result12 = cheb12(1) * amom0 + cheb12(2) * amom1
1223 80 : result24 = cheb24(1) * amom0 + cheb24(2) * amom1
1224 80 : intAbsFunc = abs(cheb24(1))
1225 : j = 1_IK
1226 960 : do i = 3_IK, 13_IK
1227 880 : amom2 = 2._RKC * cstrans * amom1 - amom0
1228 880 : ak22 = real((i - 2_IK) * (i - 2_IK), RKC)
1229 880 : if (i == (i / 2_IK) * 2_IK) amom2 = amom2 - 4._RKC / (ak22 - 1._RKC)
1230 880 : result12 = result12 + cheb12(i) * amom2
1231 880 : result24 = result24 + cheb24(i) * amom2
1232 880 : amom0 = amom1
1233 880 : amom1 = amom2
1234 880 : j = j + 2_IK
1235 960 : intAbsFunc = intAbsFunc + abs(cheb24(j - 1)) + abs(cheb24(j))
1236 : end do
1237 1040 : do i = 14_IK, 25_IK
1238 960 : amom2 = 2._RKC * cstrans * amom1 - amom0
1239 960 : ak22 = real((i - 2_IK) * (i - 2_IK), RKC)
1240 960 : if (i == (i / 2_IK) * 2_IK) amom2 = amom2 - 4._RKC / (ak22 - 1._RKC)
1241 960 : result24 = result24 + cheb24(i) * amom2
1242 960 : amom0 = amom1
1243 1040 : amom1 = amom2
1244 : end do
1245 80 : abserr = abs(result24 - result12)
1246 80 : quadQAWC = result24
1247 80 : neval = neval + 25_IK
1248 : else ! Use the default Gauss-Kronrod quadrature.
1249 160 : quadQAWC = getQuadGK(getFuncWC, lb, ub, QRULE_ARG, abserr, intAbsFunc, smoothness)
1250 160 : neval = neval + NEVAL0
1251 : end if
1252 240 : end subroutine
1253 :
1254 : ! Cauchy-Weighted function.
1255 :
1256 4184 : function getFuncWC(x) result(funcWC)
1257 : real(RKC) , intent(in) :: x
1258 : real(RKC) :: funcWC
1259 4184 : funcWC = getFunc(x) / (x - cs)
1260 4184 : end function
1261 :
1262 : #elif QAWFS_ENABLED || QAWFC_ENABLED
1263 :
1264 : subroutine getQuadWF(quadWF, lb, ub, currentLevel, ksave, abserr, intAbsFunc, smoothness, chebmom, momcom, neval)
1265 : real(RKC) , intent(in) :: lb, ub !< The current integration limits.
1266 : integer(IK) , intent(in) :: ksave !< The key which is one when the moments for the current interval have been computed.
1267 : integer(IK) , intent(in) :: currentLevel !< The length of interval `(lb,ub)` is equal to the length of the original integration interval divided by `2**currentLevel`.
1268 : !! We assume that the routine is used in an adaptive integration process, otherwise set `currentLevel = 0`.
1269 : !! `currentLevel` must be zero at the first call.
1270 : real(RKC) , intent(out) :: abserr !< estimate of the modulus of the absolute error, which should equal or exceed `abs(i-quadWF)`.
1271 : real(RKC) , intent(out) :: intAbsFunc !< Approximation to the integral of the absolute of the function `abs(f(x))`.
1272 : real(RKC) , intent(out) :: smoothness !< Approximation to the integral of `abs(f - i / (ub - lb))`.
1273 : real(RKC) , intent(inout) , contiguous :: chebmom(:, 25) !< Array of shape `(maxLenChebMom,25)` containing the modified Chebyshev moments for the first `momcom` interval lengths.
1274 : integer(IK) , intent(inout) :: momcom !< For each interval length we need to compute the Chebyshev moments.
1275 : !! `momcom` counts the number of intervals for which these moments have already been computed.
1276 : !! If `currentLevel < momcom` or `ksave = 1`, the Chebyshev moments for the interval `(lb, ub)` have already been computed.
1277 : !! Otherwise we compute them and we increase `momcom`.
1278 : integer(IK) , intent(out) :: neval !< Number of integrand evaluations.
1279 : real(RKC) :: quadWF !< approximation to the integral i
1280 :
1281 : real(RKC) , parameter :: NODE(11) = [(cos(i * acos(-1._RKC) / 24._RKC), i = 1_IK, 11_IK)]
1282 : integer(IK) :: i, iers, isym, j, k, m, noequ, noeq1
1283 : integer(IK) :: maxLenChebMom !< The upper bound on the number of Chebyshev moments which can be stored,
1284 : !< i.e. for the intervals of lengths `abs(bb-aa)*2**(-l)`, `l = 0,1,2, ..., maxLenChebMom-2`.
1285 : real(RKC) :: ac, an, ansq, invansq, as, asap, ass, conc, cons, cospar, d(25), d1(25), d2(25), estc, ests, parint, par2, par22, p2, p3, p4, sinpar, Vector(28)
1286 : real(RKC) :: halfInterval
1287 : real(RKC) :: midInterval
1288 : real(RKC) :: cheb12(13) !< Chebyshev series expansion of degree 12.
1289 : real(RKC) :: cheb24(25) !< Chebyshev series expansion of degree 24.
1290 : real(RKC) :: func(25) !< Function values at `(ub - lb) * 0.5 * cos(k * pi / 12) + (ub + lb) * 0.5`, `k = 0, ..., 24`.
1291 : real(RKC) :: result12cos !< approximation to the integral of `cos(0.5 * (ub - lb) * omega * NODE) * getFunc(0.5 * (ub - lb) * NODE + 0.5 * (ub + lb))` over `(-1, +1)`, using the Chebyshev series expansion of degree 12.
1292 : real(RKC) :: result24cos !< approximation to the integral of `cos(0.5 * (ub - lb) * omega * NODE) * getFunc(0.5 * (ub - lb) * NODE + 0.5 * (ub + lb))` over `(-1, +1)`, using the Chebyshev series expansion of degree 12.
1293 : real(RKC) :: result12sin !< approximation to the integral of `sin(0.5 * (ub - lb) * omega * NODE) * getFunc(0.5 * (ub - lb) * NODE + 0.5 * (ub + lb))` over `(-1, +1)`, using the Chebyshev series expansion of degree 12.
1294 : real(RKC) :: result24sin !< approximation to the integral of `sin(0.5 * (ub - lb) * omega * NODE) * getFunc(0.5 * (ub - lb) * NODE + 0.5 * (ub + lb))` over `(-1, +1)`, using the Chebyshev series expansion of degree 12.
1295 :
1296 : midInterval = 0.5_RKC * (ub + lb)
1297 : halfInterval = 0.5_RKC * (ub - lb)
1298 : parint = Weight%omega * halfInterval
1299 : maxLenChebMom = size(chebmom, 1_IK, IK)
1300 :
1301 : ! Compute the integral using the 15-point gauss-kronrod formula if the value of the parameter in the integrand is small.
1302 :
1303 : if (abs(parint) > 2._RKC) then
1304 :
1305 : ! Compute the integral using the generalized clenshaw-curtis method.
1306 :
1307 : conc = halfInterval * cos(midInterval * Weight%omega)
1308 : cons = halfInterval * sin(midInterval * Weight%omega)
1309 : smoothness = huge(smoothness)
1310 : neval = 25_IK
1311 :
1312 : ! Check whether the Chebyshev moments for this interval have already been computed.
1313 :
1314 : if (momcom <= currentLevel .and. ksave /= 1_IK) then ! Compute a new set of Chebyshev moments.
1315 :
1316 : m = momcom + 1_IK
1317 : par2 = parint*parint
1318 : par22 = par2 + 2._RKC
1319 : sinpar = sin(parint)
1320 : cospar = cos(parint)
1321 :
1322 : ! compute the Chebyshev moments with respect to cosine.
1323 :
1324 : Vector(1) = 2._RKC*sinpar / parint
1325 : Vector(2) = (8._RKC*cospar + (par2 + par2 - 8._RKC) * sinpar / parint) / par2
1326 : Vector(3) = (32._RKC*(par2 - 12._RKC)*cospar + (2._RKC * ((par2 - 80._RKC) * par2 + 192._RKC) * sinpar) / parint) / (par2*par2)
1327 : ac = 8._RKC * cospar
1328 : as = 24._RKC * parint * sinpar
1329 : if (abs(parint) > 24._RKC) then ! Compute the Chebyshev moments by means of forward recursion.
1330 : an = 4._RKC
1331 : do i = 4_IK, 13_IK
1332 : ansq = an * an
1333 : Vector(i) = ((ansq - 4._RKC) * (2._RKC * (par22 - ansq - ansq) * Vector(i - 1) - ac) + as - par2 * (an + 1._RKC) * (an + 2._RKC) * Vector(i - 2)) / (par2 * (an - 1._RKC) * (an - 2._RKC))
1334 : an = an + 2._RKC
1335 : end do
1336 : else ! Compute the Chebyshev moments as the solutions of a boundary value problem with 1 initial value (Vector(3)) and 1 end value (computed using an asymptotic formula).
1337 : noequ = 25_IK
1338 : noeq1 = noequ - 1_IK
1339 : an = 6._RKC
1340 : do k = 1_IK, noeq1
1341 : ansq = an * an
1342 : d(k) = -2._RKC * (ansq - 4._RKC) * (par22 - ansq - ansq)
1343 : d2(k) = (an - 1._RKC) * (an - 2._RKC) * par2
1344 : d1(k + 1) = (an + 3._RKC) * (an + 4._RKC) * par2
1345 : Vector(k + 3) = as - (ansq - 4._RKC) * ac
1346 : an = an + 2._RKC
1347 : end do
1348 : ansq = an * an
1349 : invansq = 1._RKC / ansq
1350 : d(noequ) = -2._RKC * (ansq - 4._RKC) * (par22 - ansq - ansq)
1351 : Vector(noequ + 3) = as - (ansq - 4._RKC) * ac
1352 : Vector(4) = Vector(4) - 56._RKC * par2 * Vector(3)
1353 : ass = parint * sinpar
1354 : asap = (((((210._RKC * par2 - 1._RKC) * cospar - (105._RKC * par2 - 63._RKC) * ass) * invansq - (1._RKC - 15._RKC * par2) * cospar + 15._RKC * ass) * invansq - cospar + 3._RKC * ass) * invansq - cospar) * invansq
1355 : Vector(noequ + 3) = Vector(noequ + 3) - 2._RKC * asap * par2 * (an - 1._RKC) * (an - 2._RKC)
1356 : ! Solve the tridiagonal system by means of Gaussian elimination with partial pivoting.
1357 : call dgtsl(noequ, d1, d, d2, Vector(4), iers)
1358 : end if
1359 : do j = 1_IK, 13_IK
1360 : chebmom(m, 2_IK * j - 1_IK) = Vector(j)
1361 : end do
1362 :
1363 : ! Compute the Chebyshev moments with respect to sine.
1364 :
1365 : Vector(1) = 2._RKC * (sinpar - parint * cospar) / par2
1366 : Vector(2) = (18._RKC - 48._RKC / par2) * sinpar / par2 + (-2._RKC + 48._RKC / par2) * cospar / parint
1367 : ac = -24._RKC * parint * cospar
1368 : as = -8._RKC * sinpar
1369 : if (abs(parint) > 24._RKC) then ! Compute the Chebyshev moments by means of forward recursion.
1370 : an = 3._RKC
1371 : do i = 3_IK, 12_IK
1372 : ansq = an * an
1373 : Vector(i) = ((ansq - 4._RKC) * (2._RKC * (par22 - ansq - ansq) * Vector(i - 1) + as) + ac - par2 * (an + 1._RKC) * (an + 2._RKC) * Vector(i - 2)) / (par2 * (an - 1._RKC) * (an - 2._RKC))
1374 : an = an + 2._RKC
1375 : end do
1376 : else ! Compute the Chebyshev moments as the solutions of a boundary value problem with 1 initial value (Vector(2)) and 1 end value (computed using an asymptotic formula).
1377 : an = 5._RKC
1378 : do k = 1_IK, noeq1
1379 : ansq = an * an
1380 : d(k) = -2._RKC * (ansq - 4._RKC) * (par22 - ansq - ansq)
1381 : d2(k) = (an - 1._RKC) * (an - 2._RKC) * par2
1382 : d1(k + 1) = (an + 3._RKC) * (an + 4._RKC) * par2
1383 : Vector(k + 2) = ac + (ansq - 4._RKC) * as
1384 : an = an + 2._RKC
1385 : end do
1386 : ansq = an * an
1387 : d(noequ) = -2._RKC * (ansq - 4._RKC) * (par22 - ansq - ansq)
1388 : Vector(noequ + 2) = ac + (ansq - 4._RKC) * as
1389 : Vector(3) = Vector(3) - 42._RKC * par2 * Vector(2)
1390 : ass = parint * cospar
1391 : invansq = 1._RKC / ansq
1392 : asap = (((((105._RKC * par2 - 63._RKC) * ass + (210._RKC * par2 - 1._RKC) * sinpar) * invansq + (15._RKC * par2 - 1._RKC) * sinpar - 15._RKC * ass) * invansq - 3._RKC * ass - sinpar) * invansq - sinpar) * invansq
1393 : Vector(noequ + 2) = Vector(noequ + 2) - 2._RKC * asap * par2 * (an - 1._RKC) * (an - 2._RKC)
1394 : ! Solve the tridiagonal system by means of Gaussian elimination with partial pivoting.
1395 : call dgtsl(noequ, d1, d, d2, Vector(3), iers)
1396 : end if
1397 : do j = 1_IK, 12_IK
1398 : chebmom(m, 2 * j) = Vector(j)
1399 : end do
1400 : end if
1401 : if (currentLevel < momcom) m = currentLevel + 1_IK
1402 : if (momcom < maxLenChebMom - 1_IK .and. momcom <= currentLevel) momcom = momcom + 1_IK
1403 :
1404 : ! Compute the coefficients of the Chebyshev expansions of degrees 12 and 24 of the function `f(x)`.
1405 :
1406 : func(1) = 0.5_RKC * f(midInterval + halfInterval)
1407 : func(13) = f(midInterval)
1408 : func(25) = 0.5_RKC * f(midInterval - halfInterval)
1409 : do i = 2_IK, 12_IK
1410 : isym = 26_IK - i
1411 : func(i) = f(halfInterval * NODE(i - 1) + midInterval)
1412 : func(isym) = f(midInterval - halfInterval * NODE(i - 1))
1413 : end do
1414 : call dqcheb(x, func, cheb12, cheb24)
1415 :
1416 : ! Compute the integral and error estimates.
1417 :
1418 : result12cos = cheb12(13) * chebmom(m, 13)
1419 : result12sin = 0._RKC
1420 : k = 11_IK
1421 : do j = 1_IK, 6_IK
1422 : result12cos = result12cos + cheb12(k) * chebmom(m, k)
1423 : result12sin = result12sin + cheb12(k + 1) * chebmom(m, k + 1)
1424 : k = k - 2_IK
1425 : end do
1426 : result24cos = cheb24(25) * chebmom(m, 25)
1427 : result24sin = 0._RKC
1428 : intAbsFunc = abs(cheb24(25))
1429 : k = 23_IK
1430 : do j = 1_IK, 12_IK
1431 : result24cos = result24cos + cheb24(k) * chebmom(m, k)
1432 : result24sin = result24sin + cheb24(k + 1) * chebmom(m, k + 1)
1433 : intAbsFunc = intAbsFunc + abs(cheb24(k)) + abs(cheb24(k + 1))
1434 : k = k - 2_IK
1435 : end do
1436 : estc = abs(result24cos - result12cos)
1437 : ests = abs(result24sin - result12sin)
1438 : intAbsFunc = intAbsFunc * abs(halfInterval)
1439 : if (Integr == 2_IK) then
1440 : quadWF = conc * result24sin + cons * result24cos
1441 : abserr = abs(conc * ests) + abs(cons * estc)
1442 : else
1443 : quadWF = conc * result24cos - cons * result24sin
1444 : abserr = abs(conc * estc) + abs(cons * ests)
1445 : end if
1446 : else
1447 : call dqk15w(f, dqwgtf, Weight%omega, p2, p3, p4, Integr, lb, ub, quadWF, abserr, intAbsFunc, smoothness)
1448 : neval = 15_IK
1449 : end if
1450 : end function
1451 :
1452 : function getFuncWeightedSin(x) result(funcWeightedFour)
1453 : real(RKC) , intent(in) :: x
1454 : real(RKC) :: funcWeightedFour
1455 : funcWeightedFour = sin(Weight%omega * x) * getFunc(x)
1456 : end function
1457 :
1458 : function getFuncWeightedSin(x) result(funcWeightedFour)
1459 : real(RKC) , intent(in) :: x
1460 : real(RKC) :: funcWeightedFour
1461 : funcWeightedFour = sin(Weight%omega * x) * getFunc(x)
1462 : end function
1463 :
1464 : #endif
1465 :
1466 : !%%%%%%%%%%%%%%%%%%%
1467 : #elif setChebExpan_ENABLED
1468 : !%%%%%%%%%%%%%%%%%%%
1469 :
1470 : integer(IK) :: i, j
1471 : real(RKC) :: alam, alam1, alam2, part1, part2, part3, Vector(12)
1472 : real(RKC) , parameter :: NODE(11) = [(cos(i * acos(-1._RKC) / 24._RKC), i = 1, size(NODE))]
1473 :
1474 1040 : do i = 1_IK, 12_IK
1475 960 : j = 26_IK - i
1476 960 : Vector(i) = func(i) - func(j)
1477 1040 : func(i) = func(i) + func(j)
1478 : end do
1479 80 : alam1 = Vector(1) - Vector(9)
1480 80 : alam2 = NODE(6) * (Vector(3) - Vector(7) - Vector(11))
1481 80 : cheb12(4) = alam1 + alam2
1482 80 : cheb12(10) = alam1 - alam2
1483 80 : alam1 = Vector(2) - Vector(8) - Vector(10)
1484 80 : alam2 = Vector(4) - Vector(6) - Vector(12)
1485 80 : alam = NODE(3) * alam1 + NODE(9) * alam2
1486 80 : cheb24(4) = cheb12(4) + alam
1487 80 : cheb24(22) = cheb12(4) - alam
1488 80 : alam = NODE(9) * alam1 - NODE(3) * alam2
1489 80 : cheb24(10) = cheb12(10) + alam
1490 80 : cheb24(16) = cheb12(10) - alam
1491 80 : part1 = NODE(4) * Vector(5)
1492 80 : part2 = NODE(8) * Vector(9)
1493 80 : part3 = NODE(6) * Vector(7)
1494 80 : alam1 = Vector(1) + part1 + part2
1495 80 : alam2 = NODE(2) * Vector(3) + part3 + NODE(10) * Vector(11)
1496 80 : cheb12(2) = alam1 + alam2
1497 80 : cheb12(12) = alam1 - alam2
1498 80 : alam = NODE(1) * Vector(2) + NODE(3) * Vector(4) + NODE(5) * Vector(6) + NODE(7) * Vector(8) + NODE(9) * Vector(10) + NODE(11) * Vector(12)
1499 80 : cheb24(2) = cheb12(2) + alam
1500 80 : cheb24(24) = cheb12(2) - alam
1501 80 : alam = NODE(11) * Vector(2) - NODE(9) * Vector(4) + NODE(7) * Vector(6) - NODE(5) * Vector(8) + NODE(3) * Vector(10) - NODE(1) * Vector(12)
1502 80 : cheb24(12) = cheb12(12) + alam
1503 80 : cheb24(14) = cheb12(12) - alam
1504 80 : alam1 = Vector(1) - part1 + part2
1505 80 : alam2 = NODE(10) * Vector(3) - part3 + NODE(2) * Vector(11)
1506 80 : cheb12(6) = alam1 + alam2
1507 80 : cheb12(8) = alam1 - alam2
1508 80 : alam = NODE(5) * Vector(2) - NODE(9) * Vector(4) - NODE(1) * Vector(6) - NODE(11) * Vector(8) + NODE(3) * Vector(10) + NODE(7) * Vector(12)
1509 80 : cheb24(6) = cheb12(6) + alam
1510 80 : cheb24(20) = cheb12(6) - alam
1511 80 : alam = NODE(7) * Vector(2) - NODE(3) * Vector(4) - NODE(11) * Vector(6) + NODE(1) * Vector(8) - NODE(9) * Vector(10) - NODE(5) * Vector(12)
1512 80 : cheb24(8) = cheb12(8) + alam
1513 80 : cheb24(18) = cheb12(8) - alam
1514 560 : do i = 1_IK, 6_IK
1515 480 : j = 14_IK - i
1516 480 : Vector(i) = func(i) - func(j)
1517 560 : func(i) = func(i) + func(j)
1518 : end do
1519 80 : alam1 = Vector(1) + NODE(8) * Vector(5)
1520 80 : alam2 = NODE(4) * Vector(3)
1521 80 : cheb12(3) = alam1 + alam2
1522 80 : cheb12(11) = alam1 - alam2
1523 80 : cheb12(7) = Vector(1) - Vector(5)
1524 80 : alam = NODE(2) * Vector(2) + NODE(6) * Vector(4) + NODE(10) * Vector(6)
1525 80 : cheb24(3) = cheb12(3) + alam
1526 80 : cheb24(23) = cheb12(3) - alam
1527 80 : alam = NODE(6) * (Vector(2) - Vector(4) - Vector(6))
1528 80 : cheb24(7) = cheb12(7) + alam
1529 80 : cheb24(19) = cheb12(7) - alam
1530 80 : alam = NODE(10) * Vector(2) - NODE(6) * Vector(4) + NODE(2) * Vector(6)
1531 80 : cheb24(11) = cheb12(11) + alam
1532 80 : cheb24(15) = cheb12(11) - alam
1533 320 : do i = 1_IK, 3_IK
1534 240 : j = 8_IK - i
1535 240 : Vector(i) = func(i) - func(j)
1536 320 : func(i) = func(i) + func(j)
1537 : end do
1538 80 : cheb12(5) = Vector(1) + NODE(8) * Vector(3)
1539 80 : cheb12(9) = func(1) - NODE(8) * func(3)
1540 80 : alam = NODE(4) * Vector(2)
1541 80 : cheb24(5) = cheb12(5) + alam
1542 80 : cheb24(21) = cheb12(5) - alam
1543 80 : alam = NODE(8) * func(2) - func(4)
1544 80 : cheb24(9) = cheb12(9) + alam
1545 80 : cheb24(17) = cheb12(9) - alam
1546 80 : cheb12(1) = func(1) + func(3)
1547 80 : alam = func(2) + func(4)
1548 80 : cheb24(1) = cheb12(1) + alam
1549 80 : cheb24(25) = cheb12(1) - alam
1550 80 : cheb12(13) = Vector(1) - Vector(3)
1551 80 : cheb24(13) = cheb12(13)
1552 : alam = 1._RKC / 6._RKC
1553 960 : do i = 2_IK, 12_IK
1554 960 : cheb12(i) = cheb12(i) * alam
1555 : end do
1556 80 : alam = 0.5_RKC * alam
1557 80 : cheb12(1) = cheb12(1) * alam
1558 80 : cheb12(13) = cheb12(13) * alam
1559 1920 : do i = 2_IK, 24_IK
1560 1920 : cheb24(i) = cheb24(i) * alam
1561 : end do
1562 80 : cheb24(1) = 0.5_RKC * alam * cheb24(1)
1563 80 : cheb24(25) = 0.5_RKC * alam * cheb24(25)
1564 :
1565 : #else
1566 : !%%%%%%%%%%%%%%%%%%%%%%%%
1567 : #error "Unrecognized interface."
1568 : !%%%%%%%%%%%%%%%%%%%%%%%%
1569 : #endif
1570 : #undef EPS_TABLE_START_OFFSET
1571 : #undef FINALIZE_INTEGRATION
1572 : #undef CYCLE_NEEDED
1573 : #undef TRANS_COEF
1574 : #undef QRULE_ARG
1575 : #undef HELP_ARG
|
