HepLib
Loading...
Searching...
No Matches
Lib3_GK_vec.h
Go to the documentation of this file.
1// modified version from https://github.com/drjerry/quadpackpp
2// Adaptive GaussKronrod
3#pragma once
4#include "Lib3_GaussKronrod.h"
5
6namespace {
7
8 class FtnBase {
9 public:
10 virtual int operator() (unsigned yn, Real *y, Real *e, Real_t x) = 0;
11 const unsigned n;
12 FtnBase(unsigned yn) : n(yn) { }
13 };
14
15 class Function : public FtnBase {
16 public:
17 typedef std::function<int(unsigned yn, Real *y, Real *e, Real_t x, void *fdata)> f1Type;
18 virtual int operator() (unsigned yn, Real *y, Real *e, Real_t x) override { return function_(yn, y, e, x, fdata_); }
19 Function(unsigned yn, const f1Type & function, void *fdata)
20 : FtnBase(yn), function_(function), fdata_(fdata) { }
21 ~Function() { }
22 void * fdata() { return fdata_; }
23 private:
24 f1Type function_;
25 void * fdata_;
26 };
27
28 class GK : public GaussKronrod {
29 public:
30 typedef void (*PrintHookerType) (Real*, Real*, size_t *, void *);
31 GK(size_t n, size_t m) : GaussKronrod(m), limit(n) { }
32 int QAG(FtnBase &f, Real_t a, Real_t b, Real_t epsabs, Real_t epsrel, Real *result, Real *abserr, PrintHookerType PrintHooker);
33 private:
34 size_t limit;
35 Real rescale_error(Real err, Real_t result_abs, Real_t result_asc);
36 int qag(FtnBase &f, Real_t a, Real_t b, Real *result, Real *abserr);
37 };
38
39 Real GK::rescale_error(Real err, Real_t result_abs, Real_t result_asc) {
40 err = abs(err);
41 if (result_asc != Real(0) && err != Real(0)) {
42 // cast 1.5 as Real number
43 Real exponent = Real(3)/Real(2);
44 Real scale = pow((200 * err / result_asc), exponent);
45 if (scale < Real(1)) err = result_asc * scale;
46 else err = result_asc;
47 }
48 if (result_abs > mpfr::minval() / (50 * mpfr::machine_epsilon())) {
49 Real min_err = 50 * mpfr::machine_epsilon() * result_abs;
50 if (min_err > err) err = min_err;
51 }
52 return err;
53 }
54
55 int GK::qag(FtnBase &f, Real_t a, Real_t b, Real *result, Real *abserr) {
56 const Real center = (a + b) / 2;
57 const Real half_length = (b - a) / 2;
58 const Real abs_half_length = abs(half_length);
59 unsigned yn = f.n;
60
61 std::vector<Real> f_center(yn), ef_center(yn);
62 if (f(yn, f_center.data(), ef_center.data(), center)) return FAILURE;
63
64 std::vector<Real> result_gauss(yn), result_kronrod(yn), result_abs(yn), result_asc(yn), mean(yn), err(yn);
65 for (int j = 0; j < yn; j++) {
66 result_gauss[j] = Real(0);
67 result_kronrod[j] = f_center[j] * wgk_[n_ - 1];
68 result_abs[j] = abs(result_kronrod[j]);
69 result_asc[j] = mean[j] = err[j] = Real(0);
70 if (n_ % 2 == 0) result_gauss[j] = f_center[j] * wg_[n_ / 2 - 1];
71 }
72
73 std::vector<Real> fv1(yn * n_), e_fv1(yn * n_), fv2(yn * n_), e_fv2(yn * n_);
74 if (true) { // Parallel
75 std::vector<int> RC1(n_), RC2(n_);
76 std::fill(RC1.begin(), RC1.end(), 0);
77 std::fill(RC2.begin(), RC2.end(), 0);
78 auto prec = mpfr::mpreal::get_default_prec();
79 auto rnd = mpfr::mpreal::get_default_rnd();
80 #pragma omp parallel for
81 for (int jj = 0; jj < 2*n_; jj++) {
82 int j = jj / 2, j2 = jj % 2;
83 if (mpfr::mpreal::get_default_prec() != prec)
84 mpfr::mpreal::set_default_prec(prec);
85 mpfr::mpreal::set_default_rnd(rnd);
86 Real abscissa = half_length * xgk_[j];
87 if (j2 == 0)
88 RC1[j] = f(yn, fv1.data() + yn * j, e_fv1.data() + yn * j, center - abscissa);
89 else
90 RC2[j] = f(yn, fv2.data() + yn * j, e_fv2.data() + yn * j, center + abscissa);
91 mpfr_free_cache();
92 }
93 for (int j = 0; j < n_; j++)
94 if (RC1[j] != 0 || RC2[j] != 0) return FAILURE;
95
96 for (int jj = 0; jj < yn; jj++) { // jj for y index
97 for (int j = 0; j < (n_ - 1) / 2; j++) {
98 int jtw = j * 2 + 1;
99 Real fval1 = fv1[yn * jtw + jj];
100 Real fval2 = fv2[yn * jtw + jj];
101 Real fsum = fval1 + fval2;
102 result_gauss[jj] += wg_[j] * fsum;
103 result_kronrod[jj] += wgk_[jtw] * fsum;
104 result_abs[jj] += wgk_[jtw] * (abs(fval1) + abs(fval2));
105 }
106 for (int j = 0; j < n_ / 2; j++) {
107 int jtwm1 = j * 2;
108 Real fval1 = fv1[yn * jtwm1 + jj];
109 Real fval2 = fv2[yn * jtwm1 + jj];
110 result_kronrod[jj] += wgk_[jtwm1] * (fval1 + fval2);
111 result_abs[jj] += wgk_[jtwm1] * (abs(fval1) + abs(fval2));
112 }
113 }
114 } else { // Non-Parallel
115 Real fsum, fval1, fval2;
116 for (int j = 0; j < (n_ - 1) / 2; j++) {
117 int jtw = j * 2 + 1; // j=1,2,3 jtw=2,4,6
118 Real abscissa = half_length * xgk_[jtw];
119 f(yn, fv1.data() + yn * jtw, e_fv1.data() + yn * jtw, center - abscissa);
120 f(yn, fv2.data() + yn * jtw, e_fv2.data() + yn * jtw, center + abscissa);
121 for (int jj = 0; jj < yn; jj++) { // jj for y index
122 Real_t fval1 = fv1[yn * jtw + jj];
123 Real_t fval2 = fv2[yn * jtw + jj];
124 fsum = fval1 + fval2;
125 result_gauss[jj] += wg_[j] * fsum;
126 result_kronrod[jj] += wgk_[jtw] * fsum;
127 result_abs[jj] += wgk_[jtw] * (abs(fval1) + abs(fval2));
128 }
129 }
130 for (int j = 0; j < n_ / 2; j++) {
131 int jtwm1 = j * 2;
132 Real abscissa = half_length * xgk_[jtwm1];
133 f(yn, fv1.data() + yn * jtwm1, e_fv1.data() + yn * jtwm1, center - abscissa);
134 f(yn, fv2.data() + yn * jtwm1, e_fv2.data() + yn * jtwm1, center + abscissa);
135 for (int jj = 0; jj < yn; jj++) { // jj for y index
136 Real_t fval1 = fv1[yn * jtwm1 + jj];
137 Real_t fval2 = fv2[yn * jtwm1 + jj];
138 result_kronrod[jj] += wgk_[jtwm1] * (fval1 + fval2);
139 result_abs[jj] += wgk_[jtwm1] * (abs(fval1) + abs(fval2));
140 }
141 }
142 }
143
144 for (int jj = 0; jj < yn; jj++) { // jj for y index
145 mean[jj] = result_kronrod[jj] / 2;
146 result_asc[jj] = wgk_[n_ - 1] * abs(f_center[jj] - mean[jj]);
147 for (int j = 0; j < n_ - 1; j++) {
148 result_asc[jj] += wgk_[j] * (abs(fv1[yn * j + jj] - mean[jj]) + abs(fv2[yn * j + jj] - mean[jj]));
149 }
150 /* scale by the width of the integration region */
151 err[jj] = (result_kronrod[jj] - result_gauss[jj]) * half_length;
152 result_kronrod[jj] *= half_length;
153 result_abs[jj] *= abs_half_length;
154 result_asc[jj] *= abs_half_length;
155 result[jj] = result_kronrod[jj];
156 abserr[jj] = rescale_error(err[jj], result_abs[jj], result_asc[jj]);
157 }
158
159 return SUCCESS;
160 }
161
162 inline Real max(unsigned yn, Real *x) {
163 Real max_value = x[0];
164 for (int j = 1; j < yn; j++) if (max_value < x[j]) max_value = x[j];
165 return max_value;
166 }
167
168 int GK::QAG(FtnBase &f, Real_t a, Real_t b, Real_t epsabs, Real_t epsrel, Real *result, Real *abserr, PrintHookerType PrintHooker) {
169 size_t size = 0;
170 size_t nrmax = 0;
171 size_t i_work = 0;
172 size_t maximum_level = 0;
173
174 unsigned yn = f.n;
175 size_t n = limit;
176
177 std::vector<Real> alist(n), blist(n);
178 std::vector<std::vector<Real>> rlist(n, std::vector<Real>(yn));
179 std::vector<std::vector<Real>> elist(n, std::vector<Real>(yn));
180 std::vector<size_t> order(n), level(n);
181
182 std::vector<Real> area(yn), errsum(yn), result0(yn), abserr0(yn), tolerance(yn);
183 size_t iteration = 0;
184
185 {
186 size = 0;
187 nrmax = 0;
188 i_work = 0;
189 alist[0] = a;
190 blist[0] = b;
191 for (int j = 0; j < yn; j++) rlist[0][j] = elist[0][j] = Real(0);
192 order[0] = 0;
193 level[0] = 0;
194 maximum_level = 0;
195 }
196
197 for (int j = 0; j < yn; j++) result[j] = abserr[j] = Real(0);
198
199 if (this->qag(f, a, b, result0.data(), abserr0.data())) return FAILURE;
200
201 {
202 size = 1;
203 for (int j = 0; j < yn; j++) {
204 rlist[0][j] = result0[j];
205 elist[0][j] = abserr0[j];
206 }
207 }
208
209 /* Test on accuracy */
210 for (int j = 0; j < yn; j++) tolerance[j] = max(epsabs, epsrel * abs(result0[j]));
211 bool ok = true;
212
213 for (int j = 0; j < yn; j++) {
214 if (abserr0[j] > tolerance[j]) { ok = false; break; }
215 }
216
217 if (ok) {
218 for (int j = 0; j < yn; j++) {
219 result[j] = result0[j];
220 abserr[j] = abserr0[j];
221 }
222 size_t kk = 1;
223 auto fdata = ((Function&)f).fdata();
224 if (PrintHooker) PrintHooker(result, abserr, &kk, fdata);
225 return SUCCESS;
226 }
227
228 for (int j = 0; j < yn; j++) {
229 area[j] = result0[j];
230 errsum[j] = abserr0[j];
231 }
232 iteration = 1;
233
234 do {
235 Real a1, b1, a2, b2;
236 Real a_i, b_i, r_i[yn], e_i[yn];
237 std::vector<Real> area1(yn), area2(yn), area12(yn);
238 std::vector<Real> error1(yn), error2(yn), error12(yn);
239
240 for (int j = 0; j < yn; j++) {
241 area1[j] = area2[j] = area12[j] = error1[j] = error2[j] = error12[j] = Real(0);
242 }
243
244 /* Bisect the subinterval with the largest error estimate */
245 {
246 a_i = alist[i_work];
247 b_i = blist[i_work];
248 for (int j = 0; j < yn; j++) {
249 r_i[j] = rlist[i_work][j];
250 e_i[j] = elist[i_work][j];
251 }
252 }
253
254 a1 = a_i;
255 b1 = (a_i + b_i) / Real(2);
256 a2 = b1;
257 b2 = b_i;
258
259 if (this->qag(f, a1, b1, area1.data(), error1.data())) return FAILURE;
260 if (this->qag(f, a2, b2, area2.data(), error2.data())) return FAILURE;
261
262 for (int j = 0; j < yn; j++) {
263 area12[j] = area1[j] + area2[j];
264 error12[j] = error1[j] + error2[j];
265 errsum[j] += (error12[j] - e_i[j]);
266 area[j] += area12[j] - r_i[j];
267 tolerance[j] = max(epsabs, epsrel * abs(area[j]));
268 }
269
270 {
271 const size_t i_max = i_work;
272 const size_t i_new = size;
273 const size_t new_level = level[i_max] + 1;
274
275 /* append the newly-created intervals to the list */
276 if (max(yn, error2.data()) > max(yn, error1.data())) {
277 alist[i_max] = a2; /* blist[maxerr] is already == b2 */
278 for (int j = 0; j < yn; j++) {
279 rlist[i_max][j] = area2[j];
280 elist[i_max][j] = error2[j];
281 }
282 level[i_max] = new_level;
283
284 alist[i_new] = a1;
285 blist[i_new] = b1;
286 for (int j = 0; j < yn; j++) {
287 rlist[i_new][j] = area1[j];
288 elist[i_new][j] = error1[j];
289 }
290 level[i_new] = new_level;
291 } else {
292 blist[i_max] = b1; /* alist[maxerr] is already == a1 */
293 for (int j = 0; j < yn; j++) {
294 rlist[i_max][j] = area1[j];
295 elist[i_max][j] = error1[j];
296 }
297 level[i_max] = new_level;
298
299 alist[i_new] = a2;
300 blist[i_new] = b2;
301 for (int j = 0; j < yn; j++) {
302 rlist[i_new][j] = area2[j];
303 elist[i_new][j] = error2[j];
304 }
305 level[i_new] = new_level;
306 }
307
308 size++;
309 if (new_level > maximum_level) maximum_level = new_level;
310
311 {
312 const size_t last = size - 1;
313 Real errmax, errmin;
314 int i, k, top;
315 size_t i_nrmax = nrmax;
316 size_t i_maxerr = order[i_nrmax];
317
318 /* Check whether the list contains more than two error estimates */
319 if (last < 2) {
320 order[0] = 0;
321 order[1] = 1;
322 i_work = i_maxerr;
323 } else {
324 errmax = max(yn, elist[i_maxerr].data());
325 while (i_nrmax > 0 && errmax > max(yn, elist[order[i_nrmax - 1]].data())) {
326 order[i_nrmax] = order[i_nrmax - 1];
327 i_nrmax--;
328 }
329 if (last < (limit / 2 + 2)) top = last;
330 else top = limit - last + 1;
331 i = i_nrmax + 1;
332 while (i < top && errmax < max(yn, elist[order[i]].data())) {
333 order[i - 1] = order[i];
334 i++;
335 }
336 order[i - 1] = i_maxerr;
337 errmin = max(yn, elist[last].data());
338 k = top - 1;
339 while (k > i - 2 && errmin >= max(yn, elist[order[k]].data())) {
340 order[k + 1] = order[k];
341 k--;
342 }
343 order[k + 1] = last;
344
345 /* Set i_max and e_max */
346 i_maxerr = order[i_nrmax];
347 i_work = i_maxerr;
348 nrmax = i_nrmax;
349 }
350 }
351 }
352
353 {
354 a_i = alist[i_work];
355 b_i = blist[i_work];
356 for (int j = 0; j < yn; j++) {
357 r_i[j] = rlist[i_work][j];
358 e_i[j] = elist[i_work][j];
359 }
360 }
361
362 if (PrintHooker) {
363 size_t kk = iteration;
364 auto fdata = ((Function&)f).fdata();
365 PrintHooker(area.data(), errsum.data(), &kk, fdata);
366 if (kk != iteration) {
367 for (int j = 0; j < yn; j++) {
368 result[j] = area[j];
369 abserr[j] = errsum[j];
370 }
371 return SUCCESS;
372 }
373 }
374
375 ok = true;
376 for (int j = 0; j < yn; j++) {
377 if (errsum[j] > tolerance[j]) { ok = false; break; }
378 }
379 iteration++;
380 } while (iteration < limit && !ok);
381
382 for (int j = 0; j < yn; j++) {
383 result[j] = Real(0);
384 for (size_t k = 0; k < size; k++) result[j] += rlist[k][j];
385 abserr[j] = errsum[j];
386 }
387 return SUCCESS;
388 }
389}
a
int * a
Definition Functions.cpp:234
Lib3_GaussKronrod.h
FAILURE
#define FAILURE
Definition Lib3_GaussKronrod.h:9
SUCCESS
#define SUCCESS
Definition Lib3_GaussKronrod.h:8
Generated by doxygen 1.9.8