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

   1 /* integration/qawf.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 "qelg.c"
  30
  31 int
  32 gsl_integration_qawf (gsl_function * f,
  33                       const double a,
  34                       const double epsabs,
  35                       const size_t limit,
  36                       gsl_integration_workspace * workspace,
  37                       gsl_integration_workspace * cycle_workspace,
  38                       gsl_integration_qawo_table * wf,
  39                       double *result, double *abserr)
  40 {
  41   double area, errsum;
  42   double res_ext, err_ext;
  43   double correc, total_error = 0.0, truncation_error;
  44
  45   size_t ktmin = 0;
  46   size_t iteration = 0;
  47
  48   struct extrapolation_table table;
  49
  50   double cycle;
  51   double omega = wf->omega;
  52
  53   const double p = 0.9;
  54   double factor = 1;
  55   double initial_eps, eps;
  56   int error_type = 0;
  57
  58   /* Initialize results */
  59
  60   initialise (workspace, a, a);
  61
  62   *result = 0;
  63   *abserr = 0;
  64
  65   if (limit > workspace->limit)
  66     {
  67       GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
  68     }
  69
  70   /* Test on accuracy */
  71
  72   if (epsabs <= 0)
  73     {
  74       GSL_ERROR ("absolute tolerance epsabs must be positive", GSL_EBADTOL) ;
  75     }
  76
  77   if (omega == 0.0)
  78     {
  79       if (wf->sine == GSL_INTEG_SINE)
  80         {
  81           /* The function sin(w x) f(x) is always zero for w = 0 */
  82
  83           *result = 0;
  84           *abserr = 0;
  85
  86           return GSL_SUCCESS;
  87         }
  88       else
  89         {
  90           /* The function cos(w x) f(x) is always f(x) for w = 0 */
  91
  92           int status = gsl_integration_qagiu (f, a, epsabs, 0.0,
  93                                               cycle_workspace->limit,
  94                                               cycle_workspace,
  95                                               result, abserr);
  96           return status;
  97         }
  98     }
  99
 100   if (epsabs > GSL_DBL_MIN / (1 - p))
 101     {
 102       eps = epsabs * (1 - p);
 103     }
 104   else
 105     {
 106       eps = epsabs;
 107     }
 108
 109   initial_eps = eps;
 110
 111   area = 0;
 112   errsum = 0;
 113
 114   res_ext = 0;
 115   err_ext = GSL_DBL_MAX;
 116   correc = 0;
 117
 118   cycle = (2 * floor (fabs (omega)) + 1) * M_PI / fabs (omega);
 119
 120   gsl_integration_qawo_table_set_length (wf, cycle);
 121
 122   initialise_table (&table);
 123
 124   for (iteration = 0; iteration < limit; iteration++)
 125     {
 126       double area1, error1, reseps, erreps;
 127
 128       double a1 = a + iteration * cycle;
 129       double b1 = a1 + cycle;
 130
 131       double epsabs1 = eps * factor;
 132
 133       int status = gsl_integration_qawo (f, a1, epsabs1, 0.0, limit,
 134                                          cycle_workspace, wf,
 135                                          &area1, &error1);
 136
 137       append_interval (workspace, a1, b1, area1, error1);
 138
 139       factor *= p;
 140
 141       area = area + area1;
 142       errsum = errsum + error1;
 143
 144       /* estimate the truncation error as 50 times the final term */
 145
 146       truncation_error = 50 * fabs (area1);
 147
 148       total_error = errsum + truncation_error;
 149
 150       if (total_error < epsabs && iteration > 4)
 151         {
 152           goto compute_result;
 153         }
 154
 155       if (error1 > correc)
 156         {
 157           correc = error1;
 158         }
 159
 160       if (status)
 161         {
 162           eps = GSL_MAX_DBL (initial_eps, correc * (1.0 - p));
 163         }
 164
 165       if (status && total_error < 10 * correc && iteration > 3)
 166         {
 167           goto compute_result;
 168         }
 169
 170       append_table (&table, area);
 171
 172       if (table.n < 2)
 173         {
 174           continue;
 175         }
 176
 177       qelg (&table, &reseps, &erreps);
 178
 179       ktmin++;
 180
 181       if (ktmin >= 15 && err_ext < 0.001 * total_error)
 182         {
 183           error_type = 4;
 184         }
 185
 186       if (erreps < err_ext)
 187         {
 188           ktmin = 0;
 189           err_ext = erreps;
 190           res_ext = reseps;
 191
 192           if (err_ext + 10 * correc <= epsabs)
 193             break;
 194           if (err_ext <= epsabs && 10 * correc >= epsabs)
 195             break;
 196         }
 197
 198     }
 199
 200   if (iteration == limit)
 201     error_type = 1;
 202
 203   if (err_ext == GSL_DBL_MAX)
 204     goto compute_result;
 205
 206   err_ext = err_ext + 10 * correc;
 207
 208   *result = res_ext;
 209   *abserr = err_ext;
 210
 211   if (error_type == 0)
 212     {
 213       return GSL_SUCCESS ;
 214     }
 215
 216   if (res_ext != 0.0 && area != 0.0)
 217     {
 218       if (err_ext / fabs (res_ext) > errsum / fabs (area))
 219         goto compute_result;
 220     }
 221   else if (err_ext > errsum)
 222     {
 223       goto compute_result;
 224     }
 225   else if (area == 0.0)
 226     {
 227       goto return_error;
 228     }
 229
 230   if (error_type == 4)
 231     {
 232       err_ext = err_ext + truncation_error;
 233     }
 234
 235   goto return_error;
 236
 237 compute_result:
 238
 239   *result = area;
 240   *abserr = total_error;
 241
 242 return_error:
 243
 244   if (error_type > 2)
 245     error_type--;
 246
 247   if (error_type == 0)
 248     {
 249       return GSL_SUCCESS;
 250     }
 251   else if (error_type == 1)
 252     {
 253       GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
 254     }
 255   else if (error_type == 2)
 256     {
 257       GSL_ERROR ("cannot reach tolerance because of roundoff error",
 258                  GSL_EROUND);
 259     }
 260   else if (error_type == 3)
 261     {
 262       GSL_ERROR ("bad integrand behavior found in the integration interval",
 263                  GSL_ESING);
 264     }
 265   else if (error_type == 4)
 266     {
 267       GSL_ERROR ("roundoff error detected in the extrapolation table",
 268                  GSL_EROUND);
 269     }
 270   else if (error_type == 5)
 271     {
 272       GSL_ERROR ("integral is divergent, or slowly convergent",
 273                  GSL_EDIVERGE);
 274     }
 275   else
 276     {
 277       GSL_ERROR ("could not integrate function", GSL_EFAILED);
 278     }
 279
 280 }
 281