f77 [ flags ] file(s) ... -L/usr/local/lib -lgjlC (K&R, 89, 99), C++ (98):SUBROUTINE glqrc (a, b, s, t, alpha, nquad, ierr) DOUBLE PRECISION a(0:MAXPTS), alpha, b(0:MAXPTS) DOUBLE PRECISION s(0:MAXPTS), t(0:MAXPTS) INTEGER ierr, nquad
cc [ flags ] -I/usr/local/include file(s) ... -L/usr/local/lib -lgjl
Use#include <gjl.h>to get this prototype:void glqrc(fortran_double_precision a_[], fortran_double_precision b_[], fortran_double_precision s_[], fortran_double_precision t_[], const fortran_double_precision * alpha_, const fortran_integer * nquad_, fortran_integer * ierr_);
NB: The definition of C/C++ data types fortran_ xxx, and the mapping of Fortran external names to C/C++ external names, is handled by the C/C++ header file. That way, the same function or subroutine name can be used in C, C++, and Fortran code, independent of compiler conventions for mangling of external names in these programming languages.
Compute the recursion coefficients and zeroth and first moments of the monic polynomials corresponding to the positive weight function
w(x,\alpha) = (x \(mi 1 \(mi ln(x)) * exp(\(mix) * x^\alpha
with recursion relation (n = 0, 1, 2, ...)
P_{n+1}^\alpha(x) = (x \(mi B_n^\alpha) * P_n^\alpha(x) \(mi A_n^\alpha * P_{n\(mi1}^\alpha(x)
and initial conditions
P_{\(mi1}^\alpha(x) = 0 P_{0}^\alpha(x) = 1
Except in the weight function, the superscripts indicate dependence on \alpha, NOT exponentiation.
The required moments are:
T_n^\alpha = \int_0^\infty w(x,\alpha) (P_n^\alpha(x))^2 dx S_n^\alpha = \int_0^\infty w(x,\alpha) (P_n^\alpha(x))^2 x dx
From these moments, the recursion coefficents are computed as:
A_n^\alpha = T_n^\alpha / T_{n\(mi1}^\alpha B_n^\alpha = S_n^\alpha / T_n^\alpha
On entry:
- alpha
- Power of x in the integrand (alpha > \(mi1).
- nquad
- Number of quadrature points to compute. It must be less than the limit MAXPTS defined in the header file, maxpts.inc. The default value chosen there should be large enough for any realistic application.
On return:
- a(0..nquad)
- Recursion coefficients: a(n) = A_n^\alpha.
- b(0..nquad)
- Recursion coefficients: b(n) = B_n^\alpha.
- s(0..nquad)
- First moments: s(n) = S_n^\alpha
- t(0..nquad)
- Zeroth moments: t(n) = T_n^\alpha
- ierr
- Error indicator:
= 0 (success), 1 (eigensolution could not be obtained), 2 (destructive overflow), 3 (nquad out of range), 4 (alpha out of range).
Fast Gaussian Quadrature for Two Classes of Logarithmic Weight Functionsin ACM Transactions on Mathematical Software, Volume ??, Number ??, Pages ????--???? and ????--????, 20xx, by
andNelson H. F. Beebe Center for Scientific Computing University of Utah Department of Mathematics, 110 LCB 155 S 1400 E RM 233 Salt Lake City, UT 84112-0090 Tel: +1 801 581 5254 FAX: +1 801 581 4148 Email: beebe@math.utah.edu, beebe@acm.org, beebe@computer.org WWW URL: http://www.math.utah.edu/~beebe
James S. Ball University of Utah Department of Physics Salt Lake City, UT 84112-0830 USA Tel: +1 801 581 8397 FAX: +1 801 581 6256 Email: ball@physics.utah.edu WWW URL: http://www.physics.utah.edu/people/faculty/ball.html