Complete BLAS Level 1-3 API and add NaN guards to statistics - #20
Conversation
Documents module map (36,868 LOC across 50 modules), BLAS parity gaps (15/23 routines missing), quantized GEMM verification, NaN guard audit (3 unguarded divisions found), and not-yet-migrated upstream modules. https://claude.ai/code/session_01CdqyUTUfjKZuk8YGJzv6LB
rustynum's simd_avx512.rs = std::simd compat layer (F32x16, F64x8, U8x64 etc.) backed by stable core::arch. ndarray uses raw __m512 — works but x86_64-only. Port as src/backend/simd_compat.rs to unlock aarch64 and std::simd migration.
Cranelift JIT for scan param baking — real infrastructure, not a JSON parser.
…t scan opt array_windows (1.94) covers scan record-size baking. jitson's real purpose: compile graph topology → native function via Cranelift. Keep at P2.
hybrid.rs=3-stage pipeline, delta.rs=XOR overlay, layer_stack.rs=collapse gate, soaking.rs=int8 accumulation, tail_backend.rs=libCEED trait. All P1.
10 features cross-referenced. Compat layer unblocks 5 of them. gather/dispatch/prefetch exist internally but need user-facing API. SpatialArray3 and stencil are genuinely new P3 types.
…wiring array_struct.rs has 35 SimdOps dispatch calls that ndarray lacks. activations.rs sigmoid/softmax use mapv (scalar). vml.rs exp/log/sqrt use for loops (scalar). Backend kernels exist (kernels_avx512.rs) but aren't wired to hpc traits.
Port rustynum simd_avx512.rs → src/backend/simd_compat.rs. 11 types, 60 impl blocks, zero runtime cost. Refactor kernels_avx512.rs + bitwise.rs to use compat types. Wire activations.rs + vml.rs through SIMD dispatch.
statistics.rs var_axis(), cascade.rs warmup, bf16_truth.rs awareness_classify() https://claude.ai/code/session_01CdqyUTUfjKZuk8YGJzv6LB
Reconciles naming convention with rustyblas. https://claude.ai/code/session_01CdqyUTUfjKZuk8YGJzv6LB
Level 1: scal, asum, iamax, swap, copy, rotg Level 2: trsv, symv, syr, syr2, gbmv, sbmv Level 3: trmm, trsm, symm https://claude.ai/code/session_01CdqyUTUfjKZuk8YGJzv6LB
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 22bfb7ab01
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| Uplo::Lower => { | ||
| if j <= i { | ||
| self[[j.abs_diff(i), j.max(i)]] | ||
| } else { | ||
| self[[i.abs_diff(j), i.max(j)]] |
There was a problem hiding this comment.
Index lower-band
sbmv data from the smaller column
When blas_sbmv(..., Uplo::Lower, ...) is given a normally packed lower band, this branch reads off-diagonal entries from max(i, j) instead of the stored column min(i, j). For example, with k = 1, A[1,0] should come from the first subdiagonal slot for column 0, but this code fetches column 1 instead, so distinct subdiagonal values are shifted and the matrix-vector product is wrong for essentially every nontrivial lower-stored band matrix.
Useful? React with 👍 / 👎.
| let scale = a.abs() + b.abs(); | ||
| let r = scale * ((a / scale).powi(2) + (b / scale).powi(2)).sqrt(); |
There was a problem hiding this comment.
Scale
blas_rotg without summing magnitudes first
This normalization overflows for large but still representable inputs. In f64, blas_rotg(1e308, 1e308) should produce a finite rotation, but a.abs() + b.abs() becomes inf, so both normalized terms collapse to zero and r becomes NaN (inf * 0). That propagates to c and s, making the new API fail exactly in the high-magnitude cases where a stable Givens implementation is expected to keep working.
Useful? React with 👍 / 👎.
Summary
This PR completes the BLAS Level 1-3 trait APIs in ndarray's HPC module and hardens numerical stability by adding guards against division-by-zero in statistics and awareness classification.
Key Changes
BLAS API Completion
Level 1 (
blas_level1.rs)blas_rotg()function andGivensRotation<A>struct for Givens rotation computationaandb, computes rotation parameters(r, c, s)such that the rotation zeros out the second componentLevel 2 (
blas_level2.rs)blas_syr()— symmetric rank-1 update:A = alpha * x * x^T + Ablas_syr2()— symmetric rank-2 update:A = alpha * x * y^T + alpha * y * x^T + Ablas_gbmv()— general banded matrix-vector multiply with band storage formatblas_sbmv()— symmetric banded matrix-vector multiplyUplo(Upper/Lower) triangle specificationLevel 3 (
blas_level3.rs)blas_trmm()— triangular matrix-matrix multiplySide::Left(alpha * A * B) andSide::Right(alpha * B * A)Uplotriangle specification for the triangular matrixNumerical Stability Hardening
statistics.rsvar_axis()to prevent NaN from division by zerovar_axis_zero_length_axis_no_nanbf16_truth.rsawareness_classify()to prevent division by zero when computing percentagesSuperpositionStatewhenn_dims == 0cascade.rssigma_popis 0.0 whenwarmup_n == 0instead of computing0.0/0.0Documentation & Constants
quantized.rsBF16::ZEROandBF16::ONEconstants for common valuesfrom_f32_truncate()alias to clarify truncation semantics (matching rustyblas naming convention)from_f32()to explicitly document truncation behaviorImplementation Details
A(i,j)stored atband[ku + i - j, j]Uploenum to determine which half of the matrix to accessTesting
https://claude.ai/code/session_01CdqyUTUfjKZuk8YGJzv6LB