Matrix-RHS adjoint/transpose matvec falls through to the generic dense fallback — ~3000x slowdown

WHAT: Only vector adjoint/transpose mul! are defined (src/matvec.jl:163-171). mul!(Y, A', X) / A'*X for a matrix X dispatches to LinearAlgebra's generic fallback, which materializes columns of A' via getindex. Measured (n=5000, 10 cols): 17.4 s for 20 calls vs 0.0058 s for a structured per-column loop (bit-identical) — a 3007x slowdown. WHY IT MATTERS: This is the exact '~1000x slower' cliff the matvec.jl comment warns about, fixed for vectors but not matrices. Any block/multi-RHS adjoint (block iterative least squares, lstsq's A'*res on a matrix) silently hits it, and the gap went unnoticed because the adjoint matvec is untested (see the coverage issue). FIX: Add mul!(Y::AbstractMatrix, wA::Adjoint/Transpose{<:SparseWithDenseRowColMatrix}, X::AbstractMatrix, α, β) — a matrix analogue of _adjoint_matvec!: mul!(Y, Sop(A.S), X, α, β) then the rank-r block correction Y .+= α Vᴴ(Uᴴ X) (cheap r×ncols block, r≪n), plus the Base.:* matrix wrapper. EFFORT: S-M.

Priority: high. Filed from an automated next-steps audit of the QR/lstsq work (see PR #6).

Start in src/matvec.jl:163-171 and read _adjoint_matvec! to understand the existing vector path. Add the matrix mul! path for Adjoint/Transpose and the Base.:* wrapper described in the issue, then compare matrix-RHS results and performance with the structured per-column loop; done means avoiding the generic dense fallback while preserving bit-identical results.