[lkdvos]

This PR speeds up the stable lower-weight CGC solve from #61. With dense factorizations, #61 made the per-weight cost cubic in the number of pairs. This PR brings the runtime back below that of the old, unstable code. It combines what were #64, #65 and #67, one commit per change:

1. Sparse Cholesky. B, the Gram matrix of the raising operators on the pairs |m₁, m₂⟩ of a weight, is very sparse. It is now built as a sparse R * Rᵀ from the nonzero entries of the raising operators. Blocks of at least 100 pairs with Float64 entries are factorized with CHOLMOD. Small blocks and other element types keep the dense Cholesky. This adds SparseArrays, a stdlib, as a dependency.
2. Dense per-weight blocks. The solution of every weight is kept as a dense block sols[p, α, j], together with its pairs and their indices. The right-hand sides read the parent CGCs from these blocks: one lookup in a small per-weight Dict per parent pair, instead of N123 lookups in the large SparseArray. The SparseArray is filled once at the end.
3. A single shifted system per weight. The Sylvester equation B * sols + sols * A = 2 * rhs * eqs from Numerically stable computation of lower-weight CGCs #61, with A = eqsᵀ * eqs, needed one factorization of B + μ * I for every eigenvalue μ of A. The lowering and raising equations are consistent: the true CGCs satisfy both exactly, so any positive weighting of them has the same least-squares solution. Weighting the lowering equations with c * eqs * A⁻² * eqsᵀ gives

   (B + c * I) * sols = rhs * eqs * (I + c * A⁻¹)
   

   This needs a single factorization per weight and no eigendecomposition. Here c = tr(A) / size(A, 1) is the mean eigenvalue of A.
4. Skip tiny entries. Entries with abs(v) <= TOL_PURGE are not written into the SparseArray, since purge! would remove them right afterwards. This leaves the result unchanged.

Timings

Single BLAS thread, CGCs computed directly (bypassing the caches):

case	before #61 (unstable)	main (#61, dense)	this PR
SU(2) (250,0)³	0.026 s	0.089 s	0.056 s
SU(3) (10,5,0)³	–	0.65 s	0.20 s
SU(3) (12,6,0)³	2.4 s	4.1 s	0.75 s
SU(3) (16,8,0)³	23 s	95 s	5.5 s
SU(4) (6,3,3,0)³	0.90 s	1.4 s	0.32 s

Very small cases, such as SU(3) (6,3,0)³ at about 10 ms, are slightly slower, because the Gram matrix is always built as a sparse matrix. About 60 % of the remaining time goes to the sparse Cholesky factorization.

Accuracy

Largest residual:

case	main	this PR
SU(3) (14,0,0) ⊗ (14,14,0) → (14,7,0), intertwiner	6.7e-15	1.8e-15
SU(3) (20,0,0) ⊗ (20,20,0) → (20,10,0), intertwiner / isometry	1.3e-14 / 1.7e-14	7.9e-15 / 2.4e-15
SU(3) (16,8,0)³, intertwiner / isometry	1.8e-14 / 1.2e-14	1.8e-14 / 8.0e-15
SU(4) (6,3,3,0)³, intertwiner / isometry	2.7e-15 / 1.2e-14	2.8e-15 / 2.9e-15

The full test suite passes.

Tried and left out

• Reverse Cuthill–McKee ordering instead of CHOLMOD's AMD: 3–10 % faster, but about 25 more lines of code.
• Simplicial factorization instead of supernodal: slower.
• LDLᵀ instead of Cholesky: slower.
• Skipping pairs reached only from parent pairs with vanishing CGCs: this breaks SU(2) for spins ≥ 100, because gaugefix! zeroes small highest-weight CGCs (fixed in Keep small highest-weight CGCs in gaugefix! #66).

🤖 Generated with Claude Code

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

Files with missing lines	Coverage Δ
src/SUNRepresentations.jl	100.00% <ø> (ø)
src/clebschgordan.jl	97.90% <100.00%> (+0.05%)	⬆️

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.