|
Packit |
67cb25 |
/* integration/cquad.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 2010 Pedro Gonnet
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_errno.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_integration.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include "cquad_const.c"
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Allocates a workspace for the given maximum number of intervals.
|
|
Packit |
67cb25 |
Note that if the workspace gets filled, the intervals with the
|
|
Packit |
67cb25 |
lowest error estimates are dropped. The maximum number of
|
|
Packit |
67cb25 |
intervals is therefore not the maximum number of intervals
|
|
Packit |
67cb25 |
that will be computed, but merely the size of the buffer.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_integration_cquad_workspace *
|
|
Packit |
67cb25 |
gsl_integration_cquad_workspace_alloc (const size_t n)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_integration_cquad_workspace *w;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check inputs */
|
|
Packit |
67cb25 |
if (n < 3)
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("workspace size n must be at least 3", GSL_EDOM, 0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Allocate first the workspace struct */
|
|
Packit |
67cb25 |
if ((w =
|
|
Packit |
67cb25 |
(gsl_integration_cquad_workspace *)
|
|
Packit |
67cb25 |
malloc (sizeof (gsl_integration_cquad_workspace))) == NULL)
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for workspace struct",
|
|
Packit |
67cb25 |
GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Allocate the intervals */
|
|
Packit |
67cb25 |
if ((w->ivals =
|
|
Packit |
67cb25 |
(gsl_integration_cquad_ival *)
|
|
Packit |
67cb25 |
malloc (sizeof (gsl_integration_cquad_ival) * n)) == NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (w);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for the intervals", GSL_ENOMEM,
|
|
Packit |
67cb25 |
0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Allocate the max-heap indices */
|
|
Packit |
67cb25 |
if ((w->heap = (size_t *) malloc (sizeof (size_t) * n)) == NULL)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
free (w->ivals);
|
|
Packit |
67cb25 |
free (w);
|
|
Packit |
67cb25 |
GSL_ERROR_VAL ("failed to allocate space for the heap", GSL_ENOMEM, 0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Remember the size of the workspace */
|
|
Packit |
67cb25 |
w->size = n;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Return the result */
|
|
Packit |
67cb25 |
return w;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Liberates the workspace memory.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
gsl_integration_cquad_workspace_free (gsl_integration_cquad_workspace * w)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Nothing to be done? */
|
|
Packit |
67cb25 |
if (w == NULL)
|
|
Packit |
67cb25 |
return;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Free the intervals first */
|
|
Packit |
67cb25 |
if (w->ivals != NULL)
|
|
Packit |
67cb25 |
free (w->ivals);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Free the heap */
|
|
Packit |
67cb25 |
if (w->heap != NULL)
|
|
Packit |
67cb25 |
free (w->heap);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Free the structure */
|
|
Packit |
67cb25 |
free (w);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the product of the fx with one of the inverse
|
|
Packit |
67cb25 |
Vandermonde-like matrices. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
Vinvfx (const double *fx, double *c, const int d)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
switch (d)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
case 0:
|
|
Packit |
67cb25 |
for (i = 0; i <= 4; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
c[i] = 0.0;
|
|
Packit |
67cb25 |
for (j = 0; j <= 4; j++)
|
|
Packit |
67cb25 |
c[i] += V1inv[i * 5 + j] * fx[j * 8];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
case 1:
|
|
Packit |
67cb25 |
for (i = 0; i <= 8; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
c[i] = 0.0;
|
|
Packit |
67cb25 |
for (j = 0; j <= 8; j++)
|
|
Packit |
67cb25 |
c[i] += V2inv[i * 9 + j] * fx[j * 4];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
case 2:
|
|
Packit |
67cb25 |
for (i = 0; i <= 16; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
c[i] = 0.0;
|
|
Packit |
67cb25 |
for (j = 0; j <= 16; j++)
|
|
Packit |
67cb25 |
c[i] += V3inv[i * 17 + j] * fx[j * 2];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
case 3:
|
|
Packit |
67cb25 |
for (i = 0; i <= 32; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
c[i] = 0.0;
|
|
Packit |
67cb25 |
for (j = 0; j <= 32; j++)
|
|
Packit |
67cb25 |
c[i] += V4inv[i * 33 + j] * fx[j];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Downdate the interpolation given by the n coefficients c
|
|
Packit |
67cb25 |
by removing the nodes with indices in nans. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static void
|
|
Packit |
67cb25 |
downdate (double *c, int n, int d, int *nans, int nnans)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
static const int bidx[4] = { 0, 6, 16, 34 };
|
|
Packit |
67cb25 |
double b_new[34], alpha;
|
|
Packit |
67cb25 |
int i, j;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
for (i = 0; i <= n + 1; i++)
|
|
Packit |
67cb25 |
b_new[i] = bee[bidx[d] + i];
|
|
Packit |
67cb25 |
for (i = 0; i < nnans; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
b_new[n + 1] = b_new[n + 1] / Lalpha[n];
|
|
Packit |
67cb25 |
b_new[n] = (b_new[n] + xi[nans[i]] * b_new[n + 1]) / Lalpha[n - 1];
|
|
Packit |
67cb25 |
for (j = n - 1; j > 0; j--)
|
|
Packit |
67cb25 |
b_new[j] =
|
|
Packit |
67cb25 |
(b_new[j] + xi[nans[i]] * b_new[j + 1] -
|
|
Packit |
67cb25 |
Lgamma[j + 1] * b_new[j + 2]) / Lalpha[j - 1];
|
|
Packit |
67cb25 |
for (j = 0; j <= n; j++)
|
|
Packit |
67cb25 |
b_new[j] = b_new[j + 1];
|
|
Packit |
67cb25 |
alpha = c[n] / b_new[n];
|
|
Packit |
67cb25 |
for (j = 0; j < n; j++)
|
|
Packit |
67cb25 |
c[j] -= alpha * b_new[j];
|
|
Packit |
67cb25 |
c[n] = 0;
|
|
Packit |
67cb25 |
n--;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The actual integration routine.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
gsl_integration_cquad (const gsl_function * f, double a, double b,
|
|
Packit |
67cb25 |
double epsabs, double epsrel,
|
|
Packit |
67cb25 |
gsl_integration_cquad_workspace * ws,
|
|
Packit |
67cb25 |
double *result, double *abserr, size_t * nevals)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Some constants that we will need. */
|
|
Packit |
67cb25 |
static const int n[4] = { 4, 8, 16, 32 };
|
|
Packit |
67cb25 |
static const int skip[4] = { 8, 4, 2, 1 };
|
|
Packit |
67cb25 |
static const int idx[4] = { 0, 5, 14, 31 };
|
|
Packit |
67cb25 |
static const double w = M_SQRT2 / 2;
|
|
Packit |
67cb25 |
static const int ndiv_max = 20;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Actual variables (as opposed to constants above). */
|
|
Packit |
67cb25 |
double m, h, temp;
|
|
Packit |
67cb25 |
double igral, err, igral_final, err_final, err_excess;
|
|
Packit |
67cb25 |
int nivals, neval = 0;
|
|
Packit |
67cb25 |
int i, j, d, split, t;
|
|
Packit |
67cb25 |
int nnans, nans[32];
|
|
Packit |
67cb25 |
gsl_integration_cquad_ival *iv, *ivl, *ivr;
|
|
Packit |
67cb25 |
double nc, ncdiff;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check the input arguments. */
|
|
Packit |
67cb25 |
if (f == NULL)
|
|
Packit |
67cb25 |
GSL_ERROR ("function pointer shouldn't be NULL", GSL_EINVAL);
|
|
Packit |
67cb25 |
if (result == NULL)
|
|
Packit |
67cb25 |
GSL_ERROR ("result pointer shouldn't be NULL", GSL_EINVAL);
|
|
Packit |
67cb25 |
if (ws == NULL)
|
|
Packit |
67cb25 |
GSL_ERROR ("workspace pointer shouldn't be NULL", GSL_EINVAL);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check for unreasonable accuracy demands */
|
|
Packit |
67cb25 |
if (epsabs < 0.0 || epsrel < 0.0)
|
|
Packit |
67cb25 |
GSL_ERROR ("tolerances may not be negative", GSL_EBADTOL);
|
|
Packit |
67cb25 |
if (epsabs <= 0 && epsrel < GSL_DBL_EPSILON)
|
|
Packit |
67cb25 |
GSL_ERROR ("unreasonable accuracy requirement", GSL_EBADTOL);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Create the first interval. */
|
|
Packit |
67cb25 |
iv = &(ws->ivals[0]);
|
|
Packit |
67cb25 |
m = (a + b) / 2;
|
|
Packit |
67cb25 |
h = (b - a) / 2;
|
|
Packit |
67cb25 |
nnans = 0;
|
|
Packit |
67cb25 |
for (i = 0; i <= n[3]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
iv->fx[i] = GSL_FN_EVAL (f, m + xi[i] * h);
|
|
Packit |
67cb25 |
neval++;
|
|
Packit |
67cb25 |
if (!gsl_finite (iv->fx[i]))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
nans[nnans++] = i;
|
|
Packit |
67cb25 |
iv->fx[i] = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
Vinvfx (iv->fx, &(iv->c[idx[0]]), 0);
|
|
Packit |
67cb25 |
Vinvfx (iv->fx, &(iv->c[idx[3]]), 3);
|
|
Packit |
67cb25 |
Vinvfx (iv->fx, &(iv->c[idx[2]]), 2);
|
|
Packit |
67cb25 |
for (i = 0; i < nnans; i++)
|
|
Packit |
67cb25 |
iv->fx[nans[i]] = GSL_NAN;
|
|
Packit |
67cb25 |
iv->a = a;
|
|
Packit |
67cb25 |
iv->b = b;
|
|
Packit |
67cb25 |
iv->depth = 3;
|
|
Packit |
67cb25 |
iv->rdepth = 1;
|
|
Packit |
67cb25 |
iv->ndiv = 0;
|
|
Packit |
67cb25 |
iv->igral = 2 * h * iv->c[idx[3]] * w;
|
|
Packit |
67cb25 |
nc = 0.0;
|
|
Packit |
67cb25 |
for (i = n[2] + 1; i <= n[3]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = iv->c[idx[3] + i];
|
|
Packit |
67cb25 |
nc += temp * temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ncdiff = nc;
|
|
Packit |
67cb25 |
for (i = 0; i <= n[2]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = iv->c[idx[2] + i] - iv->c[idx[3] + i];
|
|
Packit |
67cb25 |
ncdiff += temp * temp;
|
|
Packit |
67cb25 |
nc += iv->c[idx[3] + i] * iv->c[idx[3] + i];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ncdiff = sqrt (ncdiff);
|
|
Packit |
67cb25 |
nc = sqrt (nc);
|
|
Packit |
67cb25 |
iv->err = ncdiff * 2 * h;
|
|
Packit |
67cb25 |
if (ncdiff / nc > 0.1 && iv->err < 2 * h * nc)
|
|
Packit |
67cb25 |
iv->err = 2 * h * nc;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialize the heaps. */
|
|
Packit |
67cb25 |
for (i = 0; i < ws->size; i++)
|
|
Packit |
67cb25 |
ws->heap[i] = i;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Initialize some global values. */
|
|
Packit |
67cb25 |
igral = iv->igral;
|
|
Packit |
67cb25 |
err = iv->err;
|
|
Packit |
67cb25 |
nivals = 1;
|
|
Packit |
67cb25 |
igral_final = 0.0;
|
|
Packit |
67cb25 |
err_final = 0.0;
|
|
Packit |
67cb25 |
err_excess = 0.0;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Main loop. */
|
|
Packit |
67cb25 |
while (nivals > 0 && err > 0.0 &&
|
|
Packit |
67cb25 |
!(err <= fabs (igral) * epsrel || err <= epsabs)
|
|
Packit |
67cb25 |
&& !(err_final > fabs (igral) * epsrel
|
|
Packit |
67cb25 |
&& err - err_final < fabs (igral) * epsrel)
|
|
Packit |
67cb25 |
&& !(err_final > epsabs && err - err_final < epsabs))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Put our finger on the interval with the largest error. */
|
|
Packit |
67cb25 |
iv = &(ws->ivals[ws->heap[0]]);
|
|
Packit |
67cb25 |
m = (iv->a + iv->b) / 2;
|
|
Packit |
67cb25 |
h = (iv->b - iv->a) / 2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* printf
|
|
Packit |
67cb25 |
("cquad: processing ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
|
|
Packit |
67cb25 |
ws->heap[0], nivals, iv->a, iv->b, iv->igral, iv->err, iv->depth);
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
/* Should we try to increase the degree? */
|
|
Packit |
67cb25 |
if (iv->depth < 3)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Keep tabs on some variables. */
|
|
Packit |
67cb25 |
d = ++iv->depth;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Get the new (missing) function values */
|
|
Packit |
67cb25 |
for (i = skip[d]; i <= 32; i += 2 * skip[d])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
iv->fx[i] = GSL_FN_EVAL (f, m + xi[i] * h);
|
|
Packit |
67cb25 |
neval++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
nnans = 0;
|
|
Packit |
67cb25 |
for (i = 0; i <= 32; i += skip[d])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (!gsl_finite (iv->fx[i]))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
nans[nnans++] = i;
|
|
Packit |
67cb25 |
iv->fx[i] = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the new coefficients. */
|
|
Packit |
67cb25 |
Vinvfx (iv->fx, &(iv->c[idx[d]]), d);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Downdate any NaNs. */
|
|
Packit |
67cb25 |
if (nnans > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
downdate (&(iv->c[idx[d]]), n[d], d, nans, nnans);
|
|
Packit |
67cb25 |
for (i = 0; i < nnans; i++)
|
|
Packit |
67cb25 |
iv->fx[nans[i]] = GSL_NAN;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the error estimate. */
|
|
Packit |
67cb25 |
nc = 0.0;
|
|
Packit |
67cb25 |
for (i = n[d - 1] + 1; i <= n[d]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = iv->c[idx[d] + i];
|
|
Packit |
67cb25 |
nc += temp * temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ncdiff = nc;
|
|
Packit |
67cb25 |
for (i = 0; i <= n[d - 1]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = iv->c[idx[d - 1] + i] - iv->c[idx[d] + i];
|
|
Packit |
67cb25 |
ncdiff += temp * temp;
|
|
Packit |
67cb25 |
nc += iv->c[idx[d] + i] * iv->c[idx[d] + i];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ncdiff = sqrt (ncdiff);
|
|
Packit |
67cb25 |
nc = sqrt (nc);
|
|
Packit |
67cb25 |
iv->err = ncdiff * 2 * h;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the local integral. */
|
|
Packit |
67cb25 |
iv->igral = 2 * h * w * iv->c[idx[d]];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Split the interval prematurely? */
|
|
Packit |
67cb25 |
split = (nc > 0 && ncdiff / nc > 0.1);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Maximum degree reached, just split. */
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
split = 1;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Should we drop this interval? */
|
|
Packit |
67cb25 |
if ((m + h * xi[0]) >= (m + h * xi[1])
|
|
Packit |
67cb25 |
|| (m + h * xi[31]) >= (m + h * xi[32])
|
|
Packit |
67cb25 |
|| iv->err < fabs (iv->igral) * GSL_DBL_EPSILON * 10)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* printf
|
|
Packit |
67cb25 |
("cquad: dumping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
|
|
Packit |
67cb25 |
ws->heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
|
|
Packit |
67cb25 |
iv->depth);
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
/* Keep this interval's contribution */
|
|
Packit |
67cb25 |
err_final += iv->err;
|
|
Packit |
67cb25 |
igral_final += iv->igral;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Swap with the last element on the heap */
|
|
Packit |
67cb25 |
t = ws->heap[nivals - 1];
|
|
Packit |
67cb25 |
ws->heap[nivals - 1] = ws->heap[0];
|
|
Packit |
67cb25 |
ws->heap[0] = t;
|
|
Packit |
67cb25 |
nivals--;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Fix up the heap */
|
|
Packit |
67cb25 |
i = 0;
|
|
Packit |
67cb25 |
while (2 * i + 1 < nivals)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Get the kids */
|
|
Packit |
67cb25 |
j = 2 * i + 1;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* If the j+1st entry exists and is larger than the jth,
|
|
Packit |
67cb25 |
use it instead. */
|
|
Packit |
67cb25 |
if (j + 1 < nivals
|
|
Packit |
67cb25 |
&& ws->ivals[ws->heap[j + 1]].err >=
|
|
Packit |
67cb25 |
ws->ivals[ws->heap[j]].err)
|
|
Packit |
67cb25 |
j++;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Do we need to move the ith entry up? */
|
|
Packit |
67cb25 |
if (ws->ivals[ws->heap[j]].err <= ws->ivals[ws->heap[i]].err)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
t = ws->heap[j];
|
|
Packit |
67cb25 |
ws->heap[j] = ws->heap[i];
|
|
Packit |
67cb25 |
ws->heap[i] = t;
|
|
Packit |
67cb25 |
i = j;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Do we need to split this interval? */
|
|
Packit |
67cb25 |
else if (split)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Some values we will need often... */
|
|
Packit |
67cb25 |
d = iv->depth;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Generate the interval on the left */
|
|
Packit |
67cb25 |
ivl = &(ws->ivals[ws->heap[nivals++]]);
|
|
Packit |
67cb25 |
ivl->a = iv->a;
|
|
Packit |
67cb25 |
ivl->b = m;
|
|
Packit |
67cb25 |
ivl->depth = 0;
|
|
Packit |
67cb25 |
ivl->rdepth = iv->rdepth + 1;
|
|
Packit |
67cb25 |
ivl->fx[0] = iv->fx[0];
|
|
Packit |
67cb25 |
ivl->fx[32] = iv->fx[16];
|
|
Packit |
67cb25 |
for (i = skip[0]; i < 32; i += skip[0])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ivl->fx[i] =
|
|
Packit |
67cb25 |
GSL_FN_EVAL (f, (ivl->a + ivl->b) / 2 + xi[i] * h / 2);
|
|
Packit |
67cb25 |
neval++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
nnans = 0;
|
|
Packit |
67cb25 |
for (i = 0; i <= 32; i += skip[0])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (!gsl_finite (ivl->fx[i]))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
nans[nnans++] = i;
|
|
Packit |
67cb25 |
ivl->fx[i] = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
Vinvfx (ivl->fx, ivl->c, 0);
|
|
Packit |
67cb25 |
if (nnans > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
downdate (ivl->c, n[0], 0, nans, nnans);
|
|
Packit |
67cb25 |
for (i = 0; i < nnans; i++)
|
|
Packit |
67cb25 |
ivl->fx[nans[i]] = GSL_NAN;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
for (i = 0; i <= n[d]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ivl->c[idx[d] + i] = 0.0;
|
|
Packit |
67cb25 |
for (j = i; j <= n[d]; j++)
|
|
Packit |
67cb25 |
ivl->c[idx[d] + i] += Tleft[i * 33 + j] * iv->c[idx[d] + j];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ncdiff = 0.0;
|
|
Packit |
67cb25 |
for (i = 0; i <= n[0]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = ivl->c[i] - ivl->c[idx[d] + i];
|
|
Packit |
67cb25 |
ncdiff += temp * temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
for (i = n[0] + 1; i <= n[d]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = ivl->c[idx[d] + i];
|
|
Packit |
67cb25 |
ncdiff += temp * temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ncdiff = sqrt (ncdiff);
|
|
Packit |
67cb25 |
ivl->err = ncdiff * h;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check for divergence. */
|
|
Packit |
67cb25 |
ivl->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
|
|
Packit |
67cb25 |
&& ivl->c[0] / iv->c[0] > 2);
|
|
Packit |
67cb25 |
if (ivl->ndiv > ndiv_max && 2 * ivl->ndiv > ivl->rdepth)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* need copysign(INFINITY, igral) */
|
|
Packit |
67cb25 |
*result = (igral >= 0) ? GSL_POSINF : GSL_NEGINF;
|
|
Packit |
67cb25 |
if (nevals != NULL)
|
|
Packit |
67cb25 |
*nevals = neval;
|
|
Packit |
67cb25 |
return GSL_EDIVERGE;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the local integral. */
|
|
Packit |
67cb25 |
ivl->igral = h * w * ivl->c[0];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Generate the interval on the right */
|
|
Packit |
67cb25 |
ivr = &(ws->ivals[ws->heap[nivals++]]);
|
|
Packit |
67cb25 |
ivr->a = m;
|
|
Packit |
67cb25 |
ivr->b = iv->b;
|
|
Packit |
67cb25 |
ivr->depth = 0;
|
|
Packit |
67cb25 |
ivr->rdepth = iv->rdepth + 1;
|
|
Packit |
67cb25 |
ivr->fx[0] = iv->fx[16];
|
|
Packit |
67cb25 |
ivr->fx[32] = iv->fx[32];
|
|
Packit |
67cb25 |
for (i = skip[0]; i < 32; i += skip[0])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ivr->fx[i] =
|
|
Packit |
67cb25 |
GSL_FN_EVAL (f, (ivr->a + ivr->b) / 2 + xi[i] * h / 2);
|
|
Packit |
67cb25 |
neval++;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
nnans = 0;
|
|
Packit |
67cb25 |
for (i = 0; i <= 32; i += skip[0])
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
if (!gsl_finite (ivr->fx[i]))
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
nans[nnans++] = i;
|
|
Packit |
67cb25 |
ivr->fx[i] = 0.0;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
Vinvfx (ivr->fx, ivr->c, 0);
|
|
Packit |
67cb25 |
if (nnans > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
downdate (ivr->c, n[0], 0, nans, nnans);
|
|
Packit |
67cb25 |
for (i = 0; i < nnans; i++)
|
|
Packit |
67cb25 |
ivr->fx[nans[i]] = GSL_NAN;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
for (i = 0; i <= n[d]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
ivr->c[idx[d] + i] = 0.0;
|
|
Packit |
67cb25 |
for (j = i; j <= n[d]; j++)
|
|
Packit |
67cb25 |
ivr->c[idx[d] + i] += Tright[i * 33 + j] * iv->c[idx[d] + j];
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ncdiff = 0.0;
|
|
Packit |
67cb25 |
for (i = 0; i <= n[0]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = ivr->c[i] - ivr->c[idx[d] + i];
|
|
Packit |
67cb25 |
ncdiff += temp * temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
for (i = n[0] + 1; i <= n[d]; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
temp = ivr->c[idx[d] + i];
|
|
Packit |
67cb25 |
ncdiff += temp * temp;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
ncdiff = sqrt (ncdiff);
|
|
Packit |
67cb25 |
ivr->err = ncdiff * h;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Check for divergence. */
|
|
Packit |
67cb25 |
ivr->ndiv = iv->ndiv + (fabs (iv->c[0]) > 0
|
|
Packit |
67cb25 |
&& ivr->c[0] / iv->c[0] > 2);
|
|
Packit |
67cb25 |
if (ivr->ndiv > ndiv_max && 2 * ivr->ndiv > ivr->rdepth)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* need copysign(INFINITY, igral) */
|
|
Packit |
67cb25 |
*result = (igral >= 0) ? GSL_POSINF : GSL_NEGINF;
|
|
Packit |
67cb25 |
if (nevals != NULL)
|
|
Packit |
67cb25 |
*nevals = neval;
|
|
Packit |
67cb25 |
return GSL_EDIVERGE;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute the local integral. */
|
|
Packit |
67cb25 |
ivr->igral = h * w * ivr->c[0];
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Fix-up the heap: we now have one interval on top
|
|
Packit |
67cb25 |
that we don't need any more and two new, unsorted
|
|
Packit |
67cb25 |
ones at the bottom. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Flip the last interval to the top of the heap and
|
|
Packit |
67cb25 |
sift down. */
|
|
Packit |
67cb25 |
t = ws->heap[nivals - 1];
|
|
Packit |
67cb25 |
ws->heap[nivals - 1] = ws->heap[0];
|
|
Packit |
67cb25 |
ws->heap[0] = t;
|
|
Packit |
67cb25 |
nivals--;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Sift this interval back down the heap. */
|
|
Packit |
67cb25 |
i = 0;
|
|
Packit |
67cb25 |
while (2 * i + 1 < nivals - 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
j = 2 * i + 1;
|
|
Packit |
67cb25 |
if (j + 1 < nivals - 1
|
|
Packit |
67cb25 |
&& ws->ivals[ws->heap[j + 1]].err >=
|
|
Packit |
67cb25 |
ws->ivals[ws->heap[j]].err)
|
|
Packit |
67cb25 |
j++;
|
|
Packit |
67cb25 |
if (ws->ivals[ws->heap[j]].err <= ws->ivals[ws->heap[i]].err)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
t = ws->heap[j];
|
|
Packit |
67cb25 |
ws->heap[j] = ws->heap[i];
|
|
Packit |
67cb25 |
ws->heap[i] = t;
|
|
Packit |
67cb25 |
i = j;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Now grab the last interval and sift it up the heap. */
|
|
Packit |
67cb25 |
i = nivals - 1;
|
|
Packit |
67cb25 |
while (i > 0)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
j = (i - 1) / 2;
|
|
Packit |
67cb25 |
if (ws->ivals[ws->heap[j]].err < ws->ivals[ws->heap[i]].err)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
t = ws->heap[j];
|
|
Packit |
67cb25 |
ws->heap[j] = ws->heap[i];
|
|
Packit |
67cb25 |
ws->heap[i] = t;
|
|
Packit |
67cb25 |
i = j;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Otherwise, just fix-up the heap. */
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
i = 0;
|
|
Packit |
67cb25 |
while (2 * i + 1 < nivals)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
j = 2 * i + 1;
|
|
Packit |
67cb25 |
if (j + 1 < nivals
|
|
Packit |
67cb25 |
&& ws->ivals[ws->heap[j + 1]].err >=
|
|
Packit |
67cb25 |
ws->ivals[ws->heap[j]].err)
|
|
Packit |
67cb25 |
j++;
|
|
Packit |
67cb25 |
if (ws->ivals[ws->heap[j]].err <= ws->ivals[ws->heap[i]].err)
|
|
Packit |
67cb25 |
break;
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
t = ws->heap[j];
|
|
Packit |
67cb25 |
ws->heap[j] = ws->heap[i];
|
|
Packit |
67cb25 |
ws->heap[i] = t;
|
|
Packit |
67cb25 |
i = j;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* If the heap is about to overflow, remove the last two
|
|
Packit |
67cb25 |
intervals. */
|
|
Packit |
67cb25 |
while (nivals > ws->size - 2)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
iv = &(ws->ivals[ws->heap[nivals - 1]]);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* printf
|
|
Packit |
67cb25 |
("cquad: dumping ival %i (of %i) with [%e,%e] int=%e, err=%e, depth=%i\n",
|
|
Packit |
67cb25 |
ws->heap[0], nivals, iv->a, iv->b, iv->igral, iv->err,
|
|
Packit |
67cb25 |
iv->depth);
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
err_final += iv->err;
|
|
Packit |
67cb25 |
igral_final += iv->igral;
|
|
Packit |
67cb25 |
nivals--;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Collect the value of the integral and error. */
|
|
Packit |
67cb25 |
igral = igral_final;
|
|
Packit |
67cb25 |
err = err_final;
|
|
Packit |
67cb25 |
for (i = 0; i < nivals; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
igral += ws->ivals[ws->heap[i]].igral;
|
|
Packit |
67cb25 |
err += ws->ivals[ws->heap[i]].err;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Dump the contents of the heap. */
|
|
Packit |
67cb25 |
/* for (i = 0; i < nivals; i++)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
iv = &(ws->ivals[ws->heap[i]]);
|
|
Packit |
67cb25 |
printf
|
|
Packit |
67cb25 |
("cquad: ival %i (%i) with [%e,%e], int=%e, err=%e, depth=%i, rdepth=%i\n",
|
|
Packit |
67cb25 |
i, ws->heap[i], iv->a, iv->b, iv->igral, iv->err, iv->depth,
|
|
Packit |
67cb25 |
iv->rdepth);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
/* Clean up and present the results. */
|
|
Packit |
67cb25 |
*result = igral;
|
|
Packit |
67cb25 |
if (abserr != NULL)
|
|
Packit |
67cb25 |
*abserr = err;
|
|
Packit |
67cb25 |
if (nevals != NULL)
|
|
Packit |
67cb25 |
*nevals = neval;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* All is well that ends well. */
|
|
Packit |
67cb25 |
return GSL_SUCCESS;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
}
|