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