methods

Is imtqlx the absolute best?

No. It is a good mid-scale solver for the Golub-Welsch eigenproblem with the ACM 655 trick (only track first eigenvector components / weight vector). It is not state of the art for large n, and its Givens chase is inherently serial.

What ACM 655 optimizes vs gw:

Landscape**

Method Idea Best for Parallelism
rec / rec_caf Hale-Townsend style Newton on P_n large n, moderate α,β excellent (per node)
gw GW + LAPACKE dsteqr robust mid n, high α,β single rule: poor; multi-rule batch: good
algo665 GW + sequential imtqlx mid n, high α,β, faster than gw poor (one QL path)
algo665_dc Cuppen D&C + imtqlx leaves larger n, same GW problem recursive split + leaf QL
sturm / sturm_caf Sturm bisection + inv. iteration GW-level, CAF-native indices excellent (per eigenvalue)
glr Prüfer phase-by-π march + Newton full Jacobi (high α,β via Bessel first-root) sequential chain
glr_caf independent k-th starters + Newton multi-image single rule excellent (per root index)
bogaert / bogaert_caf J0-zero asymptotics + exact weights Legendre only (α=β=0, n≥21) excellent (half-node partition)
auto select_method_auto policy default public entry delegates to above

Shipped algorithms**

Algorithm Status Notes
Cuppen D&C algo665_dc / algo665_dc_caf parallel structure around imtqlx leaves
Sturm + inv. iter. sturm / sturm_caf best CAF fit for the GW eigenproblem
GLR (serial) glr Sequential Prüfer march (Δθ=π) + Newton
GLR (CAF) glr_caf Independent k-th starters + Newton (index partition)
Bogaert bogaert / bogaert_caf Gauss-Legendre only; half-node CAF; error stop if non-Legendre

Public entry (preferred)**

Call gauss_jacobi_rule(npts, alpha, beta, x, wts [, method]):

Legacy gauss_jacobi(..., method) still works (required method string) and forwards to gauss_jacobi_rule.

Auto policy (select_method_auto)**

Regime Single-image Multi-image (num_images()>1)
α=β=0 and n≥21 bogaert bogaert_caf
n≥128 and max(\ α\ ,\
n≥64 and max(\ α\ ,\
max(\ α\ ,\
else sturm sturm_caf

Policy is documented and safe, not oracle-optimal. auto selects *_caf names when more than one image is present.

Practical defaults**

Regime Prefer
Legendre, n≥21 bogaert (or omit method)
large n, mild α,β rec / rec_caf
high α or β, mid n algo665 or sturm
multi-image, high α/β sturm_caf (per-eig) or gauss_jacobi_batch_caf (per-rule)
multi-image, mild α/β rec_caf or glr_caf

Math note (Sturm)**

For tridiagonal T with diagonal a_i, off-diagonal b_i, the Sturm count (number of eigenvalues =< λ) is obtained from the LDL-style recurrence

d_1 = a_1 - λ, d_i = (a_i - λ) - b_{i-1}^2 / d_{i-1},

counting negative pivots. The k-th eigenvalue is found by bisection on that count. Weights use a few steps of inverse iteration on (T - λI) and w = mu0 * u_1^2.

See also scripts/imtqlx_parallel_math.py for why a single QL sweep cannot be reordered, and how Cuppen creates independent leaves.

SymPy accuracy reference

High-precision checks use sympy.integrals.quadrature.gauss_jacobi via scripts/sympy_accuracy.py (and pytests/test_sympy_accuracy.py). The harness calls the shipped library entry (gauss_jacobi Python/C ABI → Fortran kernels) for each cell of a regime matrix (mild Jacobi, high α/β, Legendre n≥21, auto); it does not hard-code golden node tables. Out-of-regime bogaert (α≠0 or β≠0) is asserted as a policy error, not an accuracy pass.

# after meson build of libgjp_cinterp
export GJP_CINTERP=$(find bbdir -name libgjp_cinterp.so | head -1)
python scripts/sympy_accuracy.py --table sympy_accuracy_table.txt
pytest pytests/test_sympy_accuracy.py -v

References (ookcite collection GaussJacobiQuad)

Machine-managed bibliography: docs/refs.bib (export from ookcite collection GaussJacobiQuad). DOIs validated; keys are Better-BibTeX style from ookcite.

Method / topic Primary refs (BibTeX keys)
Software cite goswamiHaozekegaussjacobiquadGaussjacobiquadI2023
gw Golub–Welsch golubCalculationGaussQuadrature1969
algo665 / imtqlx elhayAlgorithm6551987, martinImplicitqlAlgorithm1968
algo665_dc Cuppen cuppenDivideConquerMethod1980, guDivideandconquerAlgorithmSymmetric1995
rec Hale–Townsend haleFastAccurateComputation2013
bogaert bogaertIterationfreeComputationGausslegendre2014
glr Prüfer phase glaserFastAlgorithmCalculation2007

Bibliography (APA, from ookcite):

  1. Bogaert, I. (2014). Iteration-Free Computation of Gauss–Legendre Quadrature Nodes and Weights. SIAM Journal on Scientific Computing, /36/(3), A1008–A1026. https://doi.org/10.1137/140954969
  2. Cuppen, J. J. M. (1980). A divide and conquer method for the symmetric tridiagonal eigenproblem. Numerische Mathematik, /36/(2), 177–195. https://doi.org/10.1007/bf01396757
  3. Elhay, S., & Kautsky, J. (1987). Algorithm 655. ACM Transactions on Mathematical Software, /13/(4), 399–415. https://doi.org/10.1145/35078.214351
  4. Glaser, A., Liu, X., & Rokhlin, V. (2007). A Fast Algorithm for the Calculation of the Roots of Special Functions. SIAM Journal on Scientific Computing, /29/(4), 1420–1438. https://doi.org/10.1137/06067016x
  5. Golub, G. H., & Welsch, J. H. (1969). Calculation of Gauss quadrature rules. Mathematics of Computation, /23/(106), 221–230. https://doi.org/10.1090/s0025-5718-69-99647-1
  6. Goswami, R. (2023). HaoZeke/GaussJacobiQuad: GaussJacobiQuad I. Zenodo. https://doi.org/10.5281/zenodo.8285112
  7. Gu, M., & Eisenstat, S. C. (1995). A Divide-and-Conquer Algorithm for the Symmetric Tridiagonal Eigenproblem. SIAM Journal on Matrix Analysis and Applications, /16/(1), 172–191. https://doi.org/10.1137/s0895479892241287
  8. Hale, N., & Townsend, A. (2013). Fast and Accurate Computation of Gauss–Legendre and Gauss–Jacobi Quadrature Nodes and Weights. SIAM Journal on Scientific Computing, /35/(2), A652–A674. https://doi.org/10.1137/120889873
  9. Martin, R. S., & Wilkinson, J. H. (1968). The implicit QL algorithm. Numerische Mathematik, /12/(5), 377–383. https://doi.org/10.1007/bf02161360