Function initial_solve 

fn initial_solve(
    n: [usize; 3],
    p: &[f64],
    w: &[f64],
    dk: &[f64],
    crop_box: &Box3,
    max_iter: usize,
    tracer: &impl Tracer,
) -> (Vec<f64>, LsqrOutcome)

Step 3 of ilsqr: LSQR (tol 0.01) on D W D χ = D W φ with D = F⁻¹ diag(dk) F. Returns χ and how LSQR ended.

Solved in k-space, for x̃ = F χ: with the unitary U = F/√N, D W D = U* (dk·U W U*·dk) U, so the LSQR iterates for (U A U*, U b) are U times those for (A, b) and the norms are equal (exact arithmetic). The √N factors cancel in dk·F W F⁻¹·dk, and taking F b = √N·U b as the right-hand side (LSQR is linear in b and its stopping tests are scale-free) leaves χ = F⁻¹ x̃. That is 2 FFTs per operator call instead of 4; rounding differs only at the level of the transforms’ own rounding. The products with dk and W are fused into the transforms’ first and last passes, which only cover the crop box (W = 0 outside the mask).