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

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