/* integration/qag.c

 * 

 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Brian Gough

 * 

 * This program is free software; you can redistribute it and/or modify

 * it under the terms of the GNU General Public License as published by

 * the Free Software Foundation; either version 3 of the License, or (at

 * your option) any later version.

 * 

 * This program is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU

 * General Public License for more details.

 * 

 * You should have received a copy of the GNU General Public License

 * along with this program; if not, write to the Free Software

 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.

 */

#include <config.h>

#include <stdlib.h>

#include <gsl/gsl_math.h>

#include <gsl/gsl_errno.h>

#include <gsl/gsl_integration.h>

#include "initialise.c"

#include "set_initial.c"

#include "qpsrt.c"

#include "util.c"

static int

qag (const gsl_function *f,

     const double a, const double b,

     const double epsabs, const double epsrel,

     const size_t limit,

     gsl_integration_workspace * workspace,

     double * result, double * abserr,

     gsl_integration_rule * q) ;

int

gsl_integration_qag (const gsl_function *f,

                     double a, double b,

                     double epsabs, double epsrel, size_t limit,

                     int key,

                     gsl_integration_workspace * workspace,

                     double * result, double * abserr)

{

  int status ;

  gsl_integration_rule * integration_rule = gsl_integration_qk15 ;

  if (key < GSL_INTEG_GAUSS15)

    {

      key = GSL_INTEG_GAUSS15 ;

    } 

  else if (key > GSL_INTEG_GAUSS61) 

    {

      key = GSL_INTEG_GAUSS61 ;

    }

  switch (key) 

    {

    case GSL_INTEG_GAUSS15:

      integration_rule = gsl_integration_qk15 ;

      break ;

    case GSL_INTEG_GAUSS21:

      integration_rule = gsl_integration_qk21 ;

      break ;

    case GSL_INTEG_GAUSS31:

      integration_rule = gsl_integration_qk31 ; 

      break ;

    case GSL_INTEG_GAUSS41:

      integration_rule = gsl_integration_qk41 ;

      break ;      

    case GSL_INTEG_GAUSS51:

      integration_rule = gsl_integration_qk51 ;

      break ;      

    case GSL_INTEG_GAUSS61:

      integration_rule = gsl_integration_qk61 ;

      break ;      

    default:

      GSL_ERROR("value of key does specify a known integration rule", 

                GSL_EINVAL) ;

    }

  status = qag (f, a, b, epsabs, epsrel, limit,

                workspace, 

                result, abserr, 

                integration_rule) ;

  return status ;

}

static int

qag (const gsl_function * f,

     const double a, const double b,

     const double epsabs, const double epsrel,

     const size_t limit,

     gsl_integration_workspace * workspace,

     double *result, double *abserr,

     gsl_integration_rule * q)

{

  double area, errsum;

  double result0, abserr0, resabs0, resasc0;

  double tolerance;

  size_t iteration = 0;

  int roundoff_type1 = 0, roundoff_type2 = 0, error_type = 0;

  double round_off;     

  /* Initialize results */

  initialise (workspace, a, b);

  *result = 0;

  *abserr = 0;

  if (limit > workspace->limit)

    {

      GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;

    }

  if (epsabs <= 0 && (epsrel < 50 * GSL_DBL_EPSILON || epsrel < 0.5e-28))

    {

      GSL_ERROR ("tolerance cannot be achieved with given epsabs and epsrel",

                 GSL_EBADTOL);

    }

  /* perform the first integration */

  q (f, a, b, &result0, &abserr0, &resabs0, &resasc0);

  set_initial_result (workspace, result0, abserr0);

  /* Test on accuracy */

  tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (result0));

  /* need IEEE rounding here to match original quadpack behavior */

  round_off = GSL_COERCE_DBL (50 * GSL_DBL_EPSILON * resabs0);

  if (abserr0 <= round_off && abserr0 > tolerance)

    {

      *result = result0;

      *abserr = abserr0;

      GSL_ERROR ("cannot reach tolerance because of roundoff error "

                 "on first attempt", GSL_EROUND);

    }

  else if ((abserr0 <= tolerance && abserr0 != resasc0) || abserr0 == 0.0)

    {

      *result = result0;

      *abserr = abserr0;

      return GSL_SUCCESS;

    }

  else if (limit == 1)

    {

      *result = result0;

      *abserr = abserr0;

      GSL_ERROR ("a maximum of one iteration was insufficient", GSL_EMAXITER);

    }

  area = result0;

  errsum = abserr0;

  iteration = 1;

  do

    {

      double a1, b1, a2, b2;

      double a_i, b_i, r_i, e_i;

      double area1 = 0, area2 = 0, area12 = 0;

      double error1 = 0, error2 = 0, error12 = 0;

      double resasc1, resasc2;

      double resabs1, resabs2;

      /* Bisect the subinterval with the largest error estimate */

      retrieve (workspace, &a_i, &b_i, &r_i, &e_i);

      a1 = a_i; 

      b1 = 0.5 * (a_i + b_i);

      a2 = b1;

      b2 = b_i;

      q (f, a1, b1, &area1, &error1, &resabs1, &resasc1);

      q (f, a2, b2, &area2, &error2, &resabs2, &resasc2);

      area12 = area1 + area2;

      error12 = error1 + error2;

      errsum += (error12 - e_i);

      area += area12 - r_i;

      if (resasc1 != error1 && resasc2 != error2)

        {

          double delta = r_i - area12;

          if (fabs (delta) <= 1.0e-5 * fabs (area12) && error12 >= 0.99 * e_i)

            {

              roundoff_type1++;

            }

          if (iteration >= 10 && error12 > e_i)

            {

              roundoff_type2++;

            }

        }

      tolerance = GSL_MAX_DBL (epsabs, epsrel * fabs (area));

      if (errsum > tolerance)

        {

          if (roundoff_type1 >= 6 || roundoff_type2 >= 20)

            {

              error_type = 2;   /* round off error */

            }

          /* set error flag in the case of bad integrand behaviour at

             a point of the integration range */

          if (subinterval_too_small (a1, a2, b2))

            {

              error_type = 3;

            }

        }

      update (workspace, a1, b1, area1, error1, a2, b2, area2, error2);

      retrieve (workspace, &a_i, &b_i, &r_i, &e_i);

      iteration++;

    }

  while (iteration < limit && !error_type && errsum > tolerance);

  *result = sum_results (workspace);

  *abserr = errsum;

  if (errsum <= tolerance)

    {

      return GSL_SUCCESS;

    }

  else if (error_type == 2)

    {

      GSL_ERROR ("roundoff error prevents tolerance from being achieved",

                 GSL_EROUND);

    }

  else if (error_type == 3)

    {

      GSL_ERROR ("bad integrand behavior found in the integration interval",

                 GSL_ESING);

    }

  else if (iteration == limit)

    {

      GSL_ERROR ("maximum number of subdivisions reached", GSL_EMAXITER);

    }

  else

    {

      GSL_ERROR ("could not integrate function", GSL_EFAILED);

    }

}