Repository navigation
Reimplement dpnp.random using oneMKL device RNG API - #3147
Draft
ndgrigorian wants to merge 29 commits into
Draft
ndgrigorian wants to merge 29 commits into
ndgrigorian wants to merge 29 commits into
Conversation
- vm: check the 2-D binary contig dispatch table with both type ids, the row lookup could never be nullptr - vm: fix fmax/fmin docstrings copied from fmod - lapack: cast syevd eigenvalues output to RealT, fix gesvd typo and error message naming gesv - ext: fix "c-contiguos" typo, drop unused DispatchTable2 getters and unused queue/reserve overloads of the local memory size helpers - offset_utils: drop unused FourZeroOffsets_Indexer, FixedDimStridedIndexer and ThreeOffsets_FixedDimStridedIndexer - cmake: drop the dead BUILD_DPNP_EXT_SYCL branch, unused variables, unused *_headers interface targets, a no-op PIC property, the empty COMPILE_DEFINITIONS loop and commented-out include dirs; stop adding -fno-fast-math twice to elementwise_functions_type_utils.cpp
- Move the C++ guide to dpnp/libtensor/AGENTS.md and update the layout, ufunc steps and paths in the other AGENTS.md files - Add changelog entries for the moved headers, dpnp.get_include and the removed indexers; fix the get_include docstring and test it - Keep the dpnp sources first in the merged modules and instrument the shared helpers so coverage builds measure them; note the mixed FP flags - Drop the dead SYCL option of build_dpnp_tensor_ext and the no-op CMAKE_POSITION_INDEPENDENT_CODE property of dpnp_backend_c
simplify_iteration_space and elementwise_functions_type_utils.cpp are compiled with -fsycl in Coverage builds. Instrumenting their device code too made sycl-link fail with "Appending variables with different address spaces need to be linked!" for every module that links them.
Add a `_rng_impl` pybind11 extension that samples on the device with oneMKL's device Philox4x32x10 engine and SYCL ports of NumPy's distribution algorithms, and expose it as `dpnp.random.Generator`, `dpnp.random.Philox4x3210` and `dpnp.random.default_rng`. Each work-item owns a disjoint block of Philox counters and the counter is advanced on the host, so a seed gives the same stream on any device. The Generator implements random, uniform, standard_normal, normal, lognormal, standard_exponential, exponential, standard_gamma, gamma, chisquare, beta, integers and bytes; other methods raise NotImplementedError. The legacy RNG is left in place for now.
Add SYCL samplers for the rest of NumPy's Generator distributions: standard_cauchy, standard_t, wald, pareto, weibull, power, rayleigh, laplace, gumbel, logistic, triangular, vonmises, noncentral_chisquare, f and noncentral_f, and the discrete poisson, binomial, negative_binomial, geometric, logseries, zipf and hypergeometric. Discrete samplers compute in double precision when the device supports it and in single precision otherwise, and always return int64. Build dirichlet, multinomial, multivariate_hypergeometric, multivariate_normal, permutation, shuffle, permuted and choice from these kernels and dpnp functions, so no method raises NotImplementedError any more. Parameters given on the host are validated with NumPy's exceptions and messages; parameters already on the device are not checked, and invalid values give undefined results.
RandomState and the module-level functions now run on the same device kernels as Generator, with NumPy's legacy parameter checks and error messages. The module functions keep one RandomState per SYCL queue and take the queue of a dpnp array param when no device is given. Parameters given as device arrays are not checked, so the rejection samplers guard against inputs they would loop on forever: poisson, binomial and hypergeometric return INT64_MIN and beta returns NaN. API changes: device, usm_type and sycl_queue are keyword-only; randn no longer takes d0=; uniform no longer takes an integer dtype; randint(dtype=int) returns int64; shuffle requires a dpnp.ndarray.
Parameters on the device are not checked, and for p == 1 the sampler rejected every draw, so the kernel ran until the GPU watchdog reset the device. The CuPy test for an invalid p hit it. Return INT64_MIN as the other samplers do for invalid parameters.
A scalar parameter of None raised ValueError about NaN, as it was taken as an array; NumPy raises TypeError "must be real number" on its scalar path, and takes None as NaN next to array parameters, which failed in dpnp.asarray. Convert and check host parameters in NumPy's order. integers and randint did not check bounds given as host arrays, so an empty range gave garbage, and an exclusive high of 2**63 in a list overflowed. Check them as NumPy does, with its "high <= 0" message for a zero low, and shift them as exact integers. RandomState.multinomial now converts pvals before checking its shape and raises NumPy's TypeError for an n that is not an integer, and random_integers takes a scalar high as NumPy does.
dpnp.random now runs entirely on the device RNG extension, so drop the legacy implementation: the dpnp_algo_random Cython module, the backend RNG kernels, DPNPRandomState, the DPNP_FN_RNG_* function pointers, the mt19937/mcg59 engines of backend_sycl, and their C++ examples and tests. The backend library no longer links oneMKL RNG, so the conda recipes do not pin onemkl-sycl-rng. Rewrite the random reference page for the Generator and RandomState APIs, with the differences from NumPy. Enable the CuPy random tests. Skip the ones for bit generators that are not provided, for checks of parameters on the device, and for CuPy behaviors that differ from NumPy's; update the mock tests to the call signatures of the new module functions. test_generator.py stays skipped, as it tests the internals of CuPy's legacy generator.
DPC++ 2026.2 miscompiles sycl::mul_hi(uint64, uint64) on Intel GPUs when one operand is a compile-time constant: at -O2 every result is wrong, so the PCG64DXSM step (multiplying by the constant cheap multiplier) diverged from NumPy after the first element. Compute the high word from 32-bit halves instead, as oneMKL's device RNG does, in the PCG64DXSM step and in the 64-bit Lemire bounded integer.
Drop the skips of features the device RNG now provides: generator distributions, binomial with array parameters, multivariate_normal, RandomState.permutation, zero-dimensional shuffle and permutation, small and unsigned integer dtypes in shuffle and permutation, and randint over a full integer range. Adapt two upstream lines that cannot run on dpnp (`.get()` and `numpy.count_nonzero` of a dpnp.ndarray), and fix the multinomial dtype assertion, which fails for NumPy as well. Give the remaining skips their actual reasons: they test CuPy internals or CuPy-only API.
Remove the skipped_tests_cuda.tbl entries for test functions the device RNG rewrite deleted. The entries for tests that still exist stay until the device RNG has been run on a CUDA device.
The rattler-build recipe still pinned onemkl-sycl-rng, which nothing links anymore, and the agent notes said dpnp.random uses dpnp_backend_c. Also drop the Windows workaround of #2871, which only kept the oneMKL 2026.0 mt19937 engine of backend_sycl from being destroyed at DLL detach; backend_sycl no longer holds an engine.
Restore the upstream bodies of the CuPy random tests and testing helpers that the device RNG work had edited, and skip the tests that cannot run on dpnp instead: the mock tests asserting how arguments are passed, `.get()`, `numpy.count_nonzero` of a dpnp.ndarray, the multinomial dtype.kind assertion, cupy.cuda.Device and scalar Dirichlet alpha. The HIP skipif, which would fail at import, is commented out. Name the module's state dictionary `_dpnp_random_states` again, as the testing helper expects.
As in NumPy, they get and swap the bit generator the dpnp.random functions draw from, here that of the state of a queue: get takes the device and sycl_queue keywords, and set uses the queue of the bit generator. The bit generator itself is used, so the functions and a Generator of it share one stream. RandomState now takes a BitGenerator as its seed, as in NumPy. Like NumPy, RandomState.seed reseeds only the legacy bit generator, Philox4x3210, while dpnp.random.seed reseeds any bit generator in place, so get_bit_generator keeps returning the same object.
dpnp_backend_c now only provides the 1-D partition kernel behind dpnp.partition. Remove its ones, zeros, *_like and initval kernels, the non-extended partition entry with its synchronous wrapper on the global default queue, the memory helpers and their pointer adapter, the queue and verbose-mode singletons, the empty GTest suite with its DPNP_BACKEND_TESTS option, and the link to oneMKL, which nothing calls anymore. Also remove the unused Cython helpers call_origin, copy_from_origin, unwrap_array, _object_to_tuple and get_common_usm_allocation.
Generator and RandomState silently copied parameter arrays from another SYCL queue onto the generator's queue. permutation, shuffle, permuted and choice instead moved the index array to the input's queue. Every method taking arrays now calls Generator._check_queues first and raises ExecutionPlacementError on a mismatch, before any device-side parameter check. With no device given, the module functions use the common queue of their array parameters, and raise the same error before creating a state when the parameters are on incompatible queues. The temporaries of permutation, permuted and choice now take the input's or the requested USM type. Before, they were always "device", which made permutation and permuted return device memory and made choice ignore usm_type without replacement or with p. Fix TestSetBitGenerator::test_class leaving numpy.random with a class as its global bit generator. NumPy stores it before raising, which broke every later numpy.random.get_state() in the session. Add Generator, bit generator, permutation and queue-mismatch tests to test_sycl_queue.py, and USM type tests to test_usm_type.py. Disable the upstream xfail on CuPy's test_permutation_zero_dim, which now passes. Empty skipped_tests_cuda.tbl. Its entries only skipped the NotImplementedError the legacy RNG raised on CUDA. The device RNG has not been run on CUDA yet. Describe the dpnp.random layout, the kernel rules of the rng extension and the random test rules in the AGENTS.md files.
ndgrigorian
changed the title
dpnp.random using oneMKL device RNG API
Contributor
|
View rendered docs @ https://intelpython.github.io/dpnp/pull/3147/index.html |
Contributor
|
Array API standard conformance tests for dpnp=0.22.0dev0=np2py314h8d9cdd5_41 ran successfully. |
Collaborator
|
Coverage is 82.811% — feature/use-mkl-device-rng into feature/refactor-dpnp-backend-to-libtensor. No base build found for |
Add dpnp.random.MRG32k3a, L'Ecuyer's combined multiple recursive generator, on the mrg32k3a engine of the oneMKL device RNG API. It replaces the "MRG32k3a is not provided" skip of CuPy's TestBitGeneratorMRG32k3a. The engine gives values in [1, 2**32 - 209]. They are spread over the 32-bit range to make the raw words, which leaves 209 words that never occur, so the existing float and integer conversions apply unchanged. Work-item w of a launch draws from substream w, 2**76 values long, as the substreams of RngStreams. The host keeps the state at the start of the launch and moves it past the launch's substreams with cached matrix powers. The number of elements a work-item produces is now an engine trait, exported as the dict items_per_work_item. Philox4x3210 and PCG64DXSM keep 8, so their streams are unchanged. MRG32k3a uses 128: moving the engine to a substream takes a 3x3 matrix product per base-8 digit of the work-item index, which made 8 elements per work-item 2x (iGPU) to 10x (CPU) slower than Philox. With 128 it is on par with Philox on both. The tests check the kernel against a Python reference of the generator, which matches the stream of oneMKL's host MRG32K3A.
On devices without fp64 the discrete samplers compute in float32, and a host parameter that is valid in double could round onto the bound the kernel rejects: geometric(1e-300) rounded p to 0, zipf(1 + 1e-9) and logseries(1 - 1e-9) rounded a and p to 1, and poisson at NumPy's largest lam rounded it to 2**63. Every draw was then the kernel's value for an invalid parameter, INT64_MIN or 0. _discrete_sample now clips the host parameters of these samplers to float32 bounds inside the valid range when it computes in float32. The lower bound of geometric's p is the smallest normal float32, which also survives a device flushing subnormals. Parameters in device arrays are unchanged. Run test_edge_values with both discrete precisions and add test_param_near_bound for the four cases.
TestQueue named its device parameters by str(), which holds the object's address, so pytest-xdist workers collected different tests. Name them by the filter string. NumPy's legacy randint(dtype=int), random_integers and tomaxint return a C long, which is int32 on Windows, while dpnp returns int64. Compare their values without the dtype, compare tomaxint's moments with int64 samples where long is 32-bit, and document the difference in the docstrings and in random.rst. On devices without fp64 the discrete samplers compute in float32, so their values differ from the ones of fp64 devices. The device independence tests now compare discrete samples only between devices alike in fp64 support, and random.rst says so. The shuffle tests take their dtypes from get_all_dtypes() instead of hardcoding float64. Skip the CuPy tests that need fp64 on devices without it: the KS tests comparing with NumPy's float64 samples, the invalid out tests, which allocate float64 arrays or expect float32 to be the wrong type, and test_standard_gamma_generator.
BTPE computed the mode and every draw in floats, so in single precision the draws above 2**24 took only the values of floats: binomial(2**26, 0.5) never gave an odd number. Its final acceptance test summed terms of the order of |y - m| that cancel to a small result. The mode is now floor((n + 1) * r) computed exactly from the 128-bit product of n + 1 and the mantissa of r, each draw is the mode plus an offset computed in floats, and the final test uses log-factorial differences whose large terms are gathered into one log1p of an exact difference. HRUA's sum of four log-factorial differences cancelled the same way, as each is about |m - K| * log(popsize): in single precision the mean of hypergeometric(5e8, 4e8, 3e8) was 13 standard errors low at 2**20 draws. The large terms are now gathered into one log1p of the exact difference of two integer products. logfactorial_diff is split for this into the large term and logfactorial_diff_rest. zipf's rejection test computed (1 + 1/X)**(a - 1) - 1, which rounds to 0 for X above 2**24 in single precision and 2**53 in double, so every such candidate was accepted: the tail of zipf(1.1) above 2**25 was 30% too heavy in single precision. It is now expm1((a - 1) * log1p(1/X)). As NumPy's test has the same rounding above 2**53, its tail for a close to 1 is too heavy, while this now follows the Zipf distribution truncated to INT64_MAX; the docstrings and random.rst say so. In double precision, the binomial streams change only where the exact and the approximate final tests disagree, below 0.05% of the draws. test_large_params checks that draws above 2**25 are not limited to even numbers, and covers binomial(2**40, 0.5) and the large hypergeometric case. test_zipf_tail compares the tail with the exact truncated distribution. All of them fail without the fix. random.rst describes the float32 limit that remains: samples above 2**24 of geometric, logseries and zipf, and of poisson and binomial with a standard deviation above about 2**22, are multiples of a power of two.
Base automatically changed from
feature/refactor-dpnp-backend-to-libtensor
to
master
November 3, 2026 16:31
This branch has not been deployed
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.
This PR replaces the legacy
dpnp.random, which used the host-side oneMKL RNG through dpnp_backend_c, with a pybind11 extension that draws every sample on the device. It also adds NumPy's modern Generator API.Adds
dpnp.random.Generator,default_rng,BitGeneratorandSeedSequencerandom_raw,random(f32/f64),uniformandbytes)Generatormethods on deviceget/set_bit_generatorRandomStatealso receives changesRandomStateand usesdevice,usm_typeandsycl_queueare keyword-onlyrandnno longer takesd0=uniformno longer takes an integer dtyperandint(dtype=int)returns int64shufflerequires adpnp.ndarrayRemoves legacy RNG implementation as well as newly-dead code in the backend
Supesedes #1759