Version 0.4.0

.github/workflows/gpu_ci_trigger.yml

@@ -64,6 +64,7 @@ jobs:
 
           # 4. Force push (This automatically starts the GitLab Pipeline)
           git push -f gitlab HEAD:refs/heads/$TARGET_BRANCH
+          echo "PUSHED_SHA=$(git rev-parse HEAD)" >> $GITHUB_ENV
 
           # 5. Provide the direct link
           PIPELINE_URL="https://gitlab.mpcdf.mpg.de/maxlin/cunumpy/-/pipelines?ref=$TARGET_BRANCH"
@@ -72,10 +73,41 @@ jobs:
           echo "::notice::View Pipeline: $PIPELINE_URL"
 
       - name: Wait for GitLab Pipeline
-        uses: docker://gitlab/glab:latest
+        timeout-minutes: 90
         env:
           GITLAB_TOKEN: ${{ secrets.GITLAB_TOKEN }}
-          GITLAB_HOST: gitlab.mpcdf.mpg.de
-        with:
-          entrypoint: glab
-          args: ci status --live --branch ${{ env.TARGET_BRANCH }} --repo maxlin/cunumpy
+        run: |
+          API="https://gitlab.mpcdf.mpg.de/api/v4/projects/maxlin%2Fcunumpy"
+          AUTH=()
+          if [ -n "$GITLAB_TOKEN" ]; then AUTH=(--header "PRIVATE-TOKEN: $GITLAB_TOKEN"); fi
+
+          # 1. GitLab creates the pipeline asynchronously after the push: wait until
+          #    the pipeline for the pushed commit exists (asking right away finds none).
+          PIPELINE=""
+          for i in $(seq 1 60); do
+            PIPELINE=$(curl -sf "${AUTH[@]}" "$API/pipelines?ref=$TARGET_BRANCH&sha=$PUSHED_SHA&per_page=1" | jq -r '.[0].id // empty' || true)
+            if [ -n "$PIPELINE" ]; then break; fi
+            sleep 5
+          done
+          if [ -z "$PIPELINE" ]; then
+            echo "::error::No GitLab pipeline for $PUSHED_SHA on $TARGET_BRANCH after 5 minutes"
+            exit 1
+          fi
+          URL="https://gitlab.mpcdf.mpg.de/maxlin/cunumpy/-/pipelines/$PIPELINE"
+          echo "::notice::GitLab pipeline: $URL"
+
+          # 2. Follow the pipeline until it has finished.
+          while true; do
+            STATUS=$(curl -sf "${AUTH[@]}" "$API/pipelines/$PIPELINE" | jq -r '.status // empty' || true)
+            case "$STATUS" in
+              success)
+                echo "GitLab pipeline $PIPELINE succeeded: $URL"
+                exit 0 ;;
+              failed|canceled|skipped|manual|scheduled)
+                echo "::error::GitLab pipeline $PIPELINE is $STATUS: $URL"
+                exit 1 ;;
+              *)
+                echo "GitLab pipeline $PIPELINE: ${STATUS:-status not available yet}"
+                sleep 20 ;;
+            esac
+          done

.github/workflows/testing.yml

@@ -17,7 +17,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.10", "3.13"]
+        python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
 
     steps:
       # Checkout the repository

.gitlab-ci.yml

@@ -29,18 +29,19 @@ gpu_tests:
     - git --version
 
     - echo "--- Pytest Execution ---"
-    # The MPCDF image likely has a specific python environment. 
+    # The MPCDF image likely has a specific python environment.
     # We install our dependencies into the user directory or a virtualenv.
-    - python3 -m pip install --user cupy-cuda12x
-    - python3 -m pip install --user nvidia-cublas-cu12 nvidia-cufft-cu12 nvidia-curand-cu12 nvidia-cusolver-cu12 nvidia-cusparse-cu12
+    # One CUDA version only: nvhpcsdk/26 provides CUDA 13.2 (its headers are used
+    # when CuPy compiles kernels with NVRTC), so install CuPy for CUDA 13 with the
+    # CUDA 13.2 libraries and NVRTC. Mixing CUDA 12 (cupy-cuda12x) with the CUDA 13.2
+    # headers fails to compile CuPy's own kernels (e.g. CUB reductions).
+    - python3 -m pip install --user "cupy-cuda13x[ctk]" "cuda-toolkit==13.2.*"
     - python3 -m pip install --user -e .
