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

   1 /* integration/qag.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_math.h>
  23 #include <gsl/gsl_errno.h>
  24 #include <gsl/gsl_integration.h>
  25
  26 #include "initialise.c"
  27 #include "set_initial.c"
  28 #include "qpsrt.c"
  29 #include "util.c"
  30
  31 static int
  32 qag (const gsl_function *f,
  33      const double a, const double b,
  34      const double epsabs, const double epsrel,
  35      const size_t limit,
  36      gsl_integration_workspace * workspace,
  37      double * result, double * abserr,
  38      gsl_integration_rule * q) ;
  39
  40 int
  41 gsl_integration_qag (const gsl_function *f,
  42                      double a, double b,
  43                      double epsabs, double epsrel, size_t limit,
  44                      int key,
  45                      gsl_integration_workspace * workspace,
  46                      double * result, double * abserr)
  47 {
  48   int status ;
  49   gsl_integration_rule * integration_rule = gsl_integration_qk15 ;
  50
  51   if (key < GSL_INTEG_GAUSS15)
  52     {
  53       key = GSL_INTEG_GAUSS15 ;
  54     }
  55   else if (key > GSL_INTEG_GAUSS61)
  56     {
  57       key = GSL_INTEG_GAUSS61 ;
  58     }
  59
  60   switch (key)
  61     {
  62     case GSL_INTEG_GAUSS15:
  63       integration_rule = gsl_integration_qk15 ;
  64       break ;
  65     case GSL_INTEG_GAUSS21:
  66       integration_rule = gsl_integration_qk21 ;
  67       break ;
  68     case GSL_INTEG_GAUSS31:
  69       integration_rule = gsl_integration_qk31 ;
  70       break ;
  71     case GSL_INTEG_GAUSS41:
  72       integration_rule = gsl_integration_qk41 ;
  73       break ;
  74     case GSL_INTEG_GAUSS51:
  75       integration_rule = gsl_integration_qk51 ;
  76       break ;
  77     case GSL_INTEG_GAUSS61:
  78       integration_rule = gsl_integration_qk61 ;
  79       break ;
  80     default:
  81       GSL_ERROR("value of key does specify a known integration rule",
  82                 GSL_EINVAL) ;
  83     }
  84
  85   status = qag (f, a, b, epsabs, epsrel, limit,
  86                 workspace,
  87                 result, abserr,
  88                 integration_rule) ;
  89
  90   return status ;
  91 }
  92
  93 static int
  94 qag (const gsl_function * f,
  95      const double a, const double b,
  96      const double epsabs, const double epsrel,
  97      const size_t limit,
  98      gsl_integration_workspace * workspace,
  99      double *result, double *abserr,
 100      gsl_integration_rule * q)
 101 {
 102   double area, errsum;
 103   double result0, abserr0, resabs0, resasc0;
 104   double tolerance;
 105   size_t iteration = 0;
 106   int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;
 107
 108   double round_off;
 109
 110   /* Initialize results */
 111
 112   initialise (workspace, a, b);
 113
 114   *result = 0;
 115   *abserr = 0;
 116
 117   if (limit > workspace->limit)
 118     {
 119       GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
 120     }
 121
 122   if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))
 123     {
 124       GSL_ERROR ("tolerance cannot be acheived with given epsabs and epsrel",
 125                  GSL_EBADTOL);
 126     }
 127
 128   /* perform the first integration */
 129
 130   q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);
 131
 132   set_initial_result (workspace, result0, abserr0);
 133
 134   /* Test on accuracy */
 135
 136   tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));
 137
 138   /* need IEEE rounding here to match original quadpack behavior */
 139
 140   round_off = GSL_COERCE_DBL (50 * GSL_DBL_EPSILON * resabs0);
 141
 142   if (abserr0 <= round_off && abserr0 > tolerance)
 143     {
 144       *result = result0;
 145       *abserr = abserr0;
 146
 147       GSL_ERROR ("cannot reach tolerance because of roundoff error "
 148                  "on first attempt", GSL_EROUND);
 149     }
 150   else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)
 151     {
 152       *result = result0;
 153       *abserr = abserr0;
 154
 155       return GSL_SUCCESS;
 156     }
 157   else if (limit == 1)
 158     {
 159       *result = result0;
 160       *abserr = abserr0;
 161
 162       GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);
 163     }
 164
 165   area = result0;
 166   errsum = abserr0;
 167
 168   iteration = 1;
 169
 170   do
 171     {
 172       double a1, b1, a2, b2;
 173       double a_i, b_i, r_i, e_i;
 174       double area1 = 0, area2 = 0, area12 = 0;
 175       double error1 = 0, error2 = 0, error12 = 0;
 176       double resasc1, resasc2;
 177       double resabs1, resabs2;
 178
 179       /* Bisect the subinterval with the largest error estimate */
 180
 181       retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
 182
 183       a1 = a_i;
 184       b1 = 0.5 * (a_i + b_i);
 185       a2 = b1;
 186       b2 = b_i;
 187
 188       q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);
 189       q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);
 190
 191       area12 = area1 + area2;
 192       error12 = error1 + error2;
 193
 194       errsum += (error12 - e_i);
 195       area += area12 - r_i;
 196
 197       if (resasc1 != error1 && resasc2 != error2)
 198         {
 199           double delta = r_i - area12;
 200
 201           if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)
 202             {
 203               roundoff_type1++;
 204             }
 205           if (iteration >= 10 && error12 > e_i)
 206             {
 207               roundoff_type2++;
 208             }
 209         }
 210
 211       tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));
 212
 213       if (errsum > tolerance)
 214         {
 215           if (roundoff_type1 >= 6 || roundoff_type2 >= 20)
 216             {
 217               error_type = 2;   /* round off error */
 218             }
 219
 220           /* set error flag in the case of bad integrand behaviour at
 221              a point of the integration range */
 222
 223           if (subinterval_too_small (a1, a2, b2))
 224             {
 225               error_type = 3;
 226             }
 227         }
 228
 229       update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);
 230
 231       retrieve (workspace, &a_i, &b_i, &r_i, &e_i);
 232
 233       iteration++;
 234
 235     }
 236   while (iteration < limit && !error_type && errsum > tolerance);
 237
 238   *result = sum_results (workspace);
 239   *abserr = errsum;
 240
 241   if (errsum <= tolerance)
 242     {
 243       return GSL_SUCCESS;
 244     }
 245   else if (error_type == 2)
 246     {
 247       GSL_ERROR ("roundoff error prevents tolerance from being achieved",
 248                  GSL_EROUND);
 249     }
 250   else if (error_type == 3)
 251     {
 252       GSL_ERROR ("bad integrand behavior found in the integration interval",
 253                  GSL_ESING);
 254     }
 255   else if (iteration == limit)
 256     {
 257       GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);
 258     }
 259   else
 260     {
 261       GSL_ERROR ("could not integrate function", GSL_EFAILED);
 262     }
 263 }
 264
 265