Repository navigation
numpy: native det, slogdet, solve and inv, bit-exact with OpenBLAS - #218
Merged
Merged
Conversation
Python files that use numpy's core now build with a numpy written in Rust (src/dynpy/rt_np.rs, rt_npfn.rs, rt_npfmt.rs): float64/int64/bool arrays with NEP 50 promotion, strided views, masks and integer-array indexing, broadcasting, operators, ~60 ufuncs, reductions in numpy's pairwise order, sort/unique/stacking, and numpy 2.x printing byte for byte. Anything it lacks (np.random, np.linalg, other dtypes) is refused at build time and the program gets CPython's numpy through the typed build, as before. The dynamic build no longer calls CPython's numpy itself (5 of 40 corpus programs gave wrong output that way). BREVI_NUMPY=0 turns it off. tests/numpy: 11 programs whose .expected files come from real numpy; the check fails if any of them runs off Brevi's numpy. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0136hkZgpjat3dn8BGZPJVF2
src/dynpy/rt_nprand.rs: default_rng(seed) as numpy does it (SeedSequence, PCG64, the ziggurat normal and exponential with numpy's tables, Lemire bounded integers) and the legacy np.random.seed API (MT19937 init_genrand, polar gauss with its cached value, masked randint). random, integers, normal, uniform, standard_normal, exponential, choice (replace=False, p=), shuffle, permutation, rand, randn, randint, RandomState. Other distributions and ufunc methods (np.maximum.accumulate) are refused at build time, so the program gets CPython's numpy. Unseeded generators take OS entropy without file access. tests/numpy: 6 random programs, expected output from real numpy. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0136hkZgpjat3dn8BGZPJVF2
…apezoid, ufunc reduce/accumulate Float vector np.dot now follows OpenBLAS's own summation order for the CPU (AVX2 or AVX-512, with its thread split), so it matches numpy to the last digit. det/solve/inv are written but stay off (go to CPython numpy) until they match OpenBLAS's LU bit for bit. Fixes str * np.int64 raising TypeError. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0136hkZgpjat3dn8BGZPJVF2
Follows OpenBLAS's single-thread LU step by step (getrf_single recursive blocking, getf2, dgemv_n, the SkylakeX/Haswell dgemm kernels' k-loop splits, trsm_kernel_LT/LN, trsv, daxpy), so results equal numpy's to the last digit on AVX-512 and AVX2 CPUs. slogdet returns SlogdetResult; matrix_power takes negative powers. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_0136hkZgpjat3dn8BGZPJVF2
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Requested by Adam · project thread
Fourth native numpy PR. It turns on
np.linalg.det,slogdet,solveandinv, which #217 had left off because their last digit sometimes differed from numpy's.How they match exactly
numpy calls OpenBLAS's
dgetrf/dgesv. The newlk_*code inrt_nplinalg.rsfollows OpenBLAS 0.3.30's single-thread path step by step:getrf_single.c: the recursive blocking, which depends onGEMM_UNROLL_NandGEMM_Qgetf2_k.cpanel: strided ddot, thendgemv_nwith each microkernel's own summation orderdgemmkernel k-loops: the SkylakeXm4n12,m2n12andm1n12kernels split k across 2 or 4 accumulators, and Haswell's 4x1 kernel splits k 4 waystrsm_kernel_LT/LNblocking,trsv_L/Uanddaxpyfor the solveI read each part from OpenBLAS's source and checked it against the real library through ctypes. The Python model is in
/mnt/project-files/brevi-numpy/openblas-model/.Checks
OPENBLAS_CORETYPE=Haswell) against a build with Brevi forced to its AVX2 model.LinAlgError: Singular matrix, det 0.0, slogdet(0.0, -inf)), NaN and inf entries, subnormal pivots, int inputs, and the shape errors.slogdetreturns numpy'sSlogdetResultnamed tuple.matrix_powernow takes negative powers.tests/numpy/linalg_lu.py. Its output is identical under numpy's SkylakeX, Haswell and Prescott kernels, so it passes on any CI CPU.sh tools/ci_local.sh: all checks passed, 0 regressions.Known limits (documented)
OPENBLAS_NUM_THREADS=1.Generated by Claude Code