Gen laguerre rule.cpp
// Gen laguerre rule.cpp is C++ routine to evaluate nodes x_n and weight w_n for n from 1 to M of the Sonin quadrature formula to approximate \(\int_0^\infty x^a \exp(-x) f(x) \mathrm d x\) by http://people.sc.fsu.edu/~jburkardt/cpp_src/gen_laguerre_rule/gen_laguerre_rule.cpp
CODE
// The modification is not allowed, so, the misprints in the comment of the code below is kept as is
//
# include <cstdlib>
# include <cmath>
# include <iostream>
# include <fstream>
# include <iomanip>
# include <ctime>
# include <cstring>
using namespace std;
int main ( int argc, char *argv[] );
void cdgqf ( int nt, int kind, double alpha, double beta, double t[],
double wts[] );
void cgqf ( int nt, int kind, double alpha, double beta, double a, double b,
double t[], double wts[] );
double class_matrix ( int kind, int m, double alpha, double beta, double aj[],
double bj[] );
void imtqlx ( int n, double d[], double e[], double z[] );
void parchk ( int kind, int m, double alpha, double beta );
double r8_epsilon ( );
double r8_huge ( );
double r8_sign ( double x );
void r8mat_write ( string output_filename, int m, int n, double table[] );
void rule_write ( int order, string filename, double x[], double w[],
double r[] );
void scqf ( int nt, double t[], int mlt[], double wts[], int nwts, int ndx[],
double swts[], double st[], int kind, double alpha, double beta, double a,
double b );
void sgqf ( int nt, double aj[], double bj[], double zemu, double t[],
double wts[] );
void timestamp ( );
//****************************************************************************80
int main ( int argc, char *argv[] )
//****************************************************************************80
//
// Purpose:
//
// MAIN is the main program for GEN_LAGUERRE_RULE.
//
// Discussion:
//
// This program computes a standard or exponentially weighted
// Gauss-Laguerre quadrature rule and writes it to a file.
//
// The user specifies:
// * the ORDER (number of points) in the rule;
// * ALPHA, the exponent of |X|;
// * A, the left endpoint;
// * B, the scale factor in the exponential;
// * FILENAME, the root name of the output files.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 23 February 2010
//
// Author:
//
// John Burkardt
//
{
double a;
double alpha;
double b;
double beta;
string filename;
int kind;
int order;
double *r;
double *w;
double *x;
cout << "\n";
timestamp ( );
cout << "\n";
cout << "GEN_LAGUERRE_RULE\n";
cout << " C++ version\n";
cout << "\n";
cout << " Compute a generalized Gauss-Laguerre rule for approximating\n";
cout << " Integral ( a <= x < +oo ) |x-a|^ALPHA exp(-B*x(x-a)) f(x) dx\n";
cout << " of order ORDER.\n";
cout << "\n";
cout << " The user specifies ORDER, ALPHA, A, B, and FILENAME.\n";
cout << "\n";
cout << " ORDER is the number of points in the rule:\n";
cout << " ALPHA is the exponent of |X|:\n";
cout << " A is the left endpoint (typically 0).\n";
cout << " B is the exponential scale factor, typically 1:\n";
cout << " FILENAME is used to generate 3 files:\n";
cout << " * filename_w.txt - the weight file\n";
cout << " * filename_x.txt - the abscissa file.\n";
cout << " * filename_r.txt - the region file.\n";
//
// Initialize parameters;
//
beta = 0.0;
//
// Get ORDER.
//
if ( 1 < argc )
{
order = atoi ( argv[1] );
}
else
{
cout << "\n";
cout << " Enter the value of ORDER (1 or greater)\n";
cin >> order;
}
//
// Get ALPHA.
//
if ( 2 < argc )
{
alpha = atof ( argv[2] );
}
else
{
cout << "\n";
cout << " Enter the value of ALPHA.\n";
cout << " ( -1.0 < ALPHA is required.)\n";
cin >> alpha;
}
//
// Get A.
//
if ( 3 < argc )
{
a = atof ( argv[3] );
}
else
{
cout << "\n";
cout << " Enter the left endpoint A.\n";
cin >> a;
}
//
// Get B.
//
if ( 4 < argc )
{
b = atof ( argv[4] );
}
else
{
cout << "\n";
cout << " Enter the value of B\n";
cin >> b;
}
//
// Get FILENAME.
//
if ( 5 < argc )
{
filename = argv[5];
}
else
{
cout << "\n";
cout << " Enter FILENAME the \"root name\" of the quadrature files).\n";
cin >> filename;
}
//
// Input summary.
//
cout << "\n";
cout << " ORDER = " << order << "\n";
cout << " ALPHA = " << alpha << "\n";
cout << " A = " << a << "\n";
cout << " B = " << b << "\n";
cout << " FILENAME \"" << filename << "\".\n";
//
// Construct the rule.
//
w = new double[order];
x = new double[order];
kind = 5;
cgqf ( order, kind, alpha, beta, a, b, x, w );
//
// Write the rule.
//
r = new double[2];
r[0] = a;
r[1] = r8_huge ( );
rule_write ( order, filename, x, w, r );
//
// Free memory.
//
delete [] r;
delete [] w;
delete [] x;
//
// Terminate.
//
cout << "\n";
cout << "GEN_LAGUERRE_RULE:\n";
cout << " Normal end of execution.\n";
cout << "\n";
timestamp ( );
return 0;
}
//****************************************************************************80
void cdgqf ( int nt, int kind, double alpha, double beta, double t[],
double wts[] )
//****************************************************************************80
//
// Purpose:
//
// CDGQF computes a Gauss quadrature formula with default A, B and simple knots.
//
// Discussion:
//
// This routine computes all the knots and weights of a Gauss quadrature
// formula with a classical weight function with default values for A and B,
// and only simple knots.
//
// There are no moments checks and no printing is done.
//
// Use routine EIQFS to evaluate a quadrature computed by CGQFS.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 08 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.
//
// Parameters:
//
// Input, int NT, the number of knots.
//
// Input, int KIND, the rule.
// 1, Legendre, (a,b) 1.0
// 2, Chebyshev, (a,b) ((b-x)*(x-a))^(-0.5)
// 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha
// 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta
// 5, Generalized Laguerre, (a,inf) (x-a)^alpha*exp(-b*(x-a))
// 6, Generalized Hermite, (-inf,inf) |x-a|^alpha*exp(-b*(x-a)^2)
// 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha
// 8, Rational, (a,inf) (x-a)^alpha*(x+b)^beta
//
// Input, double ALPHA, the value of Alpha, if needed.
//
// Input, double BETA, the value of Beta, if needed.
//
// Output, double T[NT], the knots.
//
// Output, double WTS[NT], the weights.
//
{
double *aj;
double *bj;
double zemu;
parchk ( kind, 2 * nt, alpha, beta );
//
// Get the Jacobi matrix and zero-th moment.
//
aj = new double[nt];
bj = new double[nt];
zemu = class_matrix ( kind, nt, alpha, beta, aj, bj );
//
// Compute the knots and weights.
//
sgqf ( nt, aj, bj, zemu, t, wts );
delete [] aj;
delete [] bj;
return;
}
//****************************************************************************80
void cgqf ( int nt, int kind, double alpha, double beta, double a, double b,
double t[], double wts[] )
//****************************************************************************80
//
// Purpose:
//
// CGQF computes knots and weights of a Gauss quadrature formula.
//
// Discussion:
//
// The user may specify the interval (A,B).
//
// Only simple knots are produced.
//
// Use routine EIQFS to evaluate this quadrature formula.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 16 February 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.
//
// Parameters:
//
// Input, int NT, the number of knots.
//
// Input, int KIND, the rule.
// 1, Legendre, (a,b) 1.0
// 2, Chebyshev Type 1, (a,b) ((b-x)*(x-a))^-0.5)
// 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha
// 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta
// 5, Generalized Laguerre, (a,+oo) (x-a)^alpha*exp(-b*(x-a))
// 6, Generalized Hermite, (-oo,+oo) |x-a|^alpha*exp(-b*(x-a)^2)
// 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha
// 8, Rational, (a,+oo) (x-a)^alpha*(x+b)^beta
// 9, Chebyshev Type 2, (a,b) ((b-x)*(x-a))^(+0.5)
//
// Input, double ALPHA, the value of Alpha, if needed.
//
// Input, double BETA, the value of Beta, if needed.
//
// Input, double A, B, the interval endpoints, or
// other parameters.
//
// Output, double T[NT], the knots.
//
// Output, double WTS[NT], the weights.
//
{
int i;
int *mlt;
int *ndx;
//
// Compute the Gauss quadrature formula for default values of A and B.
//
cdgqf ( nt, kind, alpha, beta, t, wts );
//
// Prepare to scale the quadrature formula to other weight function with
// valid A and B.
//
mlt = new int[nt];
for ( i = 0; i < nt; i++ )
{
mlt[i] = 1;
}
ndx = new int[nt];
for ( i = 0; i < nt; i++ )
{
ndx[i] = i + 1;
}
scqf ( nt, t, mlt, wts, nt, ndx, wts, t, kind, alpha, beta, a, b );
delete [] mlt;
delete [] ndx;
return;
}
//****************************************************************************80
double class_matrix ( int kind, int m, double alpha, double beta, double aj[],
double bj[] )
//****************************************************************************80
//
// Purpose:
//
// CLASS_MATRIX computes the Jacobi matrix for a quadrature rule.
//
// Discussion:
//
// This routine computes the diagonal AJ and sub-diagonal BJ
// elements of the order M tridiagonal symmetric Jacobi matrix
// associated with the polynomials orthogonal with respect to
// the weight function specified by KIND.
//
// For weight functions 1-7, M elements are defined in BJ even
// though only M-1 are needed. For weight function 8, BJ(M) is
// set to zero.
//
// The zero-th moment of the weight function is returned in ZEMU.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 08 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.
//
// Parameters:
//
// Input, int KIND, the rule.
// 1, Legendre, (a,b) 1.0
// 2, Chebyshev, (a,b) ((b-x)*(x-a))^(-0.5)
// 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha
// 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta
// 5, Generalized Laguerre, (a,inf) (x-a)^alpha*exp(-b*(x-a))
// 6, Generalized Hermite, (-inf,inf) |x-a|^alpha*exp(-b*(x-a)^2)
// 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha
// 8, Rational, (a,inf) (x-a)^alpha*(x+b)^beta
//
// Input, int M, the order of the Jacobi matrix.
//
// Input, double ALPHA, the value of Alpha, if needed.
//
// Input, double BETA, the value of Beta, if needed.
//
// Output, double AJ[M], BJ[M], the diagonal and subdiagonal
// of the Jacobi matrix.
//
// Output, double CLASS_MATRIX, the zero-th moment.
//
{
double a2b2;
double ab;
double aba;
double abi;
double abj;
double abti;
double apone;
int i;
double pi = 3.14159265358979323846264338327950;
double temp;
double temp2;
double zemu;
temp = r8_epsilon ( );
parchk ( kind, 2 * m - 1, alpha, beta );
temp2 = 0.5;
if ( 500.0 * temp < fabs ( pow ( tgamma ( temp2 ), 2 ) - pi ) )
{
cout << "\n";
cout << "CLASS_MATRIX - Fatal error!\n";
cout << " Gamma function does not match machine parameters.\n";
exit ( 1 );
}
if ( kind == 1 )
{
ab = 0.0;
zemu = 2.0 / ( ab + 1.0 );
for ( i = 0; i < m; i++ )
{
aj[i] = 0.0;
}
for ( i = 1; i <= m; i++ )
{
abi = i + ab * ( i % 2 );
abj = 2 * i + ab;
bj[i-1] = sqrt ( abi * abi / ( abj * abj - 1.0 ) );
}
}
else if ( kind == 2 )
{
zemu = pi;
for ( i = 0; i < m; i++ )
{
aj[i] = 0.0;
}
bj[0] = sqrt ( 0.5 );
for ( i = 1; i < m; i++ )
{
bj[i] = 0.5;
}
}
else if ( kind == 3 )
{
ab = alpha * 2.0;
zemu = pow ( 2.0, ab + 1.0 ) * pow ( tgamma ( alpha + 1.0 ), 2 )
/ tgamma ( ab + 2.0 );
for ( i = 0; i < m; i++ )
{
aj[i] = 0.0;
}
bj[0] = sqrt ( 1.0 / ( 2.0 * alpha + 3.0 ) );
for ( i = 2; i <= m; i++ )
{
bj[i-1] = sqrt ( i * ( i + ab ) / ( 4.0 * pow ( i + alpha, 2 ) - 1.0 ) );
}
}
else if ( kind == 4 )
{
ab = alpha + beta;
abi = 2.0 + ab;
zemu = pow ( 2.0, ab + 1.0 ) * tgamma ( alpha + 1.0 )
* tgamma ( beta + 1.0 ) / tgamma ( abi );
aj[0] = ( beta - alpha ) / abi;
bj[0] = sqrt ( 4.0 * ( 1.0 + alpha ) * ( 1.0 + beta )
/ ( ( abi + 1.0 ) * abi * abi ) );
a2b2 = beta * beta - alpha * alpha;
for ( i = 2; i <= m; i++ )
{
abi = 2.0 * i + ab;
aj[i-1] = a2b2 / ( ( abi - 2.0 ) * abi );
abi = abi * abi;
bj[i-1] = sqrt ( 4.0 * i * ( i + alpha ) * ( i + beta ) * ( i + ab )
/ ( ( abi - 1.0 ) * abi ) );
}
}
else if ( kind == 5 )
{
zemu = tgamma ( alpha + 1.0 );
for ( i = 1; i <= m; i++ )
{
aj[i-1] = 2.0 * i - 1.0 + alpha;
bj[i-1] = sqrt ( i * ( i + alpha ) );
}
}
else if ( kind == 6 )
{
zemu = tgamma ( ( alpha + 1.0 ) / 2.0 );
for ( i = 0; i < m; i++ )
{
aj[i] = 0.0;
}
for ( i = 1; i <= m; i++ )
{
bj[i-1] = sqrt ( ( i + alpha * ( i % 2 ) ) / 2.0 );
}
}
else if ( kind == 7 )
{
ab = alpha;
zemu = 2.0 / ( ab + 1.0 );
for ( i = 0; i < m; i++ )
{
aj[i] = 0.0;
}
for ( i = 1; i <= m; i++ )
{
abi = i + ab * ( i % 2 );
abj = 2 * i + ab;
bj[i-1] = sqrt ( abi * abi / ( abj * abj - 1.0 ) );
}
}
else if ( kind == 8 )
{
ab = alpha + beta;
zemu = tgamma ( alpha + 1.0 ) * tgamma ( - ( ab + 1.0 ) )
/ tgamma ( - beta );
apone = alpha + 1.0;
aba = ab * apone;
aj[0] = - apone / ( ab + 2.0 );
bj[0] = - aj[0] * ( beta + 1.0 ) / ( ab + 2.0 ) / ( ab + 3.0 );
for ( i = 2; i <= m; i++ )
{
abti = ab + 2.0 * i;
aj[i-1] = aba + 2.0 * ( ab + i ) * ( i - 1 );
aj[i-1] = - aj[i-1] / abti / ( abti - 2.0 );
}
for ( i = 2; i <= m - 1; i++ )
{
abti = ab + 2.0 * i;
bj[i-1] = i * ( alpha + i ) / ( abti - 1.0 ) * ( beta + i )
/ ( abti * abti ) * ( ab + i ) / ( abti + 1.0 );
}
bj[m-1] = 0.0;
for ( i = 0; i < m; i++ )
{
bj[i] = sqrt ( bj[i] );
}
}
return zemu;
}
//****************************************************************************80
void imtqlx ( int n, double d[], double e[], double z[] )
//****************************************************************************80
//
// 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 (essentialy) of applying the orthogonal transformations
// directly to Z as they are generated.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 08 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.
//
{
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 prec;
double r;
double s;
prec = r8_epsilon ( );
if ( n == 1 )
{
return;
}
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] ) <= prec * ( fabs ( d[m-1] ) + fabs ( d[m] ) ) )
{
break;
}
}
p = d[l-1];
if ( m == l )
{
break;
}
if ( itn <= j )
{
cout << "\n";
cout << "IMTQLX - Fatal error!\n";
cout << " Iteration limit exceeded\n";
exit ( 1 );
}
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 ) * r8_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;
}
//****************************************************************************80
void parchk ( int kind, int m, double alpha, double beta )
//****************************************************************************80
//
// Purpose:
//
// PARCHK checks parameters ALPHA and BETA for classical weight functions.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 07 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.
//
// Parameters:
//
// Input, int KIND, the rule.
// 1, Legendre, (a,b) 1.0
// 2, Chebyshev, (a,b) ((b-x)*(x-a))^(-0.5)
// 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha
// 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta
// 5, Generalized Laguerre, (a,inf) (x-a)^alpha*exp(-b*(x-a))
// 6, Generalized Hermite, (-inf,inf) |x-a|^alpha*exp(-b*(x-a)^2)
// 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha
// 8, Rational, (a,inf) (x-a)^alpha*(x+b)^beta
//
// Input, int M, the order of the highest moment to
// be calculated. This value is only needed when KIND = 8.
//
// Input, double ALPHA, BETA, the parameters, if required
// by the value of KIND.
//
{
double tmp;
if ( kind <= 0 )
{
cout << "\n";
cout << "PARCHK - Fatal error!\n";
cout << " KIND <= 0.\n";
exit ( 1 );
}
//
// Check ALPHA for Gegenbauer, Jacobi, Laguerre, Hermite, Exponential.
//
if ( 3 <= kind && alpha <= -1.0 )
{
cout << "\n";
cout << "PARCHK - Fatal error!\n";
cout << " 3 <= KIND and ALPHA <= -1.\n";
exit ( 1 );
}
//
// Check BETA for Jacobi.
//
if ( kind == 4 && beta <= -1.0 )
{
cout << "\n";
cout << "PARCHK - Fatal error!\n";
cout << " KIND == 4 and BETA <= -1.0.\n";
exit ( 1 );
}
//
// Check ALPHA and BETA for rational.
//
if ( kind == 8 )
{
tmp = alpha + beta + m + 1.0;
if ( 0.0 <= tmp || tmp <= beta )
{
cout << "\n";
cout << "PARCHK - Fatal error!\n";
cout << " KIND == 8 but condition on ALPHA and BETA fails.\n";
exit ( 1 );
}
}
return;
}
//****************************************************************************80
double r8_epsilon ( )
//****************************************************************************80
//
// Purpose:
//
// R8_EPSILON returns the R8 roundoff unit.
//
// Discussion:
//
// The roundoff unit is a number R which is a power of 2 with the
// property that, to the precision of the computer's arithmetic,
// 1 < 1 + R
// but
// 1 = ( 1 + R / 2 )
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 01 September 2012
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, double R8_EPSILON, the R8 round-off unit.
//
{
const double value = 2.220446049250313E-016;
return value;
}
//****************************************************************************80
double r8_huge ( )
//****************************************************************************80
//
// Purpose:
//
// R8_HUGE returns a "huge" R8.
//
// Discussion:
//
// The value returned by this function is NOT required to be the
// maximum representable R8. This value varies from machine to machine,
// from compiler to compiler, and may cause problems when being printed.
// We simply want a "very large" but non-infinite number.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 06 October 2007
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Output, double R8_HUGE, a "huge" R8 value.
//
{
double value;
value = 1.0E+30;
return value;
}
//****************************************************************************80
double r8_sign ( double x )
//****************************************************************************80
//
// Purpose:
//
// R8_SIGN returns the sign of an R8.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 18 October 2004
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Input, double X, the number whose sign is desired.
//
// Output, double R8_SIGN, the sign of X.
//
{
double value;
if ( x < 0.0 )
{
value = -1.0;
}
else
{
value = 1.0;
}
return value;
}
//****************************************************************************80
void r8mat_write ( string output_filename, int m, int n, double table[] )
//****************************************************************************80
//
// Purpose:
//
// R8MAT_WRITE writes an R8MAT file with no header.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 29 June 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Input, string OUTPUT_FILENAME, the output filename.
//
// Input, int M, the spatial dimension.
//
// Input, int N, the number of points.
//
// Input, double TABLE[M*N], the table data.
//
{
int i;
int j;
ofstream output;
//
// Open the file.
//
output.open ( output_filename.c_str ( ) );
if ( !output )
{
cerr << "\n";
cerr << "R8MAT_WRITE - Fatal error!\n";
cerr << " Could not open the output file.\n";
return;
}
//
// Write the data.
//
for ( j = 0; j < n; j++ )
{
for ( i = 0; i < m; i++ )
{
output << " " << setw(24) << setprecision(16) << table[i+j*m];
}
output << "\n";
}
//
// Close the file.
//
output.close ( );
return;
}
//****************************************************************************80
void rule_write ( int order, string filename, double x[], double w[],
double r[] )
//****************************************************************************80
//
// Purpose:
//
// RULE_WRITE writes a quadrature rule to three files.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 18 February 2010
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// Input, int ORDER, the order of the rule.
//
// Input, double A, the left endpoint.
//
// Input, double B, the right endpoint.
//
// Input, string FILENAME, specifies the output filenames.
// "filename_w.txt", "filename_x.txt", "filename_r.txt"
// defining weights, abscissas, and region.
//
{
string filename_r;
string filename_w;
string filename_x;
int i;
int kind;
filename_w = filename + "_w.txt";
filename_x = filename + "_x.txt";
filename_r = filename + "_r.txt";
cout << "\n";
cout << " Creating quadrature files.\n";
cout << "\n";
cout << " Root file name is \"" << filename << "\".\n";
cout << "\n";
cout << " Weight file will be \"" << filename_w << "\".\n";
cout << " Abscissa file will be \"" << filename_x << "\".\n";
cout << " Region file will be \"" << filename_r << "\".\n";
r8mat_write ( filename_w, 1, order, w );
r8mat_write ( filename_x, 1, order, x );
r8mat_write ( filename_r, 1, 2, r );
return;
}
//****************************************************************************80
void scqf ( int nt, double t[], int mlt[], double wts[], int nwts, int ndx[],
double swts[], double st[], int kind, double alpha, double beta, double a,
double b )
//****************************************************************************80
//
// Purpose:
//
// SCQF scales a quadrature formula to a nonstandard interval.
//
// Discussion:
//
// The arrays WTS and SWTS may coincide.
//
// The arrays T and ST may coincide.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 16 February 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.
//
// Parameters:
//
// Input, int NT, the number of knots.
//
// Input, double T[NT], the original knots.
//
// Input, int MLT[NT], the multiplicity of the knots.
//
// Input, double WTS[NWTS], the weights.
//
// Input, int NWTS, the number of weights.
//
// Input, int NDX[NT], used to index the array WTS.
// For more details see the comments in CAWIQ.
//
// Output, double SWTS[NWTS], the scaled weights.
//
// Output, double ST[NT], the scaled knots.
//
// Input, int KIND, the rule.
// 1, Legendre, (a,b) 1.0
// 2, Chebyshev Type 1, (a,b) ((b-x)*(x-a))^(-0.5)
// 3, Gegenbauer, (a,b) ((b-x)*(x-a))^alpha
// 4, Jacobi, (a,b) (b-x)^alpha*(x-a)^beta
// 5, Generalized Laguerre, (a,+oo) (x-a)^alpha*exp(-b*(x-a))
// 6, Generalized Hermite, (-oo,+oo) |x-a|^alpha*exp(-b*(x-a)^2)
// 7, Exponential, (a,b) |x-(a+b)/2.0|^alpha
// 8, Rational, (a,+oo) (x-a)^alpha*(x+b)^beta
// 9, Chebyshev Type 2, (a,b) ((b-x)*(x-a))^(+0.5)
//
// Input, double ALPHA, the value of Alpha, if needed.
//
// Input, double BETA, the value of Beta, if needed.
//
// Input, double A, B, the interval endpoints.
//
{
double al;
double be;
int i;
int k;
int l;
double p;
double shft;
double slp;
double temp;
double tmp;
temp = r8_epsilon ( );
parchk ( kind, 1, alpha, beta );
if ( kind == 1 )
{
al = 0.0;
be = 0.0;
if ( fabs ( b - a ) <= temp )
{
cout << "\n";
cout << "SCQF - Fatal error!\n";
cout << " |B - A| too small.\n";
exit ( 1 );
}
shft = ( a + b ) / 2.0;
slp = ( b - a ) / 2.0;
}
else if ( kind == 2 )
{
al = -0.5;
be = -0.5;
if ( fabs ( b - a ) <= temp )
{
cout << "\n";
cout << "SCQF - Fatal error!\n";
cout << " |B - A| too small.\n";
exit ( 1 );
}
shft = ( a + b ) / 2.0;
slp = ( b - a ) / 2.0;
}
else if ( kind == 3 )
{
al = alpha;
be = alpha;
if ( fabs ( b - a ) <= temp )
{
cout << "\n";
cout << "SCQF - Fatal error!\n";
cout << " |B - A| too small.\n";
exit ( 1 );
}
shft = ( a + b ) / 2.0;
slp = ( b - a ) / 2.0;
}
else if ( kind == 4 )
{
al = alpha;
be = beta;
if ( fabs ( b - a ) <= temp )
{
cout << "\n";
cout << "SCQF - Fatal error!\n";
cout << " |B - A| too small.\n";
exit ( 1 );
}
shft = ( a + b ) / 2.0;
slp = ( b - a ) / 2.0;
}
else if ( kind == 5 )
{
if ( b <= 0.0 )
{
cout << "\n";
cout << "SCQF - Fatal error!\n";
cout << " B <= 0\n";
exit ( 1 );
}
shft = a;
slp = 1.0 / b;
al = alpha;
be = 0.0;
}
else if ( kind == 6 )
{
if ( b <= 0.0 )
{
cout << "\n";
cout << "SCQF - Fatal error!\n";
cout << " B <= 0.\n";
exit ( 1 );
}
shft = a;
slp = 1.0 / sqrt ( b );
al = alpha;
be = 0.0;
}
else if ( kind == 7 )
{
al = alpha;
be = 0.0;
if ( fabs ( b - a ) <= temp )
{
cout << "\n";
cout << "SCQF - Fatal error!\n";
cout << " |B - A| too small.\n";
exit ( 1 );
}
shft = ( a + b ) / 2.0;
slp = ( b - a ) / 2.0;
}
else if ( kind == 8 )
{
if ( a + b <= 0.0 )
{
cout << "\n";
cout << "SCQF - Fatal error!\n";
cout << " A + B <= 0.\n";
exit ( 1 );
}
shft = a;
slp = a + b;
al = alpha;
be = beta;
}
else if ( kind == 9 )
{
al = 0.5;
be = 0.5;
if ( fabs ( b - a ) <= temp )
{
cout << "\n";
cout << "SCQF - Fatal error!\n";
cout << " |B - A| too small.\n";
exit ( 1 );
}
shft = ( a + b ) / 2.0;
slp = ( b - a ) / 2.0;
}
p = pow ( slp, al + be + 1.0 );
for ( k = 0; k < nt; k++ )
{
st[k] = shft + slp * t[k];
l = abs ( ndx[k] );
if ( l != 0 )
{
tmp = p;
for ( i = l - 1; i <= l - 1 + mlt[k] - 1; i++ )
{
swts[i] = wts[i] * tmp;
tmp = tmp * slp;
}
}
}
return;
}
//****************************************************************************80
void sgqf ( int nt, double aj[], double bj[], double zemu, double t[],
double wts[] )
//****************************************************************************80
//
// Purpose:
//
// SGQF computes knots and weights of a Gauss Quadrature formula.
//
// Discussion:
//
// This routine computes all the knots and weights of a Gauss quadrature
// formula with simple knots from the Jacobi matrix and the zero-th
// moment of the weight function, using the Golub-Welsch technique.
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 08 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.
//
// Parameters:
//
// Input, int NT, the number of knots.
//
// Input, double AJ[NT], the diagonal of the Jacobi matrix.
//
// Input/output, double BJ[NT], the subdiagonal of the Jacobi
// matrix, in entries 1 through NT-1. On output, BJ has been overwritten.
//
// Input, double ZEMU, the zero-th moment of the weight function.
//
// Output, double T[NT], the knots.
//
// Output, double WTS[NT], the weights.
//
{
int i;
//
// Exit if the zero-th moment is not positive.
//
if ( zemu <= 0.0 )
{
cout << "\n";
cout << "SGQF - Fatal error!\n";
cout << " ZEMU <= 0.\n";
exit ( 1 );
}
//
// Set up vectors for IMTQLX.
//
for ( i = 0; i < nt; i++ )
{
t[i] = aj[i];
}
wts[0] = sqrt ( zemu );
for ( i = 1; i < nt; i++ )
{
wts[i] = 0.0;
}
//
// Diagonalize the Jacobi matrix.
//
imtqlx ( nt, t, bj, wts );
for ( i = 0; i < nt; i++ )
{
wts[i] = wts[i] * wts[i];
}
return;
}
//****************************************************************************80
void timestamp ( )
//****************************************************************************80
//
// Purpose:
//
// TIMESTAMP prints the current YMDHMS date as a time stamp.
//
// Example:
//
// 31 May 2001 09:45:54 AM
//
// Licensing:
//
// This code is distributed under the GNU LGPL license.
//
// Modified:
//
// 08 July 2009
//
// Author:
//
// John Burkardt
//
// Parameters:
//
// None
//
{
# define TIME_SIZE 40
static char time_buffer[TIME_SIZE];
const struct std::tm *tm_ptr;
size_t len;
std::time_t now;
now = std::time ( NULL );
tm_ptr = std::localtime ( &now );
len = std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );
std::cout << time_buffer << "\n";
return;
# undef TIME_SIZE
}
//
References
