inte_singular_gsl.h
Go to the documentation of this file.
1 /*
2  -------------------------------------------------------------------
3 
4  Copyright (C) 2006-2017, Jerry Gagelman and Andrew W. Steiner
5 
6  This file is part of O2scl.
7 
8  O2scl is free software; you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version.
12 
13  O2scl is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with O2scl. If not, see <http://www.gnu.org/licenses/>.
20 
21  -------------------------------------------------------------------
22 */
23 #ifndef O2SCL_GSL_INTE_SINGULAR_H
24 #define O2SCL_GSL_INTE_SINGULAR_H
25 
26 /** \file inte_singular_gsl.h
27  \brief File defining \ref o2scl::inte_singular_gsl
28 */
29 
30 #include <cmath>
31 
32 #include <gsl/gsl_integration.h>
33 
34 #include <o2scl/string_conv.h>
35 #include <o2scl/inte.h>
36 #include <o2scl/inte_kronrod_gsl.h>
37 
38 #ifndef DOXYGEN_NO_O2NS
39 namespace o2scl {
40 #endif
41 
42  /** \brief Base class for integrating a function with a
43  singularity (GSL)
44 
45  This class contains the extrapolation table mechanics and the
46  base integration function for singular integrals from GSL. The
47  casual end-user should use the classes explained in the
48  \ref inte_section section of the User's guide.
49 
50  \future Some of the functions inside this class could
51  be moved out of header files?
52  */
53  template<class func_t> class inte_singular_gsl :
54  public inte_kronrod_gsl<func_t> {
55 
56  public:
57 
58  /** \brief A structure for extrapolation for \ref inte_qags_gsl
59 
60  \future Move this to a new class, with qelg() as a method
61  */
62  typedef struct extrapolation_table {
63  /// Index of new element in the first column
64  size_t n;
65  /// Lower diagonals of the triangular epsilon table
66  double rlist2[52];
67  /// Number of calls
68  size_t nres;
69  /// Three most recent results
70  double res3la[3];
71  } extrap_table;
72 
73  protected:
74 
75  /// Initialize the table
77  table->n = 0;
78  table->nres = 0;
79  return;
80  }
81 
82  /// Append a result to the table
83  void append_table(struct extrapolation_table *table, double y) {
84  size_t n;
85  n = table->n;
86  table->rlist2[n] = y;
87  table->n++;
88  return;
89  }
90 
91  /** \brief Test if the integrand satisfies \f$ f = |f| \f$
92  */
93  inline int test_positivity(double result, double resabs) {
94  double dbl_eps=std::numeric_limits<double>::epsilon();
95  int status=(fabs(result) >= (1-50*dbl_eps)*resabs);
96  return status;
97  }
98 
99  /** \brief Determines the limit of a given sequence
100  of approximations
101 
102  For certain convergent series \f$ \sum_k a_k \f$ whose error
103  term \f$ E_n = \sum_{k=n}^\infty a_k \f$ is well behaved, it
104  is possible to find a transformation of the sequence that
105  yields a faster converging series to the same limit. This
106  method of extrapolation applies to some sequences of
107  adaptive-approximation and error-estimation for numerical
108  integration. This function implements the \f$
109  \varepsilon\f$-algorithm (\ref Wynn56, \ref Piessens83) for an
110  extrapolation table stored in \c table.
111 
112  Quadpack documentation
113  \verbatim
114  c
115  c list of major variables
116  c -----------------------
117  c e0 - the 4 elements on which the computation of a new
118  c e1 element in the epsilon table is based
119  c e2
120  c e3 e0
121  c e3 e1 new
122  c e2
123  c newelm - number of elements to be computed in the new
124  c diagonal
125  c error - error = abs(e1-e0)+abs(e2-e1)+abs(new-e2)
126  c result - the element in the new diagonal with least value
127  c of error
128  c
129  c machine dependent constants
130  c ---------------------------
131  c
132  c epmach is the largest relative spacing.
133  c oflow is the largest positive magnitude.
134  c limexp is the maximum number of elements the epsilon
135  c table can contain. if this number is reached, the upper
136  c diagonal of the epsilon table is deleted.
137  c
138  \endverbatim
139  */
140  void qelg(struct extrapolation_table *table, double *result,
141  double *abserr) {
142 
143  double *epstab = table->rlist2;
144  double *res3la = table->res3la;
145  const size_t n = table->n - 1;
146 
147  const double current = epstab[n];
148 
149  double absolute = GSL_DBL_MAX;
150  double relative = 5 * GSL_DBL_EPSILON * fabs (current);
151 
152  const size_t newelm = n / 2;
153  const size_t n_orig = n;
154  size_t n_final = n;
155  size_t i;
156 
157  const size_t nres_orig = table->nres;
158 
159  *result = current;
160  *abserr = GSL_DBL_MAX;
161 
162  if (n < 2) {
163  *result = current;
164  *abserr = GSL_MAX_DBL (absolute, relative);
165  return;
166  }
167 
168  epstab[n + 2] = epstab[n];
169  epstab[n] = GSL_DBL_MAX;
170 
171  for (i = 0; i < newelm; i++) {
172  double res = epstab[n - 2 * i + 2];
173  double e0 = epstab[n - 2 * i - 2];
174  double e1 = epstab[n - 2 * i - 1];
175  double e2 = res;
176 
177  double e1abs = fabs (e1);
178  double delta2 = e2 - e1;
179  double err2 = fabs (delta2);
180  double tol2 = GSL_MAX_DBL (fabs (e2), e1abs) * GSL_DBL_EPSILON;
181  double delta3 = e1 - e0;
182  double err3 = fabs (delta3);
183  double tol3 = GSL_MAX_DBL (e1abs, fabs (e0)) * GSL_DBL_EPSILON;
184 
185  double e3, delta1, err1, tol1, ss;
186 
187  if (err2 <= tol2 && err3 <= tol3) {
188  /* If e0, e1 and e2 are equal to within machine accuracy,
189  convergence is assumed. */
190 
191  *result = res;
192  absolute = err2 + err3;
193  relative = 5 * GSL_DBL_EPSILON * fabs (res);
194  *abserr = GSL_MAX_DBL (absolute, relative);
195  return;
196  }
197 
198  e3 = epstab[n - 2 * i];
199  epstab[n - 2 * i] = e1;
200  delta1 = e1 - e3;
201  err1 = fabs (delta1);
202  tol1 = GSL_MAX_DBL (e1abs, fabs (e3)) * GSL_DBL_EPSILON;
203 
204  /* If two elements are very close to each other, omit a part of
205  the table by adjusting the value of n */
206 
207  if (err1 <= tol1 || err2 <= tol2 || err3 <= tol3) {
208  n_final = 2 * i;
209  break;
210  }
211 
212  ss = (1 / delta1 + 1 / delta2) - 1 / delta3;
213 
214  /* Test to detect irregular behaviour in the table, and
215  eventually omit a part of the table by adjusting the value of
216  n. */
217  if (fabs (ss * e1) <= 0.0001) {
218  n_final = 2 * i;
219  break;
220  }
221  /* Compute a new element and eventually adjust the value of
222  result. */
223 
224  res = e1 + 1 / ss;
225  epstab[n - 2 * i] = res;
226 
227  {
228  const double error = err2 + fabs (res - e2) + err3;
229 
230  if (error <= *abserr) {
231  *abserr = error;
232  *result = res;
233  }
234  }
235  }
236 
237  /* Shift the table */
238 
239  {
240  const size_t limexp = 50 - 1;
241 
242  if (n_final == limexp) {
243  n_final = 2 * (limexp / 2);
244  }
245  }
246 
247  if (n_orig % 2 == 1) {
248  for (i = 0; i <= newelm; i++) {
249  epstab[1 + i * 2] = epstab[i * 2 + 3];
250  }
251  } else {
252  for (i = 0; i <= newelm; i++) {
253  epstab[i * 2] = epstab[i * 2 + 2];
254  }
255  }
256  if (n_orig != n_final) {
257  for (i = 0; i <= n_final; i++) {
258  epstab[i] = epstab[n_orig - n_final + i];
259  }
260  }
261 
262  table->n = n_final + 1;
263 
264  if (nres_orig < 3) {
265 
266  res3la[nres_orig] = *result;
267  *abserr = GSL_DBL_MAX;
268 
269  } else {
270  /* Compute error estimate */
271  *abserr = (fabs (*result - res3la[2]) + fabs (*result - res3la[1])
272  + fabs (*result - res3la[0]));
273 
274  res3la[0] = res3la[1];
275  res3la[1] = res3la[2];
276  res3la[2] = *result;
277  }
278 
279  /* In QUADPACK the variable table->nres is incremented at the top of
280  qelg, so it increases on every call. This leads to the array
281  res3la being accessed when its elements are still undefined, so I
282  have moved the update to this point so that its value more
283  useful. */
284 
285  table->nres = nres_orig + 1;
286 
287  *abserr = GSL_MAX_DBL (*abserr, 5 * GSL_DBL_EPSILON * fabs (*result));
288 
289  return;
290  }
291 
292  /// Determine if an interval is large
293  int large_interval (inte_workspace_gsl * workspace) {
294  size_t i = workspace->i ;
295  const size_t * level = workspace->level;
296 
297  if (level[i] < workspace->maximum_level) {
298  return 1;
299  } else {
300  return 0;
301  }
302  }
303 
304  /// Reset workspace to work on the interval with the largest error
305  inline void reset_nrmax (inte_workspace_gsl * workspace) {
306  workspace->nrmax = 0;
307  workspace->i = workspace->order[0] ;
308  }
309 
310  /// Increase workspace
311  int increase_nrmax (inte_workspace_gsl * workspace) {
312  int k;
313  int id = workspace->nrmax;
314  int jupbnd;
315 
316  const size_t * level = workspace->level;
317  const size_t * order = workspace->order;
318 
319  size_t limit = workspace->limit ;
320  size_t last = workspace->size - 1 ;
321 
322  if (last > (1 + limit / 2)) {
323  jupbnd = limit + 1 - last;
324  } else {
325  jupbnd = last;
326  }
327 
328  for (k = id; k <= jupbnd; k++) {
329  size_t i_max = order[workspace->nrmax];
330 
331  workspace->i = i_max ;
332 
333  if (level[i_max] < workspace->maximum_level) {
334  return 1;
335  }
336 
337  workspace->nrmax++;
338 
339  }
340  return 0;
341  }
342 
343  /** \brief Integration function
344 
345  \future Remove goto statements. Before this is done, it might
346  be best to add some tests which fail in the various ways.
347  */
348  int qags(func_t &func, const double a, const double b,
349  const double l_epsabs, const double l_epsrel,
350  double *result, double *abserr) {
351 
352  double area, errsum;
353  double res_ext, err_ext;
354  double result0, abserr0, resabs0, resasc0;
355  double tolerance;
356 
357  double ertest = 0;
358  double error_over_large_intervals = 0;
359  double reseps = 0, abseps = 0, correc = 0;
360  size_t ktmin = 0;
361  int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
362  int error_type = 0, error_type2 = 0;
363 
364  size_t iteration = 0;
365 
366  int positive_integrand = 0;
367  int extrapolate = 0;
368  int disallow_extrapolation = 0;
369 
370  struct extrapolation_table table;
371 
372  /* Initialize results */
373 
374  this->w->initialise(a,b);
375 
376  *result = 0;
377  *abserr = 0;
378 
379  size_t limit=this->w->limit;
380 
381  /* Test on accuracy */
382 
383  if (this->tol_abs <= 0 && (this->tol_rel < 50 * GSL_DBL_EPSILON ||
384  this->tol_rel < 0.5e-28)) {
385  this->last_iter=0;
386 
387  std::string estr="Tolerance cannot be achieved with given ";
388  estr+="value of tol_abs, "+dtos(l_epsabs)+", and tol_rel, "+
389  dtos(l_epsrel)+", in inte_singular_gsl::qags().";
390  O2SCL_ERR(estr.c_str(),exc_ebadtol);
391  }
392 
393  /* Perform the first integration */
394 
395  this->gauss_kronrod(func,a,b,&result0,&abserr0,&resabs0,&resasc0);
396 
397  this->w->set_initial_result (result0, abserr0);
398 
399  tolerance = GSL_MAX_DBL (this->tol_abs, this->tol_rel * fabs (result0));
400 
401  if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 &&
402  abserr0 > tolerance) {
403 
404  *result = result0;
405  *abserr = abserr0;
406 
407  this->last_iter=1;
408 
409  std::string estr="Cannot reach tolerance because of roundoff error ";
410  estr+="on first attempt in inte_singular_gsl::qags().";
411  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
412 
413  } else if ((abserr0 <= tolerance &&
414  abserr0 != resasc0) || abserr0 == 0.0) {
415 
416  *result = result0;
417  *abserr = abserr0;
418  this->last_iter=1;
419  return success;
420 
421  } else if (limit == 1) {
422 
423  *result = result0;
424  *abserr = abserr0;
425 
426  this->last_iter=1;
427  O2SCL_CONV2_RET("A maximum of 1 iteration was insufficient ",
428  "in inte_singular_gsl::qags().",
429  exc_emaxiter,this->err_nonconv);
430  }
431 
432  /* Initialization */
433 
434  initialise_table (&table);
435  append_table (&table, result0);
436 
437  area = result0;
438  errsum = abserr0;
439 
440  res_ext = result0;
441  err_ext = GSL_DBL_MAX;
442 
443  positive_integrand = this->test_positivity (result0, resabs0);
444 
445  iteration = 1;
446 
447  do {
448 
449  // Output iteration information
450  if (this->verbose>0) {
451  std::cout << this->type();
452  std::cout << " Iter: " << iteration;
453  std::cout.setf(std::ios::showpos);
454  std::cout << " Res: " << area;
455  std::cout.unsetf(std::ios::showpos);
456  std::cout << " Err: " << errsum
457  << " Tol: " << tolerance << std::endl;
458  if (this->verbose>1) {
459  char ch;
460  std::cout << "Press a key and type enter to continue. " ;
461  std::cin >> ch;
462  }
463  }
464 
465  size_t current_level;
466  double a1, b1, a2, b2;
467  double a_i, b_i, r_i, e_i;
468  double area1 = 0, area2 = 0, area12 = 0;
469  double error1 = 0, error2 = 0, error12 = 0;
470  double resasc1, resasc2;
471  double resabs1, resabs2;
472  double last_e_i;
473 
474  /* Bisect the subinterval with the largest error estimate */
475 
476  this->w->retrieve (&a_i, &b_i, &r_i, &e_i);
477 
478  current_level = this->w->level[this->w->i] + 1;
479 
480  a1 = a_i;
481  b1 = 0.5 * (a_i + b_i);
482  a2 = b1;
483  b2 = b_i;
484 
485  iteration++;
486 
487  this->gauss_kronrod(func,a1,b1,&area1,&error1,&resabs1,&resasc1);
488  this->gauss_kronrod(func,a2,b2,&area2,&error2,&resabs2,&resasc2);
489 
490  area12 = area1 + area2;
491  error12 = error1 + error2;
492  last_e_i = e_i;
493 
494  /* Improve previous approximations to the integral and test for
495  accuracy.
496 
497  We write these expressions in the same way as the original
498  QUADPACK code so that the rounding errors are the same, which
499  makes testing easier.
500  */
501 
502  errsum = errsum + error12 - e_i;
503  area = area + area12 - r_i;
504 
505  tolerance = GSL_MAX_DBL (this->tol_abs, this->tol_rel * fabs (area));
506  if (resasc1 != error1 && resasc2 != error2) {
507  double delta = r_i - area12;
508 
509  if (fabs (delta) <= 1.0e-5 * fabs (area12) &&
510  error12 >= 0.99 * e_i) {
511  if (!extrapolate) {
512  roundoff_type1++;
513  } else {
514  roundoff_type2++;
515  }
516  }
517  if (iteration > 10 && error12 > e_i) {
518  roundoff_type3++;
519  }
520  }
521 
522  // Test for roundoff and eventually set error flag
523 
524  if (roundoff_type1 + roundoff_type2 >= 10 ||
525  roundoff_type3 >= 20) {
526  /* round off error */
527  error_type = 2;
528  }
529 
530  if (roundoff_type2 >= 5) {
531  error_type2 = 1;
532  }
533 
534  // Set error flag in the case of bad integrand behaviour at
535  // a point of the integration range
536 
537  if (this->w->subinterval_too_small (a1, a2, b2)) {
538  error_type = 4;
539  }
540 
541  /* append the newly-created intervals to the list */
542 
543  this->w->update(a1,b1,area1,error1,a2,b2,area2,error2);
544 
545  if (errsum <= tolerance) {
546 
547  // Output final iteration information
548  if (this->verbose>0) {
549  std::cout << this->type();
550  std::cout << " Iter: " << iteration;
551  std::cout.setf(std::ios::showpos);
552  std::cout << " Res: " << area;
553  std::cout.unsetf(std::ios::showpos);
554  std::cout << " Err: " << errsum
555  << " Tol: " << tolerance << std::endl;
556  if (this->verbose>1) {
557  char ch;
558  std::cout << "Press a key and type enter to continue. " ;
559  std::cin >> ch;
560  }
561  }
562 
563  goto compute_result;
564  }
565 
566  if (error_type) {
567  break;
568  }
569 
570  if (iteration >= limit - 1) {
571  error_type = 1;
572  break;
573  }
574 
575  if (iteration == 2) {
576  error_over_large_intervals = errsum;
577  ertest = tolerance;
578  append_table (&table,area);
579  continue;
580  }
581 
582  if (disallow_extrapolation) {
583  continue;
584  }
585 
586  error_over_large_intervals += -last_e_i;
587 
588  if (current_level < this->w->maximum_level) {
589  error_over_large_intervals += error12;
590  }
591 
592  if (!extrapolate) {
593 
594  /* test whether the interval to be bisected next is the
595  smallest interval. */
596 
597  if (large_interval (this->w)) {
598  continue;
599  }
600 
601  extrapolate = 1;
602  this->w->nrmax = 1;
603  }
604 
605  if (!error_type2 && error_over_large_intervals > ertest) {
606  if (increase_nrmax (this->w)) {
607  continue;
608  }
609  }
610 
611  // Perform extrapolation
612 
613  append_table(&table,area);
614 
615  qelg(&table,&reseps,&abseps);
616 
617  ktmin++;
618 
619  if (ktmin > 5 && err_ext < 0.001 * errsum) {
620  error_type = 5;
621  }
622 
623  if (abseps < err_ext) {
624  ktmin = 0;
625  err_ext = abseps;
626  res_ext = reseps;
627  correc = error_over_large_intervals;
628  ertest = GSL_MAX_DBL (this->tol_abs,
629  this->tol_rel * fabs (reseps));
630  if (err_ext <= ertest) {
631  break;
632  }
633  }
634 
635  /* Prepare bisection of the smallest interval. */
636 
637  if (table.n == 1) {
638  disallow_extrapolation = 1;
639  }
640 
641  if (error_type == 5) {
642  break;
643  }
644 
645  /* work on interval with largest error */
646 
647  reset_nrmax (this->w);
648  extrapolate = 0;
649  error_over_large_intervals = errsum;
650 
651  } while (iteration < limit);
652 
653  *result = res_ext;
654  *abserr = err_ext;
655 
656  if (err_ext == GSL_DBL_MAX)
657  goto compute_result;
658 
659  if (error_type || error_type2) {
660  if (error_type2) {
661  err_ext += correc;
662  }
663 
664  if (error_type == 0)
665  error_type = 3;
666 
667  if (res_ext != 0.0 && area != 0.0) {
668  if (err_ext / fabs (res_ext) > errsum / fabs (area))
669  goto compute_result;
670  } else if (err_ext > errsum) {
671  goto compute_result;
672  } else if (area == 0.0) {
673  goto return_error;
674  }
675  }
676 
677  /* Test on divergence. */
678 
679  {
680  double max_area = GSL_MAX_DBL (fabs (res_ext),fabs (area));
681 
682  if (!positive_integrand && max_area < 0.01 * resabs0)
683  goto return_error;
684  }
685 
686  {
687  double ratio = res_ext / area;
688 
689  if (ratio < 0.01 || ratio > 100.0 || errsum > fabs (area))
690  error_type = 6;
691  }
692 
693  goto return_error;
694 
695  compute_result:
696 
697  *result = this->w->sum_results();
698  *abserr = errsum;
699 
700  return_error:
701 
702  if (error_type > 2) {
703  error_type--;
704  }
705 
706  this->last_iter=iteration;
707 
708  if (error_type == 0) {
709  return success;
710  } else if (error_type == 1) {
711  std::string estr="Maximum number of subdivisions ("+itos(iteration);
712  estr+=") reached in inte_singular_gsl::qags().";
713  O2SCL_CONV_RET(estr.c_str(),exc_emaxiter,this->err_nonconv);
714  } else if (error_type == 2) {
715  std::string estr="Roundoff error prevents tolerance ";
716  estr+="from being achieved in inte_singular_gsl::qags().";
717  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
718  } else if (error_type == 3) {
719  std::string estr="Bad integrand behavior ";
720  estr+="in inte_singular_gsl::qags().";
721  O2SCL_CONV_RET(estr.c_str(),exc_esing,this->err_nonconv);
722  } else if (error_type == 4) {
723  std::string estr="Roundoff error detected in extrapolation table ";
724  estr+="in inte_singular_gsl::qags().";
725  O2SCL_CONV_RET(estr.c_str(),exc_eround,this->err_nonconv);
726  } else if (error_type == 5) {
727  std::string estr="Integral is divergent or slowly convergent ";
728  estr+="in inte_singular_gsl::qags().";
729  O2SCL_CONV_RET(estr.c_str(),exc_ediverge,this->err_nonconv);
730  }
731 
732  std::string estr="Could not integrate function in inte_kronrod_gsl";
733  estr+="::qags() (it may have returned a non-finite result).";
734  O2SCL_ERR(estr.c_str(),exc_efailed);
735 
736  return exc_efailed;
737  }
738 
739  };
740 
741  /** \brief Integrate a function with a singularity (GSL)
742  [abstract base]
743 
744  This class contains the GSL-based integration function for
745  applying transformations to the user-defined integrand. The
746  casual end-user should use the classes explained in the
747  \ref inte_section section of the User's guide.
748  */
749  template<class func_t=funct11> class inte_transform_gsl :
750  public inte_singular_gsl<func_t> {
751 
752  public:
753 
754  /// The transformation to apply to the user-supplied function
755  virtual double transform(double t, func_t &func)=0;
756 
757  /** \brief Integration wrapper for internal transformed function
758  type
759  */
760  virtual void gauss_kronrod
761  (func_t &func, double a, double b,
762  double *result, double *abserr, double *resabs, double *resasc) {
763 
764  funct11 fmp=std::bind(std::mem_fn<double(double,func_t &)>
766  this,std::placeholders::_1,func);
767 
768  return this->gauss_kronrod_base
769  (fmp,a,b,result,abserr,resabs,resasc);
770  }
771 
772  };
773 
774 #ifndef DOXYGEN_NO_O2NS
775 }
776 #endif
777 
778 #endif
Integration workspace for the GSL integrators.
void initialise_table(struct extrapolation_table *table)
Initialize the table.
The main O<span style=&#39;position: relative; top: 0.3em; font-size: 0.8em&#39;>2</span>scl O$_2$scl names...
Definition: anneal.h:42
#define O2SCL_CONV_RET(d, n, b)
Set a "convergence" error and return the error value.
Definition: err_hnd.h:292
int qags(func_t &func, const double a, const double b, const double l_epsabs, const double l_epsrel, double *result, double *abserr)
Integration function.
std::function< double(double)> funct11
One-dimensional function typedef.
Definition: funct.h:44
Data table table class.
Definition: table.h:49
size_t * level
Numbers of subdivisions made.
void append_table(struct extrapolation_table *table, double y)
Append a result to the table.
struct o2scl::inte_singular_gsl::extrapolation_table extrap_table
A structure for extrapolation for inte_qags_gsl.
apparent singularity detected
Definition: err_hnd.h:93
exceeded max number of iterations
Definition: err_hnd.h:73
Base class for integrating a function with a singularity (GSL)
int large_interval(inte_workspace_gsl *workspace)
Determine if an interval is large.
#define O2SCL_CONV2_RET(d, d2, n, b)
Set an error and return the error value, two-string version.
Definition: err_hnd.h:298
size_t * order
Linear ordering vector for sort routine.
double rlist2[52]
Lower diagonals of the triangular epsilon table.
size_t last_iter
The most recent number of iterations taken.
Definition: inte.h:63
bool err_nonconv
If true, call the error handler if the routine does not converge or reach the desired tolerance (defa...
Definition: inte.h:81
int update(double a1, double b1, double area1, double error1, double a2, double b2, double area2, double error2)
Determine which new subinterval to add to the workspace stack and perform update. ...
size_t i
Index of current subinterval.
virtual void gauss_kronrod(func_t &func, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
Integration wrapper for user-specified function type.
generic failure
Definition: err_hnd.h:61
int subinterval_too_small(double a1, double a2, double b2)
Test whether the proposed subdivision falls before floating-point precision.
double tol_abs
The maximum absolute uncertainty in the value of the integral (default )
Definition: inte.h:73
int test_positivity(double result, double resabs)
Test if the integrand satisfies .
int initialise(double a, double b)
Initialize the workspace for an integration with limits a and b.
double res3la[3]
Three most recent results.
int set_initial_result(double result, double error)
Update the workspace with the result and error from the first integration.
std::string dtos(double x, int prec=6, bool auto_prec=false)
Convert a double to a string.
#define O2SCL_ERR(d, n)
Set an error with message d and code n.
Definition: err_hnd.h:273
size_t size
Current number of subintervals being used.
Integrate a function with a singularity (GSL) [abstract base].
Basic Gauss-Kronrod integration class (GSL)
size_t nrmax
Counter for extrapolation routine.
void qelg(struct extrapolation_table *table, double *result, double *abserr)
Determines the limit of a given sequence of approximations.
inte_workspace_gsl * w
The integration workspace.
void gauss_kronrod_base(func2_t &func, double a, double b, double *result, double *abserr, double *resabs, double *resasc)
The base Gauss-Kronrod integration function template.
double sum_results()
Add up all of the contributions to construct the final result.
double tol_rel
The maximum relative uncertainty in the value of the integral (default )
Definition: inte.h:68
void reset_nrmax(inte_workspace_gsl *workspace)
Reset workspace to work on the interval with the largest error.
int increase_nrmax(inte_workspace_gsl *workspace)
Increase workspace.
A structure for extrapolation for inte_qags_gsl.
size_t n
Index of new element in the first column.
size_t limit
Maximum number of subintervals allocated.
int retrieve(double *a, double *b, double *r, double *e) const
Retrieve the ith result from the workspace stack.
size_t maximum_level
Depth of subdivisions reached.
std::string itos(int x)
Convert an integer to a string.
user specified an invalid tolerance
Definition: err_hnd.h:77
Success.
Definition: err_hnd.h:47
int verbose
Verbosity.
Definition: inte.h:60
failed because of roundoff error
Definition: err_hnd.h:87
integral or series is divergent
Definition: err_hnd.h:95
virtual const char * type()
Return string denoting type ("inte")
Definition: inte.h:110

Documentation generated with Doxygen. Provided under the GNU Free Documentation License (see License Information).