htswanalysis/MACS/lib/gsl/gsl-1.11/integration/qaws.c

   1 /* integration/qaws.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <math.h>
  22 #include <float.h>
  23 #include <gsl/gsl_math.h>
  24 #include <gsl/gsl_errno.h>
  25 #include <gsl/gsl_integration.h>
  26
  27 #include "initialise.c"
  28 #include "append.c"
  29 #include "qpsrt.c"
  30 #include "util.c"
  31 #include "qc25s.c"
  32
  33 int
  34 gsl_integration_qaws (gsl_function * f,
  35                       const double a, const double b,
  36                       gsl_integration_qaws_table * t,
  37                       const double epsabs, const double epsrel,
  38                       const size_t limit,
  39                       gsl_integration_workspace * workspace,
  40                       double *result, double *abserr)
  41 {
  42   double area, errsum;
  43   double result0, abserr0;
  44   double tolerance;
  45   size_t iteration = 0;
  46   int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
  47
  48   /* Initialize results */
  49
  50   initialise (workspace, a, b);
  51
  52   *result = 0;
  53   *abserr = 0;
  54
  55   if (limit > workspace->limit)
  56     {
  57       GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
  58     }
  59
  60   if (b <= a)
  61     {
  62       GSL_ERROR ("limits must form an ascending sequence, a < b", GSL_EINVAL) ;
  63     }
  64
  65   if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
  66     {
  67       GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
  68                  GSL_EBADTOL);
  69     }
  70
  71   /* perform the first integration */
  72
  73   {
  74     double area1, area2;
  75     double error1, error2;
  76     int err_reliable1, err_reliable2;
  77     double a1 = a;
  78     double b1 = 0.5 * (a + b);
  79     double a2 = b1;
  80     double b2 = b;
  81
  82     qc25s (f, a, b, a1, b1, t, &area1, &error1, &err_reliable1);
  83     qc25s (f, a, b, a2, b2, t, &area2, &error2, &err_reliable2);
  84
  85     if (error1 > error2)
  86       {
  87         append_interval (workspace, a1, b1, area1, error1);
  88         append_interval (workspace, a2, b2, area2, error2);
  89       }
  90     else
  91       {
  92         append_interval (workspace, a2, b2, area2, error2);
  93         append_interval (workspace, a1, b1, area1, error1);
  94       }
  95
  96     result0 = area1 + area2;
  97     abserr0 = error1 + error2;
  98   }
  99
 100   /* Test on accuracy */
 101
 102   tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
 103
 104   /* Test on accuracy, use 0.01 relative error as an extra safety
 105      margin on the first iteration (ignored for subsequent iterations) */
 106
 107   if (abserr0 < tolerance && abserr0 < 0.01 * fabs(result0))
 108     {
 109       *result = result0;
 110       *abserr = abserr0;
 111
 112       return GSL_SUCCESS;
 113     }
 114   else if (limit == 1)
 115     {
 116       *result = result0;
 117       *abserr = abserr0;
 118
 119       GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
 120     }
 121
 122   area = result0;
 123   errsum = abserr0;
 124
 125   iteration = 2;
 126
 127   do
 128     {
 129       double a1, b1, a2, b2;
 130       double a_i, b_i, r_i, e_i;
 131       double area1 = 0, area2 = 0, area12 = 0;
 132       double error1 = 0, error2 = 0, error12 = 0;
 133       int err_reliable1, err_reliable2;
 134
 135       /* Bisect the subinterval with the largest error estimate */
 136
 137       retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
 138
 139       a1 = a_i;
 140       b1 = 0.5 * (a_i + b_i);
 141       a2 = b1;
 142       b2 = b_i;
 143
 144       qc25s (f, a, b, a1, b1, t, &area1, &error1, &err_reliable1);
 145       qc25s (f, a, b, a2, b2, t, &area2, &error2, &err_reliable2);
 146
 147       area12 = area1 + area2;
 148       error12 = error1 + error2;
 149
 150       errsum += (error12 - e_i);
 151       area += area12 - r_i;
 152
 153       if (err_reliable1 && err_reliable2)
 154         {
 155           double delta = r_i - area12;
 156
 157           if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
 158             {
 159               roundoff_type1++;
 160             }
 161           if (iteration >= 10 && error12 > e_i)
 162             {
 163               roundoff_type2++;
 164             }
 165         }
 166
 167       tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
 168
 169       if (errsum > tolerance)
 170         {
 171           if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
 172             {
 173               error_type = 2;   /* round off error */
 174             }
 175
 176           /* set error flag in the case of bad integrand behaviour at
 177              a point of the integration range */
 178
 179           if (subinterval_too_small (a1, a2, b2))
 180             {
 181               error_type = 3;
 182             }
 183         }
 184
 185       update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
 186
 187       retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
 188
 189       iteration++;
 190
 191     }
 192   while (iteration < limit && !error_type && errsum > tolerance);
 193
 194   *result = sum_results (workspace);
 195   *abserr = errsum;
 196
 197   if (errsum <= tolerance)
 198     {
 199       return GSL_SUCCESS;
 200     }
 201   else if (error_type == 2)
 202     {
 203       GSL_ERROR ("roundoff error prevents tolerance from being achieved",
 204                  GSL_EROUND);
 205     }
 206   else if (error_type == 3)
 207     {
 208       GSL_ERROR ("bad integrand behavior found in the integration interval",
 209                  GSL_ESING);
 210     }
 211   else if (iteration == limit)
 212     {
 213       GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
 214     }
 215   else
 216     {
 217       GSL_ERROR ("could not integrate function", GSL_EFAILED);
 218     }
 219
 220 }