Add cuda headers for 6d arrays

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<T>` to `CArray4D<T>` 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<T>` to `CArray16D<T>` 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

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<T>` to `Array4D<T>`; a parameter or struct field of that type
+views `Array1D<T>` to `Array16D<T>`; 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<double> markers; long long n_markers; ...

docs/source/api.md

@@ -1117,8 +1117,8 @@ cunumpy ships CUDA headers that every `CudaKernel` finds automatically;
 (`-I<dir>`).
 
 `cunumpy/array_view.cuh` defines the strided views `Array1D<T>` to
-`Array4D<T>` (4D e.g. for a 3D grid of vector components `(nx, ny, nz,
-ncomp)`): `T* data`, `long long shape[ndim]`, `long long
+`Array16D<T>` (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<T>` to `Array4D<T>` of those scalar types (see "CUDA headers and
+`Array1D<T>` to `Array16D<T>` 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<double> 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<T>` to `Array4D<T>`) take arrays of the
+pointer and view parameters (`Array1D<T>` to `Array16D<T>`) 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,

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<T>` to `Array4D<T>`. Packing checks every field like a kernel
+`Array1D<T>` to `Array16D<T>`. 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<T>` 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<T>` to
-`CArray4D<T>` instead. The view holds a pointer and the shape, no strides, and
+`CArray16D<T>` 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.
 

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<double>` | 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<T>` ... `Array4D<T>` | CuPy array of dtype `T` and that ndim, contiguous or not | wrong dtype or ndim |
+| `Array1D<T>` ... `Array16D<T>` | 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/index.cuh>` | `CUNUMPY_THREAD_1D(i, n)`, `_2D`, `_3D`, `CUNUMPY_GRID_STRIDE_1D(i, n)` |
-| `<cunumpy/array_view.cuh>` | strided views `Array1D<T>` to `Array4D<T>` |
+| `<cunumpy/array_view.cuh>` | strided views `Array1D<T>` to `Array16D<T>` |
 | `<cunumpy/atomic.cuh>` | `cunumpy_atomic_add` and indexed 2D/3D variants, see [Accumulation kernels](accumulation.md) |
 | `<cunumpy/morton.cuh>` | Morton (Z-order) keys `cunumpy_morton_key2(x, y, ...)`, `_key3`, equal to `xp.algorithms.morton_keys` on the host |
 | `<cunumpy/random.cuh>` | 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)

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

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" ]

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/reduce.cuh>`: `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)` + `<cunumpy/atomic.cuh>` |
-| N-D indexing in CUDA, non-contiguous arrays | `Array1D<T>`..`Array4D<T>` params from `<cunumpy/array_view.cuh>` |
+| N-D indexing in CUDA, non-contiguous arrays | `Array1D<T>`..`Array16D<T>` params from `<cunumpy/array_view.cuh>` |
 | 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/index.cuh>       // CUNUMPY_THREAD_1D(i, n) /_2D/_3D, CUNUMPY_GRID_STRIDE_1D(i, n) {...}
-#include <cunumpy/array_view.cuh>  // Array1D<T>..Array4D<T>: data, shape[], strides[] (elements), a(i, j), size()
+#include <cunumpy/array_view.cuh>  // Array1D<T>..Array16D<T>: data, shape[], strides[] (elements), a(i, j), size()
 #include <cunumpy/atomic.cuh>      // cunumpy_atomic_add(double*|float*, v), _2d(data, n1, i, j, v), _3d(...)
 ```
 

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<T>``
-    to ``Array4D<T>`` views passed by value) and
+    to ``Array16D<T>`` 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<T>`` to
-        ``Array4D<T>`` or ``CArray1D<T>`` to ``CArray4D<T>``, see
+        ``Array16D<T>`` or ``CArray1D<T>`` to ``CArray16D<T>``, see
         :func:`cuda_include_dir`) passed by value.
     contiguous : bool
         Whether the array view is C-contiguous (``CArray2D<T>``): 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<T> to Array4D<T> and CArray1D<T> to CArray4D<T> (cunumpy/array_view.cuh),
+# Array1D<T> to Array16D<T> and CArray1D<T> to CArray16D<T> (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<T>`` to
-        ``Array4D<T>`` of those scalar types (from ``cunumpy/array_view.cuh``,
+        ``Array16D<T>`` 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<T>`` to
-        ``Array4D<T>`` views, ``#include "cunumpy/index.cuh"`` the thread-index
+        ``Array16D<T>`` 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

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<T>`` to ``Array4D<T>``; any strides, they are passed as
+parameters (``Array1D<T>`` to ``Array16D<T>``; 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
                 )