diff --git a/CHANGELOG.md b/CHANGELOG.md index 3273399..7022e86 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -44,7 +44,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 The former function names are removed without compatibility aliases. ### Added -- C-contiguous array views `CArray1D` to `CArray4D` in +- Strided and C-contiguous CUDA array views through 16 dimensions, including + `Array5D`/`Array6D` for matrix accumulations. Higher-dimensional views work + as kernel parameters, struct fields, annotations and in CPU emulation. +- C-contiguous array views `CArray1D` to `CArray16D` in `cunumpy/array_view.cuh`. They hold a pointer and shape only, so `a(i, j)` is `data[i * shape[1] + j]`. As kernel parameters or struct fields, they reject non-contiguous arrays and never copy them. `CudaStruct.from_signature` and diff --git a/README.md b/README.md index d84a081..21e5c34 100644 --- a/README.md +++ b/README.md @@ -394,7 +394,7 @@ kernel(particles.args_markers, dt) # host or CUDA kernel Kernels ported from pyccel index arrays like `markers[ip, j]`, which needs shapes and strides rather than bare pointers. The shipped header `cunumpy/array_view.cuh` (found by every `CudaKernel`) provides the strided -views `Array1D` to `Array4D`; a parameter or struct field of that type +views `Array1D` to `Array16D`; a parameter or struct field of that type takes a CuPy array, contiguous or not, and indexes `a(i, j)`. The struct can be generated from the annotations of the pyccel argument class, so the Python class is the one definition, and written to a header that a test keeps in sync: @@ -404,7 +404,9 @@ class MarkerArguments: def __init__(self, markers: "float[:, :]", n_markers: int, valid: "bool[:]"): ... -MarkerArgs = xp.arguments.CudaStruct.from_signature(MarkerArguments.__init__, "MarkerArgs") +MarkerArgs = xp.arguments.CudaStruct.from_signature( + MarkerArguments.__init__, "MarkerArgs" +) MarkerArgs.to_header( "marker_args.cuh" ) # Array2D markers; long long n_markers; ... diff --git a/docs/source/api.md b/docs/source/api.md index acfa555..06971ac 100644 --- a/docs/source/api.md +++ b/docs/source/api.md @@ -1117,8 +1117,8 @@ cunumpy ships CUDA headers that every `CudaKernel` finds automatically; (`-I`). `cunumpy/array_view.cuh` defines the strided views `Array1D` to -`Array4D` (4D e.g. for a 3D grid of vector components `(nx, ny, nz, -ncomp)`): `T* data`, `long long shape[ndim]`, `long long +`Array16D` (for example, a 4D view can describe a 3D grid of vector +components `(nx, ny, nz, ncomp)`): `T* data`, `long long shape[ndim]`, `long long strides[ndim]` (in elements, not bytes), `operator()(i, j, ...)` returning a reference to the element, and `size()`. A kernel indexes `a(i, j)` like the pyccel kernel it is ported from indexes `a[i, j]`, without hand-passed sizes. @@ -1261,7 +1261,7 @@ changes one definition instead of every kernel signature. `CudaStruct(name, fields)` takes the fields as `(name, C type)` pairs; scalar fields, pointers to the scalar types above (or `void*`), and array views -`Array1D` to `Array4D` of those scalar types (see "CUDA headers and +`Array1D` to `Array16D` of those scalar types (see "CUDA headers and array views") are supported. * `declaration`: the C definition of the struct, to put in the CUDA source @@ -1300,7 +1300,9 @@ class MarkerArguments: # the pyccel argument class, e.g. in struphy def __init__(self, markers: "float[:, :]", n_markers: int, valid: "bool[:]"): ... -MarkerArgs = xp.arguments.CudaStruct.from_signature(MarkerArguments.__init__, "MarkerArgs") +MarkerArgs = xp.arguments.CudaStruct.from_signature( + MarkerArguments.__init__, "MarkerArgs" +) print(MarkerArgs.declaration) # struct MarkerArgs { # Array2D markers; @@ -1897,7 +1899,7 @@ and the kernel's include directories and `-D` options apply. Then the kernel is called once per thread, for every block and thread index of the launch shape. Arguments follow the signature, with NumPy arrays in place of CuPy arrays: -pointer and view parameters (`Array1D` to `Array4D`) take arrays of the +pointer and view parameters (`Array1D` to `Array16D`) take arrays of the declared dtype (and ndim), passed as contiguous copies, so any strides work, and written back into the given arrays; scalars are checked and cast like in a launch. Like NVRTC by default, the compiler may fuse `a * b + c` into an FMA, diff --git a/docs/source/kernels/arguments.md b/docs/source/kernels/arguments.md index a243a71..6c53e39 100644 --- a/docs/source/kernels/arguments.md +++ b/docs/source/kernels/arguments.md @@ -80,7 +80,7 @@ push(value, 0.1, n_threads=x.size) ``` Fields may be scalars, pointers to scalar types (or `void*`), and array views -`Array1D` to `Array4D`. Packing checks every field like a kernel +`Array1D` to `Array16D`. Packing checks every field like a kernel argument: pointers need C-contiguous CuPy arrays of the declared dtype, scalars are range-checked and cast. Adding a field means editing the one Python definition; kernels that use the struct pick it up. @@ -191,7 +191,9 @@ class MarkerArguments: def __init__(self, markers: "float[:, :]", n_markers: int, valid: "bool[:]"): ... -MarkerArgs = xp.arguments.CudaStruct.from_signature(MarkerArguments.__init__, "MarkerArgs") +MarkerArgs = xp.arguments.CudaStruct.from_signature( + MarkerArguments.__init__, "MarkerArgs" +) print(MarkerArgs.declaration) ``` @@ -240,7 +242,7 @@ extern "C" __global__ void push(MarkerArgs m, double dt) { `Array2D` takes any view, so its strides are only known at run time and a kernel cannot tell which index is the fast one. When an array is always C-contiguous (a marker array, a grid), declare it as `CArray1D` to -`CArray4D` instead. The view holds a pointer and the shape, no strides, and +`CArray16D` instead. The view holds a pointer and the shape, no strides, and `m.markers(ip, 0)` is `data[ip * shape[1] + 0]`: the last index is always the fast one, as in the row-major memory the host code uses. diff --git a/docs/source/kernels/cuda-kernel.md b/docs/source/kernels/cuda-kernel.md index 7552ffa..e0b0949 100644 --- a/docs/source/kernels/cuda-kernel.md +++ b/docs/source/kernels/cuda-kernel.md @@ -54,7 +54,7 @@ every call: | `void*` | C-contiguous CuPy array of any dtype | host arrays, views | | `double`, `float`, `complex` | Python `int`/`float`, NumPy scalars that cast safely | strings, arrays, unsafe casts (`np.float64` into `float`) | | `int`, `long long`, `size_t`, `int64_t`, ... | Python `int` and NumPy integers whose value is in range, `bool` | out-of-range values (`OverflowError`), floats | -| `Array1D` ... `Array4D` | CuPy array of dtype `T` and that ndim, contiguous or not | wrong dtype or ndim | +| `Array1D` ... `Array16D` | CuPy array of dtype `T` and that ndim, contiguous or not | wrong dtype or ndim | | a `CudaStruct` type | a value of that struct | anything else | C types map to NumPy dtypes as on 64-bit Linux: `int` is `int32`, `long` and @@ -151,7 +151,9 @@ Keeping CUDA source in `.cu` files gives editor support and lets kernels share headers: ```python -push = xp.kernels.CudaKernel.from_file("kernels/push/push_cuda.cu") # kernel name "push" +push = xp.kernels.CudaKernel.from_file( + "kernels/push/push_cuda.cu" +) # kernel name "push" ``` `from_file` derives the kernel name from the file name minus the `_cuda.cu` @@ -177,7 +179,7 @@ Pass extra include directories with `include_dirs=[...]` and NVRTC flags with | Header | Provides | | --- | --- | | `` | `CUNUMPY_THREAD_1D(i, n)`, `_2D`, `_3D`, `CUNUMPY_GRID_STRIDE_1D(i, n)` | -| `` | strided views `Array1D` to `Array4D` | +| `` | strided views `Array1D` to `Array16D` | | `` | `cunumpy_atomic_add` and indexed 2D/3D variants, see [Accumulation kernels](accumulation.md) | | `` | Morton (Z-order) keys `cunumpy_morton_key2(x, y, ...)`, `_key3`, equal to `xp.algorithms.morton_keys` on the host | | `` | counter-based random numbers `cunumpy_uniform(seed, stream, counter)`, `cunumpy_normal2(...)`, equal to `xp.rng.philox_uniform` on the host | @@ -257,7 +259,9 @@ on first use and caches it: ```python def make_matvec(ndim, dtype): - return xp.kernels.CudaKernel(generate_source(ndim, xp.cuda.ctype_of(dtype)), "matvec") + return xp.kernels.CudaKernel( + generate_source(ndim, xp.cuda.ctype_of(dtype)), "matvec" + ) matvec = xp.kernels.CudaKernelVariants(make_matvec) diff --git a/docs/source/kernels/debugging.md b/docs/source/kernels/debugging.md index 5a2493e..d4c331e 100644 --- a/docs/source/kernels/debugging.md +++ b/docs/source/kernels/debugging.md @@ -28,7 +28,7 @@ In debug mode a `CudaKernel`: * is compiled with `-lineinfo` (source lines for `compute-sanitizer` and profilers) and `-DCUNUMPY_BOUNDS_CHECK`, which turns on bounds checks in - `Array1D` to `Array4D` views (an out-of-bounds index prints the index + `Array1D` to `Array16D` views (an out-of-bounds index prints the index and shape, then traps); * synchronizes after every launch, so a failure raises at the launch that caused it, as a `RuntimeError` naming the kernel and its grid and block, with diff --git a/pyproject.toml b/pyproject.toml index 329ab9b..7dbd7ee 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ requires = [ "setuptools", "wheel" ] [project] name = "cunumpy" -version = "0.6.0" +version = "0.6.1" description = "Simple wrapper for numpy and cupy. Replace `import numpy as np` with `import cunumpy as xp`." readme = "README.md" keywords = [ "python" ] diff --git a/src/cunumpy/LLM_GUIDE.md b/src/cunumpy/LLM_GUIDE.md index a72563d..65a6a8b 100644 --- a/src/cunumpy/LLM_GUIDE.md +++ b/src/cunumpy/LLM_GUIDE.md @@ -96,7 +96,7 @@ https://max-models.github.io/cunumpy/ and in `docs/source/` of the repository. | PETSc solve on device arrays without copies | `xp.petsc.petsc_vec(array)` (CUDA/HIP petsc4py for CuPy arrays); `xp.synchronize()` around PETSc calls | | reduction inside a CUDA kernel (energy, max velocity) | ``: `cunumpy_block_sum_to(out, v)`, `cunumpy_block_min/max`, `cunumpy_warp_sum` | | kernel writes into a host buffer owned by another library | `xp.memory.DeviceMirror(host_array)` + `` | -| N-D indexing in CUDA, non-contiguous arrays | `Array1D`..`Array4D` params from `` | +| N-D indexing in CUDA, non-contiguous arrays | `Array1D`..`Array16D` params from `` | | one MPI rank per GPU | `bind_local_device()` → `from mpi4py import MPI` → `require_cuda_aware_mpi()` → `synchronize_for_mpi(...)` before each call | | timing GPU code | `with xp.profiling.timed_region("name") as t:` → `t.elapsed` | | profiler markers | `xp.profiling.nvtx_range("name")` (context manager or decorator) | @@ -307,7 +307,7 @@ Shipped CUDA headers (always on the include path): ```c #include // CUNUMPY_THREAD_1D(i, n) /_2D/_3D, CUNUMPY_GRID_STRIDE_1D(i, n) {...} -#include // Array1D..Array4D: data, shape[], strides[] (elements), a(i, j), size() +#include // Array1D..Array16D: data, shape[], strides[] (elements), a(i, j), size() #include // cunumpy_atomic_add(double*|float*, v), _2d(data, n1, i, j, v), _3d(...) ``` diff --git a/src/cunumpy/_cuda_kernel.py b/src/cunumpy/_cuda_kernel.py index 1c691d1..fafd890 100644 --- a/src/cunumpy/_cuda_kernel.py +++ b/src/cunumpy/_cuda_kernel.py @@ -123,7 +123,7 @@ def cuda_include_dir() -> str: :class:`CudaKernel` adds it to the include path automatically, so kernels can ``#include "cunumpy/array_view.cuh"`` (strided ``Array1D`` - to ``Array4D`` views passed by value) and + to ``Array16D`` views passed by value) and ``#include "cunumpy/index.cuh"`` (thread-index and grid-stride macros such as ``CUNUMPY_THREAD_1D(i, n)``), ``#include "cunumpy/atomic.cuh"`` (atomic adds) and ``#include "cunumpy/reduce.cuh"`` (warp and block reductions). @@ -188,7 +188,7 @@ class CudaParameter(NamedTuple): The struct type, for a struct passed by value. view_ndim : int | None The number of dimensions, for an array view (``Array1D`` to - ``Array4D`` or ``CArray1D`` to ``CArray4D``, see + ``Array16D`` or ``CArray1D`` to ``CArray16D``, see :func:`cuda_include_dir`) passed by value. contiguous : bool Whether the array view is C-contiguous (``CArray2D``): packed @@ -265,11 +265,11 @@ class CudaParameter(NamedTuple): _QUALIFIERS = {"const", "volatile", "__restrict__", "__restrict", "restrict"} _COMPLEX = re.compile(r"(?:(?:thrust|cuda::std)::)?complex\s*<\s*(float|double)\s*>") -# Array1D to Array4D and CArray1D to CArray4D (cunumpy/array_view.cuh), +# Array1D to Array16D and CArray1D to CArray16D (cunumpy/array_view.cuh), # T a scalar type of _CTYPES -_VIEW = re.compile(r"\b(C?)Array([1234])D\s*<((?:[^<>]|complex<[^<>]*>)+?)>") +_VIEW = re.compile(r"\b(C?)Array([1-9]|1[0-6])D\s*<((?:[^<>]|complex<[^<>]*>)+?)>") _TOKEN = re.compile( - r"C?Array[1234]D<[^<>]*(?:<[^<>]*>[^<>]*)?>|complex<(?:float|double)>" + r"C?Array(?:[1-9]|1[0-6])D<[^<>]*(?:<[^<>]*>[^<>]*)?>|complex<(?:float|double)>" r"|[A-Za-z_]\w*|\*|\[\s*\]", ) @@ -905,8 +905,8 @@ def _pyccel_ctype(annotation: Any, scalars: Mapping[str, str], what: str) -> str ctype = scalars[scalar] if ndim == 0: return ctype - if ndim > 4: - raise ValueError(f"{what}: arrays have at most 4 dimensions, got {ndim}") + if ndim > 16: + raise ValueError(f"{what}: arrays have at most 16 dimensions, got {ndim}") return f"Array{ndim}D<{ctype}>" @@ -1075,7 +1075,7 @@ class CudaStruct: ``(field name, C type)`` pairs, in order, e.g. ``("x", "double*")`` or ``("n", "int")``. Scalar fields, pointers to the scalar types of :func:`ctype_of` (or ``void*``), and array views ``Array1D`` to - ``Array4D`` of those scalar types (from ``cunumpy/array_view.cuh``, + ``Array16D`` of those scalar types (from ``cunumpy/array_view.cuh``, packed as pointer, shape and strides in elements) are supported. Examples @@ -1848,7 +1848,7 @@ class CudaKernel: The headers shipped with cunumpy (:func:`cuda_include_dir`) are always found at compile time (see :meth:`compile_options`): ``#include "cunumpy/array_view.cuh"`` gives the ``Array1D`` to - ``Array4D`` views, ``#include "cunumpy/index.cuh"`` the thread-index + ``Array16D`` views, ``#include "cunumpy/index.cuh"`` the thread-index macros, ``#include "cunumpy/atomic.cuh"`` atomic adds, ``#include "cunumpy/reduce.cuh"`` warp and block reductions. source_dir : str | Path | None diff --git a/src/cunumpy/_emulation.py b/src/cunumpy/_emulation.py index d3e2de4..56d47d0 100644 --- a/src/cunumpy/_emulation.py +++ b/src/cunumpy/_emulation.py @@ -16,7 +16,7 @@ np.testing.assert_allclose(y, 2.0 * x) Arguments follow the kernel signature: NumPy arrays for pointer and array view -parameters (``Array1D`` to ``Array4D``; any strides, they are passed as +parameters (``Array1D`` to ``Array16D``; any strides, they are passed as contiguous copies), Python or NumPy scalars for scalar parameters (cast and checked like in a launch). Arrays are written back into the given arrays. @@ -354,7 +354,7 @@ def emulate_cuda_kernel( buffer.tofile(path) ctype = param.ctype if param.dtype is not None else "unsigned char" element = ( - re.match(r"C?Array\dD<(.*)>", ctype).group(1) + re.match(r"C?Array\d+D<(.*)>", ctype).group(1) if param.view_ndim is not None else ctype ) diff --git a/src/cunumpy/cuda/include/cunumpy/array_view.cuh b/src/cunumpy/cuda/include/cunumpy/array_view.cuh index d93667f..2a299ae 100644 --- a/src/cunumpy/cuda/include/cunumpy/array_view.cuh +++ b/src/cunumpy/cuda/include/cunumpy/array_view.cuh @@ -1,6 +1,6 @@ // Array views for CUDA kernels, passed by value from Python. // -// Array1D to Array4D describe a (possibly non-contiguous) +// Array1D to Array16D describe a (possibly non-contiguous) // device array the way NumPy/CuPy do: a data pointer, a shape and strides. // Strides are in ELEMENTS, not bytes, so that `a(i, j)` is // `data[i * strides[0] + j * strides[1]]`. Elements are accessed with @@ -9,7 +9,7 @@ // // The views are created on the Python side by cunumpy (a kernel parameter or a // CudaStruct field of type `Array2D` takes a CuPy array). They are -// passed by value, so the memory layout must be exactly, for ndim = 1 to 4: +// passed by value, so the memory layout must be exactly, for ndim = 1 to 16: // // T* data; // 8 bytes // long long shape[ndim]; // ndim * 8 bytes @@ -20,12 +20,12 @@ // the end of this file check the size. Member functions do not change the // layout. // -// CArray1D to CArray4D are C-contiguous (row-major) views: a data +// CArray1D to CArray16D are C-contiguous (row-major) views: a data // pointer and a shape, no strides. `a(i, j)` is `data[i * shape[1] + j]`, so // the last index is always the fast one and the compiler knows it has unit // stride. A parameter or field of these types only takes C-contiguous arrays; // cunumpy raises for a non-contiguous view instead of copying it (a copy would -// silently drop what the kernel writes). Their layout, for ndim = 1 to 4: +// silently drop what the kernel writes). Their layout, for ndim = 1 to 16: // // T* data; // 8 bytes // long long shape[ndim]; // ndim * 8 bytes @@ -223,6 +223,101 @@ struct CArray4D { } }; +// Higher-dimensional views share an implementation; aliases keep the named +// types used by kernel signatures and Python's argument packer. +namespace cunumpy_detail { +template +struct ArrayView { + T* data; + long long shape[N]; + long long strides[N]; + + template + __device__ __forceinline__ T& operator()(Indices... indices) const { + static_assert(sizeof...(Indices) == N, "wrong number of array indices"); + const long long index[N] = {static_cast(indices)...}; + long long offset = 0; + #pragma unroll + for (int axis = 0; axis < N; ++axis) { + CUNUMPY_CHECK_INDEX(index[axis], axis, shape[axis]); + offset += index[axis] * strides[axis]; + } + return data[offset]; + } + + __device__ __forceinline__ long long size() const { + long long result = 1; + #pragma unroll + for (int axis = 0; axis < N; ++axis) result *= shape[axis]; + return result; + } +}; + +template +struct CArrayView { + T* data; + long long shape[N]; + + template + __device__ __forceinline__ T& operator()(Indices... indices) const { + static_assert(sizeof...(Indices) == N, "wrong number of array indices"); + const long long index[N] = {static_cast(indices)...}; + long long offset = 0; + #pragma unroll + for (int axis = 0; axis < N; ++axis) { + CUNUMPY_CHECK_INDEX(index[axis], axis, shape[axis]); + offset = offset * shape[axis] + index[axis]; + } + return data[offset]; + } + + __device__ __forceinline__ long long size() const { + long long result = 1; + #pragma unroll + for (int axis = 0; axis < N; ++axis) result *= shape[axis]; + return result; + } + + __host__ __device__ operator ArrayView() const { + ArrayView result{}; + result.data = data; + long long stride = 1; + #pragma unroll + for (int axis = N - 1; axis >= 0; --axis) { + result.shape[axis] = shape[axis]; + result.strides[axis] = stride; + stride *= shape[axis]; + } + return result; + } +}; +} // namespace cunumpy_detail + +template using Array5D = cunumpy_detail::ArrayView; +template using CArray5D = cunumpy_detail::CArrayView; +template using Array6D = cunumpy_detail::ArrayView; +template using CArray6D = cunumpy_detail::CArrayView; +template using Array7D = cunumpy_detail::ArrayView; +template using CArray7D = cunumpy_detail::CArrayView; +template using Array8D = cunumpy_detail::ArrayView; +template using CArray8D = cunumpy_detail::CArrayView; +template using Array9D = cunumpy_detail::ArrayView; +template using CArray9D = cunumpy_detail::CArrayView; +template using Array10D = cunumpy_detail::ArrayView; +template using CArray10D = cunumpy_detail::CArrayView; +template using Array11D = cunumpy_detail::ArrayView; +template using CArray11D = cunumpy_detail::CArrayView; +template using Array12D = cunumpy_detail::ArrayView; +template using CArray12D = cunumpy_detail::CArrayView; +template using Array13D = cunumpy_detail::ArrayView; +template using CArray13D = cunumpy_detail::CArrayView; +template using Array14D = cunumpy_detail::ArrayView; +template using CArray14D = cunumpy_detail::CArrayView; +template using Array15D = cunumpy_detail::ArrayView; +template using CArray15D = cunumpy_detail::CArrayView; +template using Array16D = cunumpy_detail::ArrayView; +template using CArray16D = cunumpy_detail::CArrayView; + // The layouts the Python side packs: pointer, shape (and strides), 8-byte // aligned. static_assert(sizeof(Array1D) == 24, "unexpected Array1D layout"); @@ -237,4 +332,53 @@ static_assert(sizeof(CArray3D) == 32, "unexpected CArray3D layout"); static_assert(sizeof(CArray4D) == 40, "unexpected CArray4D layout"); static_assert(alignof(CArray2D) == 8, "unexpected CArray2D alignment"); +static_assert(sizeof(Array5D) == 88, "unexpected Array5D layout"); +static_assert(alignof(Array5D) == 8, "unexpected Array5D alignment"); +static_assert(sizeof(CArray5D) == 48, "unexpected CArray5D layout"); +static_assert(alignof(CArray5D) == 8, "unexpected CArray5D alignment"); +static_assert(sizeof(Array6D) == 104, "unexpected Array6D layout"); +static_assert(alignof(Array6D) == 8, "unexpected Array6D alignment"); +static_assert(sizeof(CArray6D) == 56, "unexpected CArray6D layout"); +static_assert(alignof(CArray6D) == 8, "unexpected CArray6D alignment"); +static_assert(sizeof(Array7D) == 120, "unexpected Array7D layout"); +static_assert(alignof(Array7D) == 8, "unexpected Array7D alignment"); +static_assert(sizeof(CArray7D) == 64, "unexpected CArray7D layout"); +static_assert(alignof(CArray7D) == 8, "unexpected CArray7D alignment"); +static_assert(sizeof(Array8D) == 136, "unexpected Array8D layout"); +static_assert(alignof(Array8D) == 8, "unexpected Array8D alignment"); +static_assert(sizeof(CArray8D) == 72, "unexpected CArray8D layout"); +static_assert(alignof(CArray8D) == 8, "unexpected CArray8D alignment"); +static_assert(sizeof(Array9D) == 152, "unexpected Array9D layout"); +static_assert(alignof(Array9D) == 8, "unexpected Array9D alignment"); +static_assert(sizeof(CArray9D) == 80, "unexpected CArray9D layout"); +static_assert(alignof(CArray9D) == 8, "unexpected CArray9D alignment"); +static_assert(sizeof(Array10D) == 168, "unexpected Array10D layout"); +static_assert(alignof(Array10D) == 8, "unexpected Array10D alignment"); +static_assert(sizeof(CArray10D) == 88, "unexpected CArray10D layout"); +static_assert(alignof(CArray10D) == 8, "unexpected CArray10D alignment"); +static_assert(sizeof(Array11D) == 184, "unexpected Array11D layout"); +static_assert(alignof(Array11D) == 8, "unexpected Array11D alignment"); +static_assert(sizeof(CArray11D) == 96, "unexpected CArray11D layout"); +static_assert(alignof(CArray11D) == 8, "unexpected CArray11D alignment"); +static_assert(sizeof(Array12D) == 200, "unexpected Array12D layout"); +static_assert(alignof(Array12D) == 8, "unexpected Array12D alignment"); +static_assert(sizeof(CArray12D) == 104, "unexpected CArray12D layout"); +static_assert(alignof(CArray12D) == 8, "unexpected CArray12D alignment"); +static_assert(sizeof(Array13D) == 216, "unexpected Array13D layout"); +static_assert(alignof(Array13D) == 8, "unexpected Array13D alignment"); +static_assert(sizeof(CArray13D) == 112, "unexpected CArray13D layout"); +static_assert(alignof(CArray13D) == 8, "unexpected CArray13D alignment"); +static_assert(sizeof(Array14D) == 232, "unexpected Array14D layout"); +static_assert(alignof(Array14D) == 8, "unexpected Array14D alignment"); +static_assert(sizeof(CArray14D) == 120, "unexpected CArray14D layout"); +static_assert(alignof(CArray14D) == 8, "unexpected CArray14D alignment"); +static_assert(sizeof(Array15D) == 248, "unexpected Array15D layout"); +static_assert(alignof(Array15D) == 8, "unexpected Array15D alignment"); +static_assert(sizeof(CArray15D) == 128, "unexpected CArray15D layout"); +static_assert(alignof(CArray15D) == 8, "unexpected CArray15D alignment"); +static_assert(sizeof(Array16D) == 264, "unexpected Array16D layout"); +static_assert(alignof(Array16D) == 8, "unexpected Array16D alignment"); +static_assert(sizeof(CArray16D) == 136, "unexpected CArray16D layout"); +static_assert(alignof(CArray16D) == 8, "unexpected CArray16D alignment"); + #endif // CUNUMPY_ARRAY_VIEW_CUH diff --git a/tests/unit/test_cuda_kernel.py b/tests/unit/test_cuda_kernel.py index 011fcf2..0ef8c8c 100644 --- a/tests/unit/test_cuda_kernel.py +++ b/tests/unit/test_cuda_kernel.py @@ -705,7 +705,7 @@ def test_parse_view_parameters(): with pytest.raises(ValueError, match="array views"): parse_cuda_signature("__global__ void f(Array2D a) {}", "f") with pytest.raises(ValueError, match="unsupported type"): - parse_cuda_signature("__global__ void f(Array5D a) {}", "f") + parse_cuda_signature("__global__ void f(Array17D a) {}", "f") def test_view_parameters_pack_pointer_shape_and_strides(): @@ -1088,7 +1088,7 @@ def missing(x, n: int): def unknown(x: "str[:]"): pass - def too_many(x: "float[:, :, :, :, :]"): + def too_many(x: "float[:, :, :, :, :, :, :, :, :, :, :, :, :, :, :, :, :]"): pass def unparsable(x: "float[:](order=F)"): # noqa: F821 (deliberately unparsable) @@ -1101,7 +1101,7 @@ def unparsable(x: "float[:](order=F)"): # noqa: F821 (deliberately unparsable) CudaStruct.from_signature(missing, "A") with pytest.raises(ValueError, match="unsupported scalar type 'str'"): CudaStruct.from_signature(unknown, "A") - with pytest.raises(ValueError, match="at most 4 dimensions"): + with pytest.raises(ValueError, match="at most 16 dimensions"): CudaStruct.from_signature(too_many, "A") with pytest.raises(ValueError, match="cannot parse the annotation"): CudaStruct.from_signature(unparsable, "A") @@ -2263,9 +2263,11 @@ def init(self, e: "float[:, :, :, :]", n: int): ... from_annotations = CudaStruct.from_signature(init, "Grid") assert from_annotations.fields[0].ctype == "Array4D" - def too_many(self, e: "float[:, :, :, :, :]"): ... + def too_many( + self, e: "float[:, :, :, :, :, :, :, :, :, :, :, :, :, :, :, :, :]" + ): ... - with pytest.raises(ValueError, match="at most 4 dimensions"): + with pytest.raises(ValueError, match="at most 16 dimensions"): CudaStruct.from_signature(too_many, "Grid") @@ -2381,3 +2383,78 @@ def test_shared_memory_above_the_default_is_opted_in(recorded): with pytest.raises(ValueError, match="exceeds the 100000 bytes"): kernel(1.0, x, y, 1, n_threads=1, shared_mem=100_001) assert [s for *_, s in raw.launches] == [40_000, 80_000, 60_000] + + +@pytest.mark.parametrize("ndim", range(5, 17)) +@pytest.mark.parametrize("contiguous", [False, True]) +def test_high_dimensional_views_pack_and_generate_structs(ndim, contiguous): + ctype = f"{'C' if contiguous else ''}Array{ndim}D" + source = f"__global__ void f({ctype} a) {{}}" + (param,) = parse_cuda_signature(source, "f") + assert (param.view_ndim, param.contiguous) == (ndim, contiguous) + shape = (2, 3) + (1,) * (ndim - 3) + (4,) + array = FakeDeviceArray(np.float64, shape=shape) + if not contiguous: + array.strides = tuple(-2 * s for s in array.strides) + (packed,) = CudaKernel(source, "f").prepare_args(array) + assert packed["data"] == array.data.ptr + assert packed["shape"].tolist() == list(shape) + assert packed.dtype.itemsize == 8 * (1 + ndim * (1 if contiguous else 2)) + if not contiguous: + assert packed["strides"].tolist() == [s // 8 for s in array.strides] + else: + assert packed.dtype.names == ("data", "shape") + with pytest.raises(TypeError, match="must be C-contiguous"): + CudaKernel(source, "f").prepare_args( + FakeDeviceArray(np.float64, shape=shape, strides=(16,) * ndim) + ) + with pytest.raises(TypeError, match=f"must be a {ndim}D array"): + CudaKernel(source, "f").prepare_args(FakeDeviceArray(np.float64)) + with pytest.raises(TypeError, match="dtype"): + CudaKernel(source, "f").prepare_args(FakeDeviceArray(np.float32, shape=shape)) + + def init(a): ... + + init.__annotations__ = {"a": "float[" + ", ".join([":"] * ndim) + "]"} + struct = CudaStruct.from_signature(init, "HighDim", contiguous=contiguous) + assert struct.fields[0].ctype == ctype + assert struct.dtype.fields["a"][0] == packed.dtype + assert f"{ctype} a;" in struct.declaration + value = struct(a=array) + assert value.packed["a"]["shape"].tolist() == list(shape) + + +@pytest.mark.parametrize("ndim", [5, 6, 10, 16]) +@pytest.mark.parametrize("contiguous", [False, True]) +def test_high_dimensional_struct_layout_and_execution_on_gpu(ndim, contiguous): + _skip_without_cupy() + import cupy as cp + + annotation = "float[" + ", ".join([":"] * ndim) + "]" + struct = CudaStruct.from_pyccel_class( + f'class Grid:\n def __init__(self, a: "{annotation}", n: int):\n' + " self.a = a\n self.n = n\n", + "Grid", + contiguous=contiguous, + ) + struct.verify_layout() + indices = ", ".join(["0"] * (ndim - 1) + ["i"]) + source = ( + struct.to_header() + + f""" + #include + #include + extern "C" __global__ void accumulate(Grid g) {{ + CUNUMPY_THREAD_1D(i, g.n); + cunumpy_atomic_add(&g.a({indices}), 2.0); + }} + """ + ) + base = cp.zeros((1,) * (ndim - 1) + (6,)) + a = base if contiguous else base[..., ::2] + CudaKernel(source, "accumulate", structs=(struct,))( + struct(a=a, n=a.size), n_threads=a.size + ) + expected = np.zeros(base.shape) + expected[..., slice(None) if contiguous else slice(None, None, 2)] = 2 + np.testing.assert_array_equal(cp.asnumpy(base), expected) diff --git a/tests/unit/test_emulation.py b/tests/unit/test_emulation.py index f9e5482..c5f28f6 100644 --- a/tests/unit/test_emulation.py +++ b/tests/unit/test_emulation.py @@ -387,3 +387,72 @@ def test_bounds_checks_from_the_view_header(): n_threads=4, options=("-DCUNUMPY_BOUNDS_CHECK",), ) + + +@pytest.mark.parametrize("ndim", range(5, 17)) +@pytest.mark.parametrize("contiguous", [False, True]) +@pytest.mark.parametrize("backend", ["emulation", "cuda"]) +def test_high_dimensional_view_indexing_and_conversion(ndim, contiguous, backend): + # Exercise every axis with unequal extents, including singleton dimensions. + # Reading via a strided device helper also tests CArray -> Array conversion. + ctype = f"{'C' if contiguous else ''}Array{ndim}D" + indices = ", ".join(f"i[{axis}]" for axis in range(ndim)) + source = f""" + #include + #include + __device__ double read(Array{ndim}D a, const long long* i) {{ + return a({indices}); + }} + extern "C" __global__ void update({ctype} a) {{ + CUNUMPY_THREAD_1D(flat, a.size()); + long long i[{ndim}], remaining = flat; + for (int axis = {ndim} - 1; axis >= 0; --axis) {{ + i[axis] = remaining % a.shape[axis]; + remaining /= a.shape[axis]; + }} + a({indices}) = read(a, i) + flat + 1; + }} + """ + shape = (2, 3) + (1,) * (ndim - 3) + (4,) + base = np.arange(48.0).reshape(shape[:-1] + (8,)) + if contiguous: + a = base[..., :4].copy() + else: + a = base[..., ::-2].swapaxes(0, 1) + expected = a.copy() + np.arange(a.size).reshape(a.shape) + 1 + kernel = CudaKernel(source, "update", options=("-DCUNUMPY_BOUNDS_CHECK",)) + if backend == "emulation": + emulate_cuda_kernel(kernel, a, n_threads=a.size) + np.testing.assert_array_equal(a, expected) + if not contiguous: + np.testing.assert_array_equal( + base[..., 0::2], np.arange(48.0).reshape(base.shape)[..., 0::2] + ) + else: + import cunumpy as xp + + if not xp.cupy_available(): + pytest.skip("CuPy not installed or not functional") + import cupy as cp + + device_base = cp.asarray(base) + device = cp.asarray(a) if contiguous else device_base[..., ::-2].swapaxes(0, 1) + kernel(device, n_threads=device.size) + np.testing.assert_array_equal(cp.asnumpy(device), expected) + + +@pytest.mark.parametrize("contiguous", [False, True]) +def test_high_dimensional_bounds_check(contiguous): + ctype = f"{'C' if contiguous else ''}Array16D" + indices = ", ".join(["0"] * 15 + ["1"]) + source = f""" + #include + extern "C" __global__ void invalid({ctype} a) {{ a({indices}) = 1; }} + """ + with pytest.raises(RuntimeError, match="crashed"): + emulate_cuda_kernel( + CudaKernel(source, "invalid"), + np.zeros((1,) * 16), + n_threads=1, + options=("-DCUNUMPY_BOUNDS_CHECK",), + )