HepLib
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 
6 namespace {
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.1