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

   1 /* integration/qagp.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 <stdlib.h>
  22 #include <gsl/gsl_errno.h>
  23 #include <gsl/gsl_integration.h>
  24
  25 static int
  26 qagp (const gsl_function *f,
  27       const double *pts, const size_t npts,
  28       const double epsabs, const double epsrel, const size_t limit,
  29       gsl_integration_workspace * workspace,
  30       double *result, double *abserr,
  31       gsl_integration_rule * q);
  32
  33 #include "initialise.c"
  34 #include "qpsrt.c"
  35 #include "util.c"
  36 #include "append.c"
  37 #include "reset.c"
  38 #include "qelg.c"
  39 #include "qpsrt2.c"
  40 #include "ptsort.c"
  41 #include "positivity.c"
  42
  43 int
  44 gsl_integration_qagp (const gsl_function *f,
  45                       double * pts, size_t npts,
  46                       double epsabs, double epsrel, size_t limit,
  47                       gsl_integration_workspace * workspace,
  48                       double * result, double * abserr)
  49 {
  50   int status = qagp (f, pts, npts,
  51                      epsabs, epsrel, limit,
  52                      workspace,
  53                      result, abserr,
  54                      &gsl_integration_qk21) ;
  55
  56   return status ;
  57 }
  58
  59
  60 static int
  61 qagp (const gsl_function * f,
  62       const double *pts, const size_t npts,
  63       const double epsabs, const double epsrel,
  64       const size_t limit,
  65       gsl_integration_workspace * workspace,
  66       double *result, double *abserr,
  67       gsl_integration_rule * q)
  68 {
  69   double area, errsum;
  70   double res_ext, err_ext;
  71   double result0, abserr0, resabs0;
  72   double tolerance;
  73
  74   double ertest = 0;
  75   double error_over_large_intervals = 0;
  76   double reseps = 0, abseps = 0, correc = 0;
  77   size_t ktmin = 0;
  78   int roundoff_type1 = 0, roundoff_type2 = 0, roundoff_type3 = 0;
  79   int error_type = 0, error_type2 = 0;
  80
  81   size_t iteration = 0;
  82
  83   int positive_integrand = 0;
  84   int extrapolate = 0;
  85   int disallow_extrapolation = 0;
  86
  87   struct extrapolation_table table;
  88
  89   const size_t nint = npts - 1; /* number of intervals */
  90
  91   size_t *ndin = workspace->level; /* temporarily alias ndin to level */
  92
  93   size_t i;
  94
  95   /* Initialize results */
  96
  97   *result = 0;
  98   *abserr = 0;
  99
 100   /* Test on validity of parameters */
 101
 102   if (limit > workspace->limit)
 103     {
 104       GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
 105     }
 106
 107   if (npts > workspace->limit)
 108     {
 109       GSL_ERROR ("npts exceeds size of workspace", GSL_EINVAL);
 110     }
 111
 112   if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
 113     {
 114       GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
 115                  GSL_EBADTOL);
 116     }
 117
 118   /* Check that the integration range and break points are an
 119      ascending sequence */
 120
 121   for (i = 0; i < nint; i++)
 122     {
 123       if (pts[i + 1] < pts[i])
 124         {
 125           GSL_ERROR ("points are not in an ascending sequence", GSL_EINVAL);
 126         }
 127     }
 128
 129   /* Perform the first integration */
 130
 131   result0 = 0;
 132   abserr0 = 0;
 133   resabs0 = 0;
 134
 135   initialise (workspace, 0.0, 0.0) ;
 136
 137   for (i = 0; i < nint; i++)
 138     {
 139       double area1, error1, resabs1, resasc1;
 140       const double a1 = pts[i];
 141       const double b1 = pts[i + 1];
 142
 143       q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
 144
 145       result0 = result0 + area1;
 146       abserr0 = abserr0 + error1;
 147       resabs0 = resabs0 + resabs1;
 148
 149       append_interval (workspace, a1, b1, area1, error1);
 150
 151       if (error1 == resasc1 && error1 != 0.0)
 152         {
 153           ndin[i] = 1;
 154         }
 155       else
 156         {
 157           ndin[i] = 0;
 158         }
 159     }
 160
 161   /* Compute the initial error estimate */
 162
 163   errsum = 0.0;
 164
 165   for (i = 0; i < nint; i++)
 166     {
 167       if (ndin[i])
 168         {
 169           workspace->elist[i] = abserr0;
 170         }
 171
 172       errsum = errsum + workspace->elist[i];
 173
 174     }
 175
 176   for (i = 0; i < nint; i++)
 177     {
 178       workspace->level[i] = 0;
 179     }
 180
 181   /* Sort results into order of decreasing error via the indirection
 182      array order[] */
 183
 184   sort_results (workspace);
 185
 186   /* Test on accuracy */
 187
 188   tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
 189
 190   if (abserr0 <= 100 * GSL_DBL_EPSILON * resabs0 && abserr0 > tolerance)
 191     {
 192       *result = result0;
 193       *abserr = abserr0;
 194
 195       GSL_ERROR ("cannot reach tolerance because of roundoff error"
 196                  "on first attempt", GSL_EROUND);
 197     }
 198   else if (abserr0 <= tolerance)
 199     {
 200       *result = result0;
 201       *abserr = abserr0;
 202
 203       return GSL_SUCCESS;
 204     }
 205   else if (limit == 1)
 206     {
 207       *result = result0;
 208       *abserr = abserr0;
 209
 210       GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
 211     }
 212
 213   /* Initialization */
 214
 215   initialise_table (&table);
 216   append_table (&table, result0);
 217
 218   area = result0;
 219
 220   res_ext = result0;
 221   err_ext = GSL_DBL_MAX;
 222
 223   error_over_large_intervals = errsum;
 224   ertest = tolerance;
 225
 226   positive_integrand = test_positivity (result0, resabs0);
 227
 228   iteration = nint - 1;
 229
 230   do
 231     {
 232       size_t current_level;
 233       double a1, b1, a2, b2;
 234       double a_i, b_i, r_i, e_i;
 235       double area1 = 0, area2 = 0, area12 = 0;
 236       double error1 = 0, error2 = 0, error12 = 0;
 237       double resasc1, resasc2;
 238       double resabs1, resabs2;
 239       double last_e_i;
 240
 241       /* Bisect the subinterval with the largest error estimate */
 242
 243       retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
 244
 245       current_level = workspace->level[workspace->i] + 1;
 246
 247       a1 = a_i;
 248       b1 = 0.5 * (a_i + b_i);
 249       a2 = b1;
 250       b2 = b_i;
 251
 252       iteration++;
 253
 254       q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
 255       q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
 256
 257       area12 = area1 + area2;
 258       error12 = error1 + error2;
 259       last_e_i = e_i;
 260
 261       /* Improve previous approximations to the integral and test for
 262          accuracy.
 263
 264          We write these expressions in the same way as the original
 265          QUADPACK code so that the rounding errors are the same, which
 266          makes testing easier. */
 267
 268       errsum = errsum + error12 - e_i;
 269       area = area + area12 - r_i;
 270
 271       tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
 272
 273       if (resasc1 != error1 && resasc2 != error2)
 274         {
 275           double delta = r_i - area12;
 276
 277           if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
 278             {
 279               if (!extrapolate)
 280                 {
 281                   roundoff_type1++;
 282                 }
 283               else
 284                 {
 285                   roundoff_type2++;
 286                 }
 287             }
 288
 289           if (i > 10 && error12 > e_i)
 290             {
 291               roundoff_type3++;
 292             }
 293         }
 294
 295       /* Test for roundoff and eventually set error flag */
 296
 297       if (roundoff_type1 + roundoff_type2 >= 10 || roundoff_type3 >= 20)
 298         {
 299           error_type = 2;       /* round off error */
 300         }
 301
 302       if (roundoff_type2 >= 5)
 303         {
 304           error_type2 = 1;
 305         }
 306
 307       /* set error flag in the case of bad integrand behaviour at
 308          a point of the integration range */
 309
 310       if (subinterval_too_small (a1, a2, b2))
 311         {
 312           error_type = 4;
 313         }
 314
 315       /* append the newly-created intervals to the list */
 316
 317       update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
 318
 319       if (errsum <= tolerance)
 320         {
 321           goto compute_result;
 322         }
 323
 324       if (error_type)
 325         {
 326           break;
 327         }
 328
 329       if (iteration >= limit - 1)
 330         {
 331           error_type = 1;
 332           break;
 333         }
 334
 335       if (disallow_extrapolation)
 336         {
 337           continue;
 338         }
 339
 340       error_over_large_intervals += -last_e_i;
 341
 342       if (current_level < workspace->maximum_level)
 343         {
 344           error_over_large_intervals += error12;
 345         }
 346
 347       if (!extrapolate)
 348         {
 349           /* test whether the interval to be bisected next is the
 350              smallest interval. */
 351           if (large_interval (workspace))
 352             continue;
 353
 354           extrapolate = 1;
 355           workspace->nrmax = 1;
 356         }
 357
 358       /* The smallest interval has the largest error.  Before
 359          bisecting decrease the sum of the errors over the larger
 360          intervals (error_over_large_intervals) and perform
 361          extrapolation. */
 362
 363       if (!error_type2 && error_over_large_intervals > ertest)
 364         {
 365           if (increase_nrmax (workspace))
 366             continue;
 367         }
 368
 369       /* Perform extrapolation */
 370
 371       append_table (&table, area);
 372
 373       if (table.n < 3)
 374         {
 375           goto skip_extrapolation;
 376         }
 377
 378       qelg (&table, &reseps, &abseps);
 379
 380       ktmin++;
 381
 382       if (ktmin > 5 && err_ext < 0.001 * errsum)
 383         {
 384           error_type = 5;
 385         }
 386
 387       if (abseps < err_ext)
 388         {
 389           ktmin = 0;
 390           err_ext = abseps;
 391           res_ext = reseps;
 392           correc = error_over_large_intervals;
 393           ertest = GSL_MAX_DBL (epsabs, epsrel * fabs (reseps));
 394           if (err_ext <= ertest)
 395             break;
 396         }
 397
 398       /* Prepare bisection of the smallest interval. */
 399
 400       if (table.n == 1)
 401         {
 402           disallow_extrapolation = 1;
 403         }
 404
 405       if (error_type == 5)
 406         {
 407           break;
 408         }
 409
 410     skip_extrapolation:
 411
 412       reset_nrmax (workspace);
 413       extrapolate = 0;
 414       error_over_large_intervals = errsum;
 415
 416     }
 417   while (iteration < limit);
 418
 419   *result = res_ext;
 420   *abserr = err_ext;
 421
 422   if (err_ext == GSL_DBL_MAX)
 423     goto compute_result;
 424
 425   if (error_type || error_type2)
 426     {
 427       if (error_type2)
 428         {
 429           err_ext += correc;
 430         }
 431
 432       if (error_type == 0)
 433         error_type = 3;
 434
 435       if (result != 0 && area != 0)
 436         {
 437           if (err_ext / fabs (res_ext) > errsum / fabs (area))
 438             goto compute_result;
 439         }
 440       else if (err_ext > errsum)
 441         {
 442           goto compute_result;
 443         }
 444       else if (area == 0.0)
 445         {
 446           goto return_error;
 447         }
 448     }
 449
 450   /*  Test on divergence. */
 451
 452   {
 453     double max_area = GSL_MAX_DBL (fabs (res_ext), fabs (area));
 454
 455     if (!positive_integrand && max_area < 0.01 * resabs0)
 456       goto return_error;
 457   }
 458
 459   {
 460     double ratio = res_ext / area;
 461
 462     if (ratio < 0.01 || ratio > 100 || errsum > fabs (area))
 463       error_type = 6;
 464   }
 465
 466   goto return_error;
 467
 468 compute_result:
 469
 470   *result = sum_results (workspace);
 471   *abserr = errsum;
 472
 473 return_error:
 474
 475   if (error_type > 2)
 476     error_type--;
 477
 478   if (error_type == 0)
 479     {
 480       return GSL_SUCCESS;
 481     }
 482   else if (error_type == 1)
 483     {
 484       GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
 485     }
 486   else if (error_type == 2)
 487     {
 488       GSL_ERROR ("cannot reach tolerance because of roundoff error",
 489                  GSL_EROUND);
 490     }
 491   else if (error_type == 3)
 492     {
 493       GSL_ERROR ("bad integrand behavior found in the integration interval",
 494                  GSL_ESING);
 495     }
 496   else if (error_type == 4)
 497     {
 498       GSL_ERROR ("roundoff error detected in the extrapolation table",
 499                  GSL_EROUND);
 500     }
 501   else if (error_type == 5)
 502     {
 503       GSL_ERROR ("integral is divergent, or slowly convergent",
 504                  GSL_EDIVERGE);
 505     }
 506   else
 507     {
 508       GSL_ERROR ("could not integrate function", GSL_EFAILED);
 509     }
 510 }