Function solve_tridiagonal 

fn solve_tridiagonal<R: Scalar>(
    sub: &[f64],
    diag: &[f64],
    sup: &[f64],
    d: &[R],
) -> Vec<R>

Solves a tridiagonal system A x = d by the Thomas algorithm, with the matrix passive and the right-hand side in the pricing mode.

sub, diag and sup are the three bands — sub[0] and sup[n - 1] are unused — and every one of them is f64, because a spline’s matrix depends on the knot spacing alone. That is the whole reason an interpolating spline can record exact adjoints: the elimination’s divisors are passive constants, so the solve is a linear map from the ordinates to the coefficients and the mode’s arithmetic differentiates the ordinates and nothing else. The forward sweep’s f64 half is computed once into sup_p and never reaches the tape.

A scheme whose bands depended on the ordinates could not use this, and that is exactly the family that has to declare an adjoint refusal.