-    
+
     # Add the user bin to PATH for pytest
     - export PATH="$HOME/.local/bin:$PATH"
-
-    # Try to find libcublas and other libraries in the HPC environment
-    - export LD_LIBRARY_PATH=$(find /mpcdf/soft /opt/nvidia -name libcublas.so.12 -exec dirname {} \; 2>/dev/null | head -n 1):$LD_LIBRARY_PATH
-
+
     - export ARRAY_BACKEND=cupy
+    - python3 -c "import cupy; cupy.show_config()"
 
     - pytest -xvs .

README.md

@@ -123,6 +123,21 @@ with xp.use_backend("cupy"):
     result = xp.to_numpy(filtered)  # one transfer for a CPU-only consumer
 ```
 
+To verify that a block, such as a time step, makes no transfer at all, count
+them: `count_transfers()` records every `to_numpy()`, `to_cupy()` and
+`to_cunumpy()` call that actually copies, every `PyccelKernel` call that
+converts device arrays, and every `Kernel` fallback to the host kernel, with
+the call site of each. `assert_no_transfers()` raises with that report if
+anything was counted. Only transfers made through CuNumpy are seen; raw
+`cupy.ndarray.get()` or `cupy.asarray()` calls need a profiler such as `nsys`.
+
+```python
+with xp.count_transfers() as counter:
+    propagator(dt)
+
+assert counter.total == 0, counter.report()
+```
+
 ## Random numbers and dtypes
 
 `get_rng(seed)` returns a random generator for the active backend. NumPy and
@@ -159,6 +174,34 @@ active CuPy device and `None` on NumPy. `set_device_for_rank(rank)` is a
 round-robin convenience for MPI layouts where local ranks map contiguously to
 GPUs. If your scheduler uses a different mapping, select the device directly.
 
+For MPI programs with one rank per GPU, the startup sequence is:
+
+1. `bind_local_device()` selects the GPU from the node-local rank that the MPI
+   launcher exports (`local_rank()`) and creates its CUDA context. It runs
+   before MPI is initialized because a CUDA-aware MPI binds to the device that
+   is current at `MPI_Init`; without it, every rank of a node would use
+   device 0.
+2. `from mpi4py import MPI` initializes MPI.
+3. `require_cuda_aware_mpi()` (or `mpi_is_cuda_aware(comm)`) checks, with one
+   tiny device `Sendrecv` on every rank, that the MPI library can pass device
+   buffers at all. Passing CuPy arrays to a plain MPI build segfaults or
+   silently sends garbage; the check turns that into a clear error at
+   startup. It is a no-op on the NumPy backend.
+4. `synchronize_for_mpi(*buffers)` before every MPI call with device buffers:
+   kernels run asynchronously, and MPI would otherwise send a buffer a kernel
+   is still writing, without an error.
+
+```python
+xp.set_backend("cupy")
+xp.bind_local_device()  # before MPI_Init
+from mpi4py import MPI  # MPI_Init
+
+xp.require_cuda_aware_mpi()  # once, on all ranks
+
+xp.synchronize_for_mpi(send, recv)
+MPI.COMM_WORLD.Sendrecv(send, dest, recvbuf=recv, source=source)
+```
+
 CuPy caches released allocations in memory pools. This can make process-level
 GPU memory appear occupied after arrays go out of scope. `free_memory()` asks
 CuPy to release currently free cached blocks; it does not free memory still
@@ -179,6 +222,23 @@ xp.synchronize()
 result = xp.to_numpy(transformed)
 ```
 
+Because GPU work is asynchronous, a wall-clock timer around a kernel launch
+measures the launch, not the kernel. `timed_region(name)` synchronizes the
+device before reading the clock (on NumPy it is a plain timer), and
+`nvtx_range(name)` marks a region so it shows up in `nsys`/Nsight; both are
+no-ops or plain timers on NumPy, and `nvtx_range` also works as a decorator:
+
+```python
+with xp.timed_region("fft") as timing:
+    transformed = xp.fft.fft(device)
+print(timing.elapsed, timing.synced)
+
+
+@xp.nvtx_range("step")
+def step(dt):
+    ...
+```
+
 ## Use NumPy-only kernels with CuPy arrays
 
 `PyccelKernel` adapts a callable that expects NumPy arrays. When conversion is
@@ -213,6 +273,197 @@ The wrapper can also traverse arrays nested in lists, tuples, dictionaries,
 and selected application objects; see the full [API reference](docs/source/api.md)
 for `object_modules`, `is_array`, aliasing, and output declarations.
 
+## Write CUDA kernels next to host kernels
+
+`CudaKernel` wraps a CUDA C kernel (compiled with NVRTC through
+`cupy.RawKernel`) so that it is called with the same arguments as the host
+kernel it mirrors, plus the number of threads. Arrays are never copied: they
+must be C-contiguous CuPy arrays. The `extern "C" __global__` signature is
+parsed once and every call is checked against it: Python scalars are cast to
+the declared C types, and a wrong argument count, an array of the wrong dtype
+or a non-contiguous view, or a scalar that does not fit its type raises instead
+of silently producing wrong values.
+
+`Kernel` pairs a host kernel with its CUDA kernel and calls the one matching
+the active backend, so kernels can be ported to CUDA one at a time:
+
+```python
+import cunumpy as xp
+
+AXPY = r"""
+extern "C" __global__
+void axpy(double a, const double* x, double* y, int n) {
+    int i = blockDim.x * blockIdx.x + threadIdx.x;
+    if (i < n) y[i] += a * x[i];
+}
+"""
+
+
+def axpy(a, x, y, n):  # host version, e.g. compiled with Pyccel
+    for i in range(n):
+        y[i] += a * x[i]
+
+
+kernel = xp.Kernel(axpy, xp.CudaKernel(AXPY, "axpy"))
+
+with xp.use_backend("cupy"):
+    x = xp.arange(1000, dtype=xp.float64)
+    y = xp.zeros(1000)
+    kernel(2.0, x, y, 1000, n_threads=1000)  # runs the CUDA kernel
+```
+
+On the CuPy backend, a `Kernel` without CUDA kernel raises
+`NotImplementedError` (or, with `missing_cuda="fallback"`, runs the host kernel
+through `PyccelKernel`, with host copies; `host_options` configure that
+`PyccelKernel`). `KernelCatalog.from_package()` collects kernel pairs from a
+package with one folder per kernel (`name/name_kernels.py` and
+`name/name_cuda.cu`), and `catalog.compile_all()` compiles all CUDA kernels at
+setup.
+
+Groups of arguments can be passed as one: objects implementing
+`__cuda_args__()` (see `CudaArguments`) are flattened into several kernel
+arguments, and `CudaStruct` defines a C struct once (its C `declaration` and
+the matching memory layout) and packs values into it, which the kernel takes
+as one parameter:
+
+```python
+Vec = xp.CudaStruct("Vec", [("data", "double*"), ("n", "int")])
+scale = xp.CudaKernel(
+    Vec.declaration
+    + r"""
+    extern "C" __global__ void scale(Vec v, double a) {
+        int i = blockDim.x * blockIdx.x + threadIdx.x;
+        if (i < v.n) v.data[i] *= a;
+    }""",
+    "scale",
+    structs=[Vec],
+)
+scale(Vec(data=y, n=y.size), 0.5, n_threads=y.size)
+```
+
+When building such argument objects, `xp.as_device_array(value, dtype,
+ndim=None)` applies the "reference or copy once" rule: a CuPy array that
+already has the dtype and is C-contiguous is returned as it is, anything else
+(a tuple such as `degree = (3, 3, 3)`, a host array, another dtype, a
+non-contiguous view) becomes one device copy. Call it once when the object is
+built, not per kernel call; on the NumPy backend it raises, so host data is
+never copied to the device implicitly.
+When the host kernel takes such a group as one object too (e.g. a Pyccel class
+holding NumPy arrays), give the group both forms with `KernelArguments`:
+`__host_args__()` returns the object for the host kernel, `__cuda_args__()`
+the flattened device arguments. `Kernel` and `PyccelKernel` resolve
+`__host_args__()` on the host path and `CudaKernel` flattens `__cuda_args__()`
+on the CUDA path, so the call site is the same on both backends and each form
+can be built lazily on first access (a CPU run never builds device arguments):
+
+```python
+class ParticleArguments(xp.KernelArguments):
+    def __init__(self, markers):
+        self.markers = markers
+        self._host = None
+
+    def __host_args__(self):
+        if self._host is None:
+            self._host = MarkerArguments(self.markers)  # Pyccel class
+        return self._host
+
+    def __cuda_args__(self):
+        return (self.markers, self.markers.shape[0])
+
+
+kernel(particles.kernel_args, dt, n_threads=n)  # 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 `Array3D<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:
+
+```python
+class MarkerArguments:
+    def __init__(self, markers: "float[:, :]", n_markers: int, valid: "bool[:]"):
+        ...
+
+MarkerArgs = xp.CudaStruct.from_signature(MarkerArguments.__init__, "MarkerArgs")
+MarkerArgs.to_header("marker_args.cuh")  # Array2D<double> markers; long long n_markers; ...
+push = xp.CudaKernel(
+    r"""
+    #include "marker_args.cuh"
+    #include <cunumpy/index.cuh>
+    extern "C" __global__ void push(MarkerArgs m, double dt) {
+        CUNUMPY_THREAD_1D(ip, m.n_markers);
+        if (m.valid(ip)) m.markers(ip, 0) += dt * m.markers(ip, 3);
+    }""",
+    "push",
+    structs=[MarkerArgs],
+    include_dirs=["."],
+)
+push(MarkerArgs(markers=markers, n_markers=markers.shape[0], valid=valid), 0.1,
+     n_threads=markers.shape[0])
+```
+
+Launches can be 1D to 3D (`n_threads=(nx, ny)`, `block_size=(16, 16)`) or use
+an explicit `grid`, with dynamic shared memory (`shared_mem`) and a `stream`.
+C++ function templates are instantiated with `template_args`, and
+`CudaKernelVariants` caches kernels whose source is generated per variant
+(e.g. per dimension and dtype). See the [API reference](docs/source/api.md) for
+details.
+
+Kernels run asynchronously, so a CUDA error (an illegal memory access, say)
+normally surfaces at a later `.get()` or MPI call, far from the kernel that
+caused it. In debug mode, enabled with `xp.set_cuda_debug(True)`, the
+context manager `xp.cuda_debug()`, `CudaKernel(..., debug=True)` or the
+environment variable `CUNUMPY_CUDA_DEBUG=1`, kernels are compiled with
+`-lineinfo` and `-DCUNUMPY_BOUNDS_CHECK` and every launch is synchronized, so
+the error is raised as a `RuntimeError` naming the kernel and its launch shape.
+To find the faulting line and out-of-bounds accesses that do not crash, the
+next step is NVIDIA's memory checker:
+`CUNUMPY_CUDA_DEBUG=1 compute-sanitizer python -m pytest ...`.
+
+## Test kernel pairs
+
+`cunumpy.testing` helps to test the ports with pytest. `assert_kernels_agree`
+builds the arguments on both backends, runs the host and the CUDA kernel and
+compares the arrays they wrote; with `catalog.parity_cases()`, one
+parametrised test covers every ported kernel of a catalog. `BACKENDS` and
+`requires_cupy` parametrize tests over the backends, skipping CuPy without a
+GPU, and `device_function_kernel` wraps a `__device__` helper in an elementwise
+kernel so it can be checked against its host version without writing a test
+kernel:
+
+```python
+import pytest
+from cunumpy.testing import assert_kernels_agree
+
+
+def make_args(backend, seed):
+    x = xp.to_cunumpy(np.random.default_rng(seed).random(1000))
+    return (x, 2.0, x.size)
+
+
+@pytest.mark.parametrize("name, kernel", catalog.parity_cases())
+def test_parity(name, kernel):
+    assert_kernels_agree(kernel, make_args, n_threads=1000)
+
+Accumulation kernels often write into a buffer that another library owns on
+the host (a stencil vector's `_data`, exchanged over MPI). `DeviceMirror`
+pairs that NumPy array with a device copy: `mirror.device` is the CuPy array
+on the GPU and the host array itself on the CPU, `to_host()` copies back in
+place (the host array keeps its identity) and `zero()` clears the buffer, so
+the one transfer per accumulation is explicit. The shipped header
+`cunumpy/atomic.cuh` (found automatically, see `cuda_include_dir()`) provides
+`cunumpy_atomic_add()` and 2D/3D indexed variants for the many-threads-to-one-cell
+writes:
+
+```python
+mirror = xp.DeviceMirror(vector._data)
+mirror.zero()
+accumulate(markers, mirror.device, n_threads=n_markers)
+mirror.to_host()  # vector._data holds the result on both backends
+```
+
 ## Pyodide
 
 CuNumpy supports the NumPy backend in Pyodide. It does not provide CuPy/CUDA