/* integration/fixed.c
*
* Copyright (C) 2017 Patrick Alken
*
* 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.
*/
/* the code in this module performs fixed-point quadrature calculations for
* integrands and is based on IQPACK */
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_integration.h>
static int fixed_compute(const double a, const double b, const double alpha, const double beta,
gsl_integration_fixed_workspace * w);
static int imtqlx ( const int n, double d[], double e[], double z[] );
gsl_integration_fixed_workspace *
gsl_integration_fixed_alloc(const gsl_integration_fixed_type * type, const size_t n,
const double a, const double b, const double alpha, const double beta)
{
int status;
gsl_integration_fixed_workspace *w;
/* check inputs */
if (n < 1)
{
GSL_ERROR_VAL ("workspace size n must be at least 1", GSL_EDOM, 0);
}
w = calloc(1, sizeof(gsl_integration_fixed_workspace));
if (w == NULL)
{
GSL_ERROR_VAL ("unable to allocate workspace", GSL_ENOMEM, 0);
}
w->weights = malloc(n * sizeof(double));
if (w->weights == NULL)
{
gsl_integration_fixed_free(w);
GSL_ERROR_VAL ("unable to allocate weights", GSL_ENOMEM, 0);
}
w->x = malloc(n * sizeof(double));
if (w->x == NULL)
{
gsl_integration_fixed_free(w);
GSL_ERROR_VAL ("unable to allocate x", GSL_ENOMEM, 0);
}
w->diag = malloc(n * sizeof(double));
if (w->diag == NULL)
{
gsl_integration_fixed_free(w);
GSL_ERROR_VAL ("unable to allocate diag", GSL_ENOMEM, 0);
}
w->subdiag = malloc(n * sizeof(double));
if (w->subdiag == NULL)
{
gsl_integration_fixed_free(w);
GSL_ERROR_VAL ("unable to allocate subdiag", GSL_ENOMEM, 0);
}
w->n = n;
w->type = type;
/* compute quadrature weights and nodes */
status = fixed_compute(a, b, alpha, beta, w);
if (status)
{
gsl_integration_fixed_free(w);
GSL_ERROR_VAL ("error in integration parameters", GSL_EDOM, 0);
}
return w;
}
void
gsl_integration_fixed_free(gsl_integration_fixed_workspace * w)
{
if (w->weights)
free(w->weights);
if (w->x)
free(w->x);
if (w->diag)
free(w->diag);
if (w->subdiag)
free(w->subdiag);
free(w);
}
size_t
gsl_integration_fixed_n(const gsl_integration_fixed_workspace * w)
{
return w->n;
}
double *
gsl_integration_fixed_nodes(const gsl_integration_fixed_workspace * w)
{
return w->x;
}
double *
gsl_integration_fixed_weights(const gsl_integration_fixed_workspace * w)
{
return w->weights;
}
int
gsl_integration_fixed(const gsl_function * func, double * result,
const gsl_integration_fixed_workspace * w)
{
const size_t n = w->n;
size_t i;
double sum = 0.0;
for (i = 0; i < n; ++i)
{
double fi = GSL_FN_EVAL(func, w->x[i]);
sum += w->weights[i] * fi;
}
*result = sum;
return GSL_SUCCESS;
}
/*
fixed_compute()
Compute quadrature weights and nodes
*/
static int
fixed_compute(const double a, const double b, const double alpha, const double beta,
gsl_integration_fixed_workspace * w)
{
int s;
const size_t n = w->n;
gsl_integration_fixed_params params;
size_t i;
params.a = a;
params.b = b;
params.alpha = alpha;
params.beta = beta;
/* check input parameters */
s = (w->type->check)(n, ¶ms);
if (s)
return s;
/* initialize Jacobi matrix */
s = (w->type->init)(n, w->diag, w->subdiag, ¶ms);
if (s)
return s;
if (params.zemu <= 0.0)
{
GSL_ERROR("zeroth moment must be positive", GSL_EINVAL);
}
for ( i = 0; i < n; i++ )
{
w->x[i] = w->diag[i];
}
w->weights[0] = sqrt (params.zemu);
for ( i = 1; i < n; i++ )
{
w->weights[i] = 0.0;
}
/* diagonalize the Jacobi matrix */
s = imtqlx (n, w->x, w->subdiag, w->weights);
if (s)
return s;
for (i = 0; i < n; i++)
{
w->weights[i] = w->weights[i] * w->weights[i];
}
/*
* The current weights and nodes are valid for a = 0, b = 1.
* Now scale them for arbitrary a,b
*/
{
double p = pow ( params.slp, params.al + params.be + 1.0 );
size_t k;
for ( k = 0; k < n; k++ )
{
w->x[k] = params.shft + params.slp * w->x[k];
w->weights[k] = w->weights[k] * p;
}
}
return GSL_SUCCESS;
}
/******************************************************************************/
/*
Purpose:
IMTQLX diagonalizes a symmetric tridiagonal matrix.
Discussion:
This routine is a slightly modified version of the EISPACK routine to
perform the implicit QL algorithm on a symmetric tridiagonal matrix.
The authors thank the authors of EISPACK for permission to use this
routine.
It has been modified to produce the product Q' * Z, where Z is an input
vector and Q is the orthogonal matrix diagonalizing the input matrix.
The changes consist (essentially) of applying the orthogonal transformations
directly to Z as they are generated.
Licensing:
This code is distributed under the GNU LGPL license.
Modified:
11 January 2010
Author:
Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
C version by John Burkardt.
Reference:
Sylvan Elhay, Jaroslav Kautsky,
Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
Interpolatory Quadrature,
ACM Transactions on Mathematical Software,
Volume 13, Number 4, December 1987, pages 399-415.
Roger Martin, James Wilkinson,
The Implicit QL Algorithm,
Numerische Mathematik,
Volume 12, Number 5, December 1968, pages 377-383.
Parameters:
Input, int N, the order of the matrix.
Input/output, double D(N), the diagonal entries of the matrix.
On output, the information in D has been overwritten.
Input/output, double E(N), the subdiagonal entries of the
matrix, in entries E(1) through E(N-1). On output, the information in
E has been overwritten.
Input/output, double Z(N). On input, a vector. On output,
the value of Q' * Z, where Q is the matrix that diagonalizes the
input symmetric tridiagonal matrix.
*/
static int
imtqlx ( const int n, double d[], double e[], double z[] )
{
double b;
double c;
double f;
double g;
int i;
int ii;
int itn = 30;
int j;
int k;
int l;
int m;
int mml;
double p;
double r;
double s;
if ( n == 1 )
{
return GSL_SUCCESS;
}
e[n-1] = 0.0;
for ( l = 1; l <= n; l++ )
{
j = 0;
for ( ; ; )
{
for ( m = l; m <= n; m++ )
{
if ( m == n )
{
break;
}
if ( fabs ( e[m-1] ) <= GSL_DBL_EPSILON * ( fabs ( d[m-1] ) + fabs ( d[m] ) ) )
{
break;
}
}
p = d[l-1];
if ( m == l )
{
break;
}
if ( itn <= j )
{
return GSL_EMAXITER;
}
j = j + 1;
g = ( d[l] - p ) / ( 2.0 * e[l-1] );
r = sqrt ( g * g + 1.0 );
g = d[m-1] - p + e[l-1] / ( g + fabs ( r ) * GSL_SIGN ( g ) );
s = 1.0;
c = 1.0;
p = 0.0;
mml = m - l;
for ( ii = 1; ii <= mml; ii++ )
{
i = m - ii;
f = s * e[i-1];
b = c * e[i-1];
if ( fabs ( g ) <= fabs ( f ) )
{
c = g / f;
r = sqrt ( c * c + 1.0 );
e[i] = f * r;
s = 1.0 / r;
c = c * s;
}
else
{
s = f / g;
r = sqrt ( s * s + 1.0 );
e[i] = g * r;
c = 1.0 / r;
s = s * c;
}
g = d[i] - p;
r = ( d[i-1] - g ) * s + 2.0 * c * b;
p = s * r;
d[i] = g + p;
g = c * r - b;
f = z[i];
z[i] = s * z[i-1] + c * f;
z[i-1] = c * z[i-1] - s * f;
}
d[l-1] = d[l-1] - p;
e[l-1] = g;
e[m-1] = 0.0;
}
}
/*
Sorting.
*/
for ( ii = 2; ii <= m; ii++ )
{
i = ii - 1;
k = i;
p = d[i-1];
for ( j = ii; j <= n; j++ )
{
if ( d[j-1] < p )
{
k = j;
p = d[j-1];
}
}
if ( k != i )
{
d[k-1] = d[i-1];
d[i-1] = p;
p = z[i-1];
z[i-1] = z[k-1];
z[k-1] = p;
}
}
return GSL_SUCCESS;
}