From cbc7d3a947a27e6c06a9c566aa4cf3056e5aedf3 Mon Sep 17 00:00:00 2001 From: Max Date: Thu, 1 Oct 2026 15:56:45 +0200 Subject: [PATCH] Extend docs and add llm guide --- CHANGELOG.md | 2 + README.md | 41 ++- docs/source/ai-assistants.md | 4 + docs/source/api.md | 8 +- docs/source/best-practices.md | 70 +++++ docs/source/examples/index.md | 20 ++ docs/source/examples/particle-pusher.md | 324 +++++++++++++++++++++ docs/source/examples/portable-script.md | 94 ++++++ docs/source/guides/backends.md | 122 ++++++++ docs/source/guides/data-movement.md | 153 ++++++++++ docs/source/guides/gpu-devices.md | 97 ++++++ docs/source/guides/mpi.md | 123 ++++++++ docs/source/guides/portable-code.md | 171 +++++++++++ docs/source/guides/profiling.md | 94 ++++++ docs/source/index.md | 80 ++++- docs/source/installation.md | 82 ++++++ docs/source/kernels/accumulation.md | 115 ++++++++ docs/source/kernels/arguments.md | 211 ++++++++++++++ docs/source/kernels/cuda-kernel.md | 266 +++++++++++++++++ docs/source/kernels/debugging.md | 95 ++++++ docs/source/kernels/dispatch.md | 160 ++++++++++ docs/source/kernels/overview.md | 78 +++++ docs/source/kernels/pyccel-kernel.md | 117 ++++++++ docs/source/kernels/testing.md | 204 +++++++++++++ docs/source/quickstart.md | 309 +++++--------------- docs/source/troubleshooting.md | 94 ++++++ pyproject.toml | 2 +- src/cunumpy/LLM_GUIDE.md | 372 ++++++++++++++++++++++++ 28 files changed, 3252 insertions(+), 256 deletions(-) create mode 100644 docs/source/ai-assistants.md create mode 100644 docs/source/best-practices.md create mode 100644 docs/source/examples/index.md create mode 100644 docs/source/examples/particle-pusher.md create mode 100644 docs/source/examples/portable-script.md create mode 100644 docs/source/guides/backends.md create mode 100644 docs/source/guides/data-movement.md create mode 100644 docs/source/guides/gpu-devices.md create mode 100644 docs/source/guides/mpi.md create mode 100644 docs/source/guides/portable-code.md create mode 100644 docs/source/guides/profiling.md create mode 100644 docs/source/installation.md create mode 100644 docs/source/kernels/accumulation.md create mode 100644 docs/source/kernels/arguments.md create mode 100644 docs/source/kernels/cuda-kernel.md create mode 100644 docs/source/kernels/debugging.md create mode 100644 docs/source/kernels/dispatch.md create mode 100644 docs/source/kernels/overview.md create mode 100644 docs/source/kernels/pyccel-kernel.md create mode 100644 docs/source/kernels/testing.md create mode 100644 docs/source/troubleshooting.md create mode 100644 src/cunumpy/LLM_GUIDE.md diff --git a/CHANGELOG.md b/CHANGELOG.md index d5e5da3..cfd0d1f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 - `KernelCatalog.from_package(..., include_dirs=None)` is now an explicit keyword; by default the source root of the top-level package (the directory containing it) is an include directory of every CUDA kernel, in addition to the kernel's own folder, so kernels can `#include "my_pkg/common.cuh"`. ### Added +- Documentation restructured into getting started, user guide (backends, backend-agnostic code, data movement, devices, MPI, profiling), kernel porting guides (`PyccelKernel`, `CudaKernel`, `Kernel`/`KernelCatalog`, argument objects, accumulation, debugging, testing), worked examples, best practices and troubleshooting pages. +- `cunumpy/LLM_GUIDE.md`: a self-contained guide to the API and its rules for AI coding assistants, shipped as package data and rendered in the documentation. - `xp.CudaKernel`: Wraps a CUDA C kernel (`cupy.RawKernel`, compiled lazily with NVRTC) so it can be called with the same arguments as the host kernel it mirrors, plus `n_threads`. The `extern "C" __global__` signature is parsed once and every call is checked against it: argument count, array dtypes (host arrays raise, they are never copied), and scalars (Python scalars are cast to the declared C types with range checks; lossy or mismatching scalars raise instead of reaching the kernel as silently wrong values). Supports `block_size`, NVRTC `options`, `include_dirs`, `shared_mem`, `stream`, `CudaKernel.from_file()` (`_cuda.cu`), `compile()` and `prepare_args()`; `check_signature=False` skips the checks. - `xp.CudaArguments`: Base class for argument objects that are flattened into several CUDA kernel arguments; any object with a `__cuda_args__()` method is flattened. - `xp.parse_cuda_signature(source, name)` and `xp.CudaParameter`: Parse the parameters of a `__global__` function. diff --git a/README.md b/README.md index 090abdb..90a9000 100644 --- a/README.md +++ b/README.md @@ -372,6 +372,7 @@ class ParticleArguments(xp.KernelArguments): 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 @@ -446,6 +447,7 @@ def make_args(backend, seed): @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` @@ -473,6 +475,39 @@ installation example and compatibility notes. ## Documentation -The [user guide](docs/source/quickstart.md) explains common workflows. The -[API reference](docs/source/api.md) documents each helper and its behavior. -The [Pyodide guide](docs/source/pyodide.md) covers WebAssembly usage. +The full documentation lives in [`docs/source`](docs/source/index.md) and is +published at : + +* Getting started: [installation](docs/source/installation.md) and a + [quickstart](docs/source/quickstart.md) with a map of which guide covers what. +* User guide: [choosing a backend](docs/source/guides/backends.md), + [backend-agnostic code](docs/source/guides/portable-code.md), + [data movement](docs/source/guides/data-movement.md), + [devices, memory and streams](docs/source/guides/gpu-devices.md), + [MPI with one rank per GPU](docs/source/guides/mpi.md), + [timing and profiling](docs/source/guides/profiling.md). +* Porting kernels: [overview](docs/source/kernels/overview.md), + [`PyccelKernel`](docs/source/kernels/pyccel-kernel.md), + [`CudaKernel`](docs/source/kernels/cuda-kernel.md), + [`Kernel` and `KernelCatalog`](docs/source/kernels/dispatch.md), + [argument objects and structs](docs/source/kernels/arguments.md), + [accumulation kernels](docs/source/kernels/accumulation.md), + [debugging](docs/source/kernels/debugging.md), + [testing](docs/source/kernels/testing.md). +* [Worked examples](docs/source/examples/index.md), + [best practices](docs/source/best-practices.md), + [troubleshooting](docs/source/troubleshooting.md), + [Pyodide](docs/source/pyodide.md) and the + [API reference](docs/source/api.md). + +### For AI coding assistants + +[`src/cunumpy/LLM_GUIDE.md`](src/cunumpy/LLM_GUIDE.md) is a compact, +self-contained guide to the API and its rules for LLM-based coding assistants. +It ships inside the installed package, so an assistant working in a project that +depends on CuNumpy can read it from `site-packages/cunumpy/LLM_GUIDE.md`, or +locate it with: + +```bash +python -c "import cunumpy, pathlib; print(pathlib.Path(cunumpy.__file__).parent / 'LLM_GUIDE.md')" +``` diff --git a/docs/source/ai-assistants.md b/docs/source/ai-assistants.md new file mode 100644 index 0000000..c97db33 --- /dev/null +++ b/docs/source/ai-assistants.md @@ -0,0 +1,4 @@ + + +```{include} ../../src/cunumpy/LLM_GUIDE.md +``` diff --git a/docs/source/api.md b/docs/source/api.md index 645a1de..2305a91 100644 --- a/docs/source/api.md +++ b/docs/source/api.md @@ -241,6 +241,8 @@ too. An exception raised inside the block propagates as it is: def test_time_step_stays_on_the_device(): with xp.assert_no_transfers(): propagator(dt) +``` + ### `as_device_array(value, dtype=None, ndim=None, *, name=None)` The "reference or copy once" rule for building CUDA argument objects @@ -432,7 +434,7 @@ does not necessarily indicate a leak. ### `cuda_include_dir()` Returns the directory (as `str`) of the CUDA headers shipped with CuNumpy, -currently `cunumpy/atomic.cuh`. `CudaKernel` adds it to its NVRTC options as +`cunumpy/array_view.cuh`, `cunumpy/atomic.cuh` and `cunumpy/index.cuh`. `CudaKernel` adds it to its NVRTC options as `-I` automatically (and only once), so kernel sources can write `#include ` without configuration. Use it to pass the same headers to other compilers. @@ -668,11 +670,9 @@ kernels["shift"](x, 1.0, x.size, n_threads=x.size) fix it for this kernel, see "Debugging" below. Properties: `name`, `expression` (`name`, or the template instantiation such -as `"scale"`), `source`, `block_size`, `options`, `structs`, -`template_args`, `signature`, `is_compiled`, `debug`. as `"scale"`), `source`, `block_size`, `options`, `include_dirs`, `source_dir`, `included_headers`, `structs`, `template_args`, `signature`, -`is_compiled`. +`is_compiled`, `debug`. ### Included headers and the compile cache diff --git a/docs/source/best-practices.md b/docs/source/best-practices.md new file mode 100644 index 0000000..4ba86a5 --- /dev/null +++ b/docs/source/best-practices.md @@ -0,0 +1,70 @@ +# Best practices + +A condensed checklist. Each item links to the guide with the reasoning. + +## Structure + +* Import as `import cunumpy as xp` and always call `xp.(...)`; never + `from cunumpy import `. ([Backend-agnostic + code](guides/portable-code.md)) +* Choose the backend once, in the entry point, with `ARRAY_BACKEND` or + `set_backend()`. Library code never calls `set_backend()`. ([Choosing a + backend](guides/backends.md)) +* Read back `xp.get_backend()` after requesting CuPy; log it with + `xp.device_count()` and `xp.__version__`. +* Use `use_backend()` for scoped switches (tests, CPU reference computations), + never from several threads at once. +* Functions that receive arrays follow them with `get_array_module()`; functions + that combine arrays check them with `assert_same_backend()`. + +## Data + +* Create arrays with `xp.*` so they are born on the right device. Use `numpy` + directly only for host-only data and dtypes. +* Convert explicitly, at boundaries: `to_cunumpy()` on input, `to_numpy()` for + output, plotting and host-only libraries. ([Data + movement](guides/data-movement.md)) +* Keep the time loop free of transfers and of host synchronization (`float()`, + `.item()`, `print`, `if` on device values). Verify with + `assert_no_transfers()` in a test. +* Give dtypes explicitly for arrays that go to kernels. +* Generate random test data on the host with a seed, then convert; NumPy and + CuPy generators differ. + +## Kernels + +* Start with `KernelCatalog.from_package(..., missing_cuda="fallback")`, port + kernels by profile order, switch to `"raise"` when done. ([Porting + kernels](kernels/overview.md)) +* Keep the CUDA kernel's argument list identical to the host kernel's; put both + in one folder. +* Declare `outputs` on host kernels so the fallback copies back only what was + written, and never forget an argument that is written. +* Pass device arrays to `CudaKernel`; build argument objects once with + `as_device_array()`. ([Kernel arguments](kernels/arguments.md)) +* Use `long long` indices, `CUNUMPY_THREAD_1D` guards, `const` on inputs, and + `cunumpy_atomic_add` for scatter writes. +* Use array views (`Array2D`) instead of hand-computed offsets for + multi-dimensional data. +* Generate structs from the host argument class (`CudaStruct.from_signature`) + and test that committed headers are up to date. +* Compile at setup (`catalog.compile_all(jobs=...)`). +* Keep signature checks on; disable them (`check_signature=False`) only for + tiny kernels in hot loops after they are tested. + +## Verification + +* Parametrize tests with `BACKENDS` or the `backend` fixture; they run + everywhere and use the GPU where there is one. ([Testing + kernels](kernels/testing.md)) +* One `assert_kernels_agree` test over `catalog.parity_cases()`. +* Debug crashes with `CUNUMPY_CUDA_DEBUG=1`, then `compute-sanitizer`. + ([Debugging](kernels/debugging.md)) +* Time with `timed_region()`, profile with `nvtx_range()` and `nsys`. + ([Profiling](guides/profiling.md)) + +## MPI + +* `set_backend("cupy")`, `bind_local_device()`, then `from mpi4py import MPI`, + then `require_cuda_aware_mpi()`. ([MPI](guides/mpi.md)) +* `synchronize_for_mpi(*buffers)` before every MPI call with device buffers. diff --git a/docs/source/examples/index.md b/docs/source/examples/index.md new file mode 100644 index 0000000..aace660 --- /dev/null +++ b/docs/source/examples/index.md @@ -0,0 +1,20 @@ +# Worked examples + +Complete programs that combine the pieces described in the guides. + +* [A portable diffusion solver](portable-script.md): array-level code only, + one script for CPU and GPU, with a command-line backend switch, timing and + output at the boundaries. +* [Porting a particle-in-cell code](particle-pusher.md): a small simulation + whose kernels are ported to CUDA one at a time with a `KernelCatalog`, + verified with parity tests, and checked for transfers. + +For an MPI program with one rank per GPU, see the start-up sequence and halo +exchange in [Multi-GPU programs with MPI](../guides/mpi.md). + +```{toctree} +:hidden: + +portable-script +particle-pusher +``` diff --git a/docs/source/examples/particle-pusher.md b/docs/source/examples/particle-pusher.md new file mode 100644 index 0000000..b1f3f9d --- /dev/null +++ b/docs/source/examples/particle-pusher.md @@ -0,0 +1,324 @@ +# Porting a particle-in-cell code + +This example follows a small 1D electrostatic particle-in-cell (PIC) code from +a CPU-only program to a GPU port, kernel by kernel, the way a larger code base +would be ported. Every step leaves a working program. + +A PIC time step has four phases: + +1. **deposit**: each particle adds its charge to the grid cell it is in + (scatter, needs atomics on the GPU); +2. **field solve**: compute the electric field from the charge density (here + with an FFT, array-level code); +3. **accelerate**: each particle reads the field in its cell and updates its + velocity (gather); +4. **push**: each particle moves. + +Phases 1, 3 and 4 are loops over particles: kernels. Phase 2 is array code +that runs on either backend as it is. + +## The layout + +```text +pic/ +├── __init__.py +├── simulation.py +└── kernels/ + ├── __init__.py # the KernelCatalog + ├── deposit/deposit_kernels.py + ├── accelerate/accelerate_kernels.py + └── push/push_kernels.py +tests/ +└── test_kernels.py +``` + +Each kernel has its own folder, so its CUDA version can be added next to it +later (`push/push_cuda.cu`). The `kernels` folder and every kernel folder are +packages (an `__init__.py` in each, possibly empty). + +## Step 1: host kernels + +The host kernels are written in the Pyccel style: annotated loops that Pyccel +can compile to fast native code. They also run unchanged as plain Python, +which is what this example does. + +```python +# pic/kernels/push/push_kernels.py +from math import fmod + + +def push(x: "float[:]", v: "float[:]", n: int, dt: float, length: float): + for i in range(n): + x[i] = fmod(x[i] + dt * v[i], length) + if x[i] < 0.0: + x[i] += length +``` + +```python +# pic/kernels/deposit/deposit_kernels.py +def deposit(x: "float[:]", rho: "float[:]", n: int, dx: float, n_cells: int, weight: float): + for i in range(n): + cell = min(int(x[i] / dx), n_cells - 1) + rho[cell] += weight / dx +``` + +```python +# pic/kernels/accelerate/accelerate_kernels.py +def accelerate(x: "float[:]", v: "float[:]", e: "float[:]", n: int, dx: float, + n_cells: int, qm_dt: float): + for i in range(n): + cell = min(int(x[i] / dx), n_cells - 1) + v[i] += qm_dt * e[cell] +``` + +The catalog collects them. `outputs` tells the host wrapper which argument each +kernel writes, so the GPU fallback copies back only that one: + +```python +# pic/kernels/__init__.py +import cunumpy as xp + +OUTPUTS = {"push": (0,), "deposit": (1,), "accelerate": (1,)} + +catalog = xp.KernelCatalog.from_package( + __name__, + missing_cuda="fallback", + host_options=lambda name: {"outputs": OUTPUTS[name]}, +) +``` + +The simulation calls the kernels through the catalog and does the field solve +with array operations: + +```python +# pic/simulation.py +import numpy as np + +import cunumpy as xp + +from .kernels import catalog + + +class Simulation: + def __init__(self, n_particles=100_000, n_cells=64, length=2 * np.pi, seed=0): + rng = np.random.default_rng(seed) # host data: identical on both backends + x = rng.uniform(0.0, length, n_particles) + v = rng.normal(0.0, 1.0, n_particles) + 0.1 * np.sin(x) + + self.n, self.n_cells, self.length = n_particles, n_cells, length + self.dx = length / n_cells + self.weight = length / n_particles # mean density 1 + self.x = xp.to_cunumpy(x) + self.v = xp.to_cunumpy(v) + self.rho = xp.zeros(n_cells) + self.e = xp.zeros(n_cells) + self.k = xp.to_cunumpy(2 * np.pi * np.fft.rfftfreq(n_cells, d=self.dx)) + + def solve_field(self): + rho_k = xp.fft.rfft(self.rho - xp.mean(self.rho)) + e_k = xp.zeros_like(rho_k) + e_k[1:] = -1j * rho_k[1:] / self.k[1:] # E = -d(phi)/dx, phi'' = -rho + self.e[:] = xp.fft.irfft(e_k, n=self.n_cells) + + def step(self, dt): + self.rho[:] = 0.0 + catalog["deposit"](self.x, self.rho, self.n, self.dx, self.n_cells, + self.weight, n_threads=self.n) + self.solve_field() + catalog["accelerate"](self.x, self.v, self.e, self.n, self.dx, + self.n_cells, -dt, n_threads=self.n) + catalog["push"](self.x, self.v, self.n, dt, self.length, n_threads=self.n) + + def field_energy(self): + return 0.5 * float(xp.sum(self.e**2)) * self.dx +``` + +This runs on the CPU, and, thanks to `missing_cuda="fallback"`, already on the +GPU backend too: + +```python +import cunumpy as xp +from pic.kernels import catalog +from pic.simulation import Simulation + +xp.set_backend("cupy") +print(catalog.summary()) # CUDA kernels: 0 of 3 (missing: accelerate, deposit, push) + +sim = Simulation() +with xp.count_transfers() as counter: + sim.step(0.1) +print(counter.report()) +``` + +```text +6 transfer(s) through cunumpy (0 to_host, 0 to_device, 3 kernel_conversion, 3 fallback) + ... +``` + +Every kernel call copies its arrays to the host and back. The results are +correct, which is the baseline for the port, but slow. + +## Step 2: port the first kernel + +Add `push/push_cuda.cu` next to the host kernel. It mirrors the host kernel's +argument list one to one: + +```c +// pic/kernels/push/push_cuda.cu +#include + +extern "C" __global__ +void push(double* x, const double* v, long long n, double dt, double length) { + CUNUMPY_THREAD_1D(i, n); + double xi = fmod(x[i] + dt * v[i], length); + x[i] = xi < 0.0 ? xi + length : xi; +} +``` + +Nothing else changes. The catalog finds the file on the next start, and +`catalog["push"]` launches the CUDA kernel on the GPU backend: + +```text +CUDA kernels: 1 of 3 (missing: accelerate, deposit) +``` + +The Python `int` and `float` arguments are cast to `long long` and `double` as +declared; passing, say, a NumPy array for `x` would raise a `TypeError` instead +of being copied. + +## Step 3: test the port + +Before porting more, make sure the CUDA kernel computes what the host kernel +computes. One parametrised test covers every ported kernel: + +```python +# tests/test_kernels.py +import numpy as np +import pytest + +import cunumpy as xp +from cunumpy.testing import assert_kernels_agree +from pic.kernels import catalog + +N, N_CELLS, LENGTH = 10_000, 64, 2 * np.pi +DX = LENGTH / N_CELLS + + +def particles(seed): + rng = np.random.default_rng(seed) + x = xp.to_cunumpy(rng.uniform(0.0, LENGTH, N)) + v = xp.to_cunumpy(rng.normal(0.0, 1.0, N)) + return x, v + + +def push_args(backend, seed): + x, v = particles(seed) + return (x, v, N, 0.1, LENGTH) + + +def deposit_args(backend, seed): + x, _ = particles(seed) + return (x, xp.zeros(N_CELLS), N, DX, N_CELLS, LENGTH / N) + + +def accelerate_args(backend, seed): + x, v = particles(seed) + e = xp.to_cunumpy(np.sin(np.linspace(0.0, LENGTH, N_CELLS))) + return (x, v, e, N, DX, N_CELLS, -0.1) + + +MAKE_ARGS = {"push": push_args, "deposit": deposit_args, "accelerate": accelerate_args} + + +@pytest.mark.parametrize("name, kernel", catalog.parity_cases()) +def test_parity(name, kernel): + # deposit sums with atomics in arbitrary order: allow round-off + assert_kernels_agree(kernel, MAKE_ARGS[name], n_threads=N, rtol=1e-10) +``` + +Without a GPU, the test is skipped; on a GPU runner it compares every ported +kernel, and each newly added `.cu` file is tested automatically. + +## Step 4: port the remaining kernels + +The deposit writes many particles into few cells, so it needs atomic adds: + +```c +// pic/kernels/deposit/deposit_cuda.cu +#include +#include + +extern "C" __global__ +void deposit(const double* x, double* rho, long long n, double dx, + long long n_cells, double weight) { + CUNUMPY_THREAD_1D(i, n); + long long cell = min((long long)(x[i] / dx), n_cells - 1); + cunumpy_atomic_add(&rho[cell], weight / dx); +} +``` + +```c +// pic/kernels/accelerate/accelerate_cuda.cu +#include + +extern "C" __global__ +void accelerate(const double* x, double* v, const double* e, long long n, + double dx, long long n_cells, double qm_dt) { + CUNUMPY_THREAD_1D(i, n); + long long cell = min((long long)(x[i] / dx), n_cells - 1); + v[i] += qm_dt * e[cell]; +} +``` + +Now the parity test covers all three, and the step should not transfer +anything. Turn that into a test: + +```python +from cunumpy.testing import requires_cupy +from pic.simulation import Simulation + + +@requires_cupy +def test_step_stays_on_device(): + with xp.use_backend("cupy"): + sim = Simulation(n_particles=N) + sim.step(0.1) # warm-up: compiles the kernels + with xp.assert_no_transfers(): + sim.step(0.1) +``` + +## Step 5: lock it in and measure + +With every kernel ported, switch the catalog to `missing_cuda="raise"`, so a +kernel added later without a CUDA version fails immediately on the GPU instead +of silently copying. Compile everything at start-up and time the loop: + +```python +import cunumpy as xp +from pic.kernels import catalog +from pic.simulation import Simulation + +xp.set_backend("cupy") +if xp.cupy_backend: + catalog.compile_all(jobs=4) + +sim = Simulation(n_particles=1_000_000) +with xp.timed_region("100 steps") as timing: + for _ in range(100): + with xp.nvtx_range("step"): + sim.step(0.05) +print(f"{timing.elapsed / 100 * 1e3:.2f} ms/step, field energy {sim.field_energy():.4e}") +``` + +The same script, without `set_backend("cupy")`, runs the Pyccel-compiled (or, +as here, pure Python) host kernels on the CPU. + +## Where to go from here + +* The kernels take five to seven loose arguments. In a real code, group the + particle data with [`KernelArguments` or a `CudaStruct`](../kernels/arguments.md). +* If the charge density belongs to a host library (a distributed vector + exchanged over MPI), deposit into a [`DeviceMirror`](../kernels/accumulation.md). +* For several GPUs, split the particles over MPI ranks, bind one GPU per rank + and reduce `rho` across ranks ([Multi-GPU programs with + MPI](../guides/mpi.md)). diff --git a/docs/source/examples/portable-script.md b/docs/source/examples/portable-script.md new file mode 100644 index 0000000..492947f --- /dev/null +++ b/docs/source/examples/portable-script.md @@ -0,0 +1,94 @@ +# A portable diffusion solver + +This script solves the 2D heat equation with an explicit finite-difference +scheme. It uses only array operations, so it needs no kernels: the same file +runs on NumPy and CuPy. It shows the habits that keep such a script fast on the +GPU: + +* the backend is chosen once, at start-up, and reported; +* initial data is created on the host with a fixed seed, then moved to the + device once; +* the time loop has no transfers and no host synchronization; +* diagnostics and output convert explicitly, every `--every` steps only; +* timing waits for the device. + +```python +"""heat.py: explicit 2D diffusion, on CPU or GPU.""" + +import argparse + +import numpy as np + +import cunumpy as xp + + +def laplacian(u, dx): + """5-point Laplacian with periodic boundaries, on u's own backend.""" + array_xp = xp.get_array_module(u) + return ( + array_xp.roll(u, 1, axis=0) + + array_xp.roll(u, -1, axis=0) + + array_xp.roll(u, 1, axis=1) + + array_xp.roll(u, -1, axis=1) + - 4.0 * u + ) / dx**2 + + +def initial_condition(n, seed): + """Built on the host, so CPU and GPU runs start from identical data.""" + rng = np.random.default_rng(seed) + x = np.linspace(0.0, 1.0, n, endpoint=False) + gauss = np.exp(-((x[:, None] - 0.5) ** 2 + (x[None, :] - 0.5) ** 2) / 0.01) + return gauss + 0.01 * rng.standard_normal((n, n)) + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument("--gpu", action="store_true") + parser.add_argument("--n", type=int, default=512) + parser.add_argument("--steps", type=int, default=2000) + parser.add_argument("--every", type=int, default=500) + args = parser.parse_args() + + xp.set_backend("cupy" if args.gpu else "numpy") + print(f"backend={xp.get_backend()} devices={xp.device_count()} cunumpy={xp.__version__}") + + dx = 1.0 / args.n + dt = 0.2 * dx**2 # stable for the explicit scheme + u = xp.to_cunumpy(initial_condition(args.n, seed=0)) # one host-to-device copy + + with xp.timed_region("time loop") as timing: + for step in range(1, args.steps + 1): + u = u + dt * laplacian(u, dx) + if step % args.every == 0: + total = float(xp.sum(u)) * dx**2 # one scalar to the host + print(f"step {step:5d}: integral = {total:.6f}") + + print(f"{timing.elapsed:.3f} s ({timing.elapsed / args.steps * 1e6:.1f} us/step)") + np.save("heat_final.npy", xp.to_numpy(u)) # one device-to-host copy + + +if __name__ == "__main__": + main() +``` + +Run it: + +```bash +python heat.py # NumPy +python heat.py --gpu # CuPy, if available; falls back to NumPy otherwise +``` + +The integral printed every 500 steps is conserved by the periodic scheme, so it +is also a quick check that both backends compute the same thing. + +## Variations + +* **Check for transfers in a test.** Wrap a few steps in + `xp.assert_no_transfers()` to make sure nobody adds a `to_numpy()` to the loop + later. +* **Profile.** Mark the update with `xp.nvtx_range("update")` and run + `nsys profile -t cuda,nvtx python heat.py --gpu`. +* **Use it as a library.** `laplacian()` follows its input via + `get_array_module()`, so other code can call it with NumPy or CuPy arrays no + matter which backend is active. diff --git a/docs/source/guides/backends.md b/docs/source/guides/backends.md new file mode 100644 index 0000000..a39a432 --- /dev/null +++ b/docs/source/guides/backends.md @@ -0,0 +1,122 @@ +# Choosing a backend + +CuNumpy has two backends: `"numpy"` (CPU) and `"cupy"` (NVIDIA GPU). The +*active backend* decides which library `xp.` calls, and therefore +where newly created arrays live. + +## Select the backend at start-up + +The backend is NumPy unless the environment variable `ARRAY_BACKEND=cupy` is +set when CuNumpy is first imported: + +```bash +ARRAY_BACKEND=cupy python simulate.py +``` + +This is the least intrusive option for scripts and batch jobs: the code does +not change, and a job script decides whether it runs on a GPU. The variable is +read once, at import; setting it later in `os.environ` has no effect. + +To choose from inside the program, for example from a command-line flag, call +`set_backend()` once, early, before arrays are created: + +```python +import argparse + +import cunumpy as xp + +parser = argparse.ArgumentParser() +parser.add_argument("--gpu", action="store_true") +args = parser.parse_args() + +xp.set_backend("cupy" if args.gpu else "numpy") +print(f"running on {xp.get_backend()}") +``` + +## Check what you actually got + +Requesting CuPy is a request, not a guarantee. If CuPy is not installed, or +installed but not functional (no driver, no visible GPU, a CUDA version +mismatch), CuNumpy falls back to NumPy without raising. Always read the +effective backend back when it matters: + +```python +xp.set_backend("cupy") +if xp.get_backend() != "cupy": + raise SystemExit("this run needs a GPU, but CuPy is not usable") +``` + +`xp.cupy_available()` answers the question without changing the backend. +`xp.numpy_backend` and `xp.cupy_backend` are booleans for the active backend, +handy in conditionals: + +```python +if xp.cupy_backend: + xp.bind_local_device() +``` + +## Switch temporarily + +`use_backend()` selects a backend for a block and restores the previous one on +exit, also if the block raises. It is the right tool for tests, notebooks, and +CPU reference computations inside a GPU program: + +```python +xp.set_backend("cupy") + +with xp.use_backend("numpy"): + reference = xp.linspace(0.0, 1.0, 100) # NumPy array + +print(xp.get_backend()) # 'cupy' again +``` + +The selection is a single process-wide setting. Do not switch it from several +threads or async tasks at once: one task's `use_backend()` changes the backend +seen by all others. + +## The active backend is not the array's backend + +Changing the backend never moves existing arrays. A program can have NumPy +active while CuPy arrays are alive, and vice versa: + +```python +xp.set_backend("cupy") +on_gpu = xp.arange(4) + +xp.set_backend("numpy") +print(xp.get_backend()) # 'numpy': what xp.* creates now +print(xp.get_array_backend(on_gpu)) # 'cupy': where this array lives +``` + +Two families of functions answer the two questions: + +| Question | Functions | +| --- | --- | +| Which library does `xp.*` use now? | `get_backend()`, `numpy_backend`, `cupy_backend` | +| Where does this array live? | `get_array_backend(a)`, `is_cpu(a)`, `is_gpu(a)`, `get_array_module(a)` | +| Do these arrays live in the same place? | `same_backend(*arrays)`, `assert_same_backend(*arrays)` | + +Code that creates arrays uses the first family; code that receives arrays +uses the second (see [Writing backend-agnostic code](portable-code.md)). + +## Recommended patterns + +* **Choose once, at the top.** Select the backend in the entry point (the + `main()` of a script, a configuration loader), not inside library functions. + Library code should not call `set_backend()`. +* **Library code follows its inputs.** Use `get_array_module(array)` or check + with `assert_same_backend()` instead of reading the global setting. +* **Use `use_backend()` for scoped work.** It is exception-safe and makes the + scope obvious; a bare `set_backend()` in the middle of a function changes + the state for everything that runs afterwards. +* **Log the effective backend.** Print `xp.get_backend()`, `xp.device_count()` + and `xp.__version__` in the run's output so results can be traced to the + hardware they ran on. + +## What is behind `xp` + +`xp.` resolves `name` on `array_api_compat.numpy` or +`array_api_compat.cupy`, depending on the active backend. These are thin +compatibility layers over the real libraries, so the arrays are plain +`numpy.ndarray` and `cupy.ndarray` objects. [Why CuNumpy uses +`array-api-compat`](../array-api-compat.md) explains the layer in more detail. diff --git a/docs/source/guides/data-movement.md b/docs/source/guides/data-movement.md new file mode 100644 index 0000000..09a65b5 --- /dev/null +++ b/docs/source/guides/data-movement.md @@ -0,0 +1,153 @@ +# Moving data between host and device + +On a GPU system, NumPy arrays live in host (CPU) memory and CuPy arrays live in +device (GPU) memory. Copying between the two goes over PCIe or NVLink and is +often slower than the computation itself. CuNumpy therefore never moves data +implicitly: every transfer is a visible function call. + +## The three conversion functions + +| Function | Returns | Copies when | +| --- | --- | --- | +| `xp.to_numpy(a)` | `numpy.ndarray` (host) | `a` is a CuPy array (device to host) | +| `xp.to_cupy(a)` | `cupy.ndarray` (device) | `a` is not a CuPy array (host to device) | +| `xp.to_cunumpy(a)` | an array of the *active* backend | `a` lives on the other backend | + +* `to_numpy()` accepts anything `numpy.asarray` accepts (lists, tuples, NumPy + arrays, views). A NumPy array is returned as it is, without a copy. +* `to_cupy()` raises `ImportError` when CuPy or CUDA is not usable. A CuPy + array is returned as it is. +* `to_cunumpy()` is the portable choice: on the NumPy backend it behaves like + `to_numpy()`, on CuPy like `to_cupy()`. +* None of them modify the source array or change the active backend. + +## Transfer at boundaries, not in loops + +Load or generate data, move it to the device once, run the whole computation +there, and bring back only what the host needs: + +```python +import numpy as np + +import cunumpy as xp + +xp.set_backend("cupy") + +signal = xp.to_cunumpy(np.load("signal.npy")) # one host-to-device copy +spectrum = xp.abs(xp.fft.rfft(signal)) ** 2 +peak = int(xp.argmax(spectrum)) # one tiny device-to-host copy +np.save("spectrum.npy", xp.to_numpy(spectrum)) # one device-to-host copy +``` + +Typical boundaries are file I/O, plotting, calls into host-only libraries +(SciPy, h5py, matplotlib), and diagnostics that run every N steps rather than +every step. + +An anti-pattern to avoid: + +```python +for step in range(n_steps): + state = xp.to_cupy(state) # copied up every step + state = advance(state, dt) + state = xp.to_numpy(state) # and down again + if step % 100 == 0: + write_output(state) +``` + +Move the conversions out of the loop and convert inside the `if` only. + +## Count the transfers + +The classic performance bug of a GPU port is a transfer that sneaks into the +time loop. `count_transfers()` records every copy made through CuNumpy in a +block, with the file and line that caused it: + +```python +with xp.count_transfers() as counter: + for _ in range(10): + step(state, dt) + +print(counter.total) +print(counter.report()) +``` + +```text +20 transfer(s) through cunumpy (10 to_host, 10 to_device, 0 kernel_conversion, 0 fallback) + to_host (10): + /home/me/sim/diagnostics.py:42: to_numpy(shape=(100000,), dtype=float64) (x10) + to_device (10): + /home/me/sim/step.py:17: to_cupy(shape=(100000,), dtype=float64) (x10) +``` + +Four kinds of events are recorded: `to_host`, `to_device`, +`kernel_conversion` (a [`PyccelKernel`](../kernels/pyccel-kernel.md) that copied +device arrays to the host and back) and `fallback` (a +[`Kernel`](../kernels/dispatch.md) without CUDA version running its host kernel +on the GPU backend). Calls that do not copy, such as `to_numpy()` of a NumPy +array, are not counted, so a `count_transfers()` block on the NumPy backend +reports zero. + +In tests, `assert_no_transfers()` turns this into a check that fails with the +report: + +```python +def test_step_stays_on_device(): + state = make_state() + with xp.assert_no_transfers(): + step(state, dt) +``` + +The counter only sees transfers made through CuNumpy. Raw `cupy.asarray(host)`, +`device_array.get()`, `float(device_scalar)` and conversions inside other +libraries are invisible to it; use `nsys` to find those (see [Timing and +profiling](profiling.md)). + +## Build device arguments once: `as_device_array` + +Objects that hold arrays for CUDA kernels (see [Kernel arguments and +structs](../kernels/arguments.md)) should reference existing device arrays and +copy only what is not already in the right form. `as_device_array(value, dtype, +ndim=None, *, name=None)` implements that rule: + +* a C-contiguous CuPy array with the requested dtype is returned unchanged (the + same object, so kernels write into the caller's array); +* anything else (a tuple such as `(3, 3, 3)`, a host array, another dtype, a + non-contiguous view) becomes one C-contiguous device copy. + +```python +import numpy as np + +degree = xp.as_device_array((3, 3, 3), np.int32, ndim=1, name="degree") +markers = xp.as_device_array(markers, np.float64, ndim=2, name="markers") +``` + +Call it when the argument object is built, never per kernel call. On the NumPy +backend it raises `RuntimeError`, so host data is never copied to a device by +accident. + +## Buffers owned by another library: `DeviceMirror` + +Sometimes the array a GPU kernel must write into belongs to a host library (a +distributed vector exchanged over MPI, a buffer a solver keeps using). +`DeviceMirror(host_array)` pairs such an array with a device copy and makes the +two transfers explicit: `mirror.device` is the array kernels write into, and +`mirror.to_host()` copies the result back into the original host array in +place. On the NumPy backend `mirror.device` *is* the host array and the copies +are no-ops. See [Accumulation kernels](../kernels/accumulation.md) for the full +pattern. + +## Pinned memory + +Host-to-device copies from page-locked ("pinned") host memory are faster and +can overlap with computation on a stream. `xp.pin_memory(host_array)` returns a +pinned copy: + +```python +pinned = xp.pin_memory(np.load("snapshot.npy")) +with xp.stream(): + device = xp.to_cupy(pinned) +``` + +Pinned memory is a limited system resource; use it for large, repeatedly +transferred buffers after a profile shows transfers matter. +`pin_memory()` requires CuPy. diff --git a/docs/source/guides/gpu-devices.md b/docs/source/guides/gpu-devices.md new file mode 100644 index 0000000..4a5f2f6 --- /dev/null +++ b/docs/source/guides/gpu-devices.md @@ -0,0 +1,97 @@ +# Devices, memory and streams + +These helpers wrap the parts of CuPy's device API that a portable program +needs. All of them are safe to call on the NumPy backend, where they are no-ops +or return neutral values, so the calls can stay in the code unconditionally. + +| Function | CuPy backend | NumPy backend | +| --- | --- | --- | +| `cupy_available()` | `True` if CuPy imports and works | same check (backend-independent) | +| `device_count()` | number of visible GPUs | number of visible GPUs, `0` without CuPy | +| `set_device(i)` | makes device `i` current | no-op | +| `memory_info()` | `(free_bytes, total_bytes)` of the current device | `None` | +| `free_memory()` | releases cached blocks of CuPy's pools | no-op | +| `synchronize()` | waits for the current device | no-op | +| `stream()` | yields a new non-blocking stream | yields `None` | +| `pin_memory(a)` | pinned host copy | raises `ImportError` without CuPy | + +## Select a GPU + +On a multi-GPU workstation, pick the device before creating arrays: + +```python +xp.set_backend("cupy") +print("GPUs:", xp.device_count()) +xp.set_device(1) # arrays created from now on live on GPU 1 +values = xp.zeros(10**6) +``` + +The usual alternative is to restrict visibility from outside the process, +which also works for libraries that do not know about CuNumpy: + +```bash +CUDA_VISIBLE_DEVICES=1 ARRAY_BACKEND=cupy python simulate.py +``` + +For MPI programs with one rank per GPU, use `bind_local_device()` instead +(see [Multi-GPU programs with MPI](mpi.md)). + +## Memory + +```python +free, total = xp.memory_info() +print(f"{free / 2**30:.1f} of {total / 2**30:.1f} GiB free") +``` + +`memory_info()` reports the CUDA runtime's view of the whole device, including +memory held by other processes and by CuPy's caches. + +CuPy does not return freed memory to the driver. It keeps released blocks in a +memory pool and reuses them for the next allocation, which makes allocation +cheap. As a consequence, `nvidia-smi` shows memory as used after arrays went +out of scope. That is not a leak. `free_memory()` returns the currently +unused cached blocks to the driver, for example before handing the GPU to +another library or between phases with very different memory needs: + +```python +del large_temporary +xp.free_memory() +``` + +It cannot free memory still referenced by live arrays. If memory keeps +growing, look for arrays kept alive by lists, caches or closures. + +## Asynchronous execution and synchronization + +A CuPy operation returns as soon as the work is queued; the GPU runs it later. +This is what makes GPUs fast, and it has three practical consequences: + +1. Reading a value on the host (`to_numpy()`, `float()`, `print`) waits for all + queued work that produces it. This happens automatically. +2. Timing with `time.perf_counter()` around a launch measures the launch, not + the work. Use `xp.timed_region()` (see [Timing and profiling](profiling.md)). +3. Libraries that read device memory without CuPy's knowledge, most + importantly MPI, need an explicit `xp.synchronize()` (or + `xp.synchronize_for_mpi()`) before they access a buffer. + +## Streams + +Work on one stream runs in order; work on different streams may overlap. +`xp.stream()` creates a non-blocking stream and makes it current for the block: + +```python +with xp.stream() as s: + device = xp.to_cupy(pinned_host) # copy and compute queued on s + result = xp.fft.fft(device) + +xp.synchronize() # wait before using the result elsewhere +``` + +On NumPy, `stream()` yields `None`, so do not call methods on the yielded value +in portable code; use `xp.synchronize()` instead of `s.synchronize()`. +CUDA kernels take a `stream=` argument as well (see [Writing CUDA +kernels](../kernels/cuda-kernel.md)). + +Streams help when independent pieces of work (a transfer and an unrelated +kernel, or several small kernels) can overlap. Most array-level code gains +nothing from them; reach for streams after a profile shows idle gaps. diff --git a/docs/source/guides/mpi.md b/docs/source/guides/mpi.md new file mode 100644 index 0000000..10c440e --- /dev/null +++ b/docs/source/guides/mpi.md @@ -0,0 +1,123 @@ +# Multi-GPU programs with MPI + +The common layout for GPU simulations is one MPI rank per GPU. Getting it right +takes four steps, each of which fails silently or with a segfault if skipped. +CuNumpy provides one helper per step. + +## The start-up sequence + +```python +import cunumpy as xp + +xp.set_backend("cupy") +xp.bind_local_device() # 1. pick this rank's GPU, before MPI_Init + +from mpi4py import MPI # 2. MPI_Init happens here + +xp.require_cuda_aware_mpi() # 3. fail clearly if MPI cannot take GPU buffers + +comm = MPI.COMM_WORLD +send = xp.full(1000, comm.rank, dtype=xp.float64) +recv = xp.empty_like(send) + +xp.synchronize_for_mpi(send, recv) # 4. before every MPI call on device buffers +comm.Sendrecv(send, dest=(comm.rank + 1) % comm.size, + recvbuf=recv, source=(comm.rank - 1) % comm.size) +``` + +The same file runs on the NumPy backend: `bind_local_device()` returns `None`, +`require_cuda_aware_mpi()` does nothing, and `synchronize_for_mpi()` ignores +host buffers. + +### 1. Bind the device before `MPI_Init` + +Without device selection every rank on a node uses GPU 0. A CUDA-aware MPI also +binds to whatever device is current when `MPI_Init` runs, so the device must be +chosen *before* `from mpi4py import MPI`. At that point MPI cannot be asked for +the rank yet, but launchers export the node-local rank in environment +variables, which `xp.local_rank()` reads (Open MPI, MVAPICH2, Intel MPI/MPICH, +PMI, Cray PALS, Slurm and `LOCAL_RANK`). + +`xp.bind_local_device()` selects device `local_rank() % device_count()`, +creates its CUDA context, and returns the device id. If the launcher already +restricts each rank to one GPU (`CUDA_VISIBLE_DEVICES` per rank, or Slurm's +`--gpus-per-task=1`), every process sees one device and selects it. + +`set_device_for_rank(rank)` is an older alternative that needs the global MPI +rank and assumes ranks are numbered contiguously per node. Prefer +`bind_local_device()`. + +### 2. Initialize MPI + +Importing `mpi4py.MPI` initializes MPI by default. Do it after step 1. + +### 3. Check that MPI is CUDA-aware + +Only an MPI library built with CUDA support can send and receive CuPy arrays +directly. A plain build reads the device pointer as a host address: the result +is a segfault, or garbage values without any error. `require_cuda_aware_mpi()` +performs a tiny `Sendrecv` of a device array on every rank, verifies the +received values, and raises a `RuntimeError` with hints on getting a +CUDA-aware build if it fails. + +The check is collective: call it once, on all ranks, at start-up. The boolean +variant `mpi_is_cuda_aware(comm)` lets a program choose a fallback instead, +such as staging buffers through the host: + +```python +if xp.cupy_backend and not xp.mpi_is_cuda_aware(comm): + stage_through_host = True +``` + +A segfault *inside* the check means the same as `False`: the MPI library is not +CUDA-aware. + +### 4. Synchronize before every MPI call + +CuPy kernels run asynchronously, and MPI knows nothing about CUDA streams. A +buffer that a kernel is still writing would be sent as it is at that moment. +`synchronize_for_mpi(*buffers)` waits for the current stream if any argument is +a CuPy array, and costs nothing otherwise: + +```python +def exchange_halo(field, comm, left, right): + # field has shape (nx + 2, ny): one ghost row on each side + send_l, send_r = field[1], field[-2] + recv_l, recv_r = xp.empty_like(send_l), xp.empty_like(send_r) + xp.synchronize_for_mpi(send_l, send_r) + comm.Sendrecv(send_l, dest=left, recvbuf=recv_r, source=right) + comm.Sendrecv(send_r, dest=right, recvbuf=recv_l, source=left) + field[0], field[-1] = recv_l, recv_r +``` + +Note that MPI buffers must be contiguous: rows of a C-ordered array are, +columns are not (copy them with `xp.ascontiguousarray()` first). No +synchronization is needed after MPI returns; kernels launched afterwards see +the received data. + +## Launching + +```bash +# Open MPI, 4 GPUs on one node +mpirun -n 4 python simulate.py --gpu + +# Slurm, 2 nodes with 4 GPUs each +srun --nodes=2 --ntasks-per-node=4 --gpus-per-task=1 python simulate.py --gpu +``` + +Print the binding once at start-up to catch mapping errors early: + +```python +device = xp.bind_local_device() +from mpi4py import MPI +print(f"rank {MPI.COMM_WORLD.rank}: local rank {xp.local_rank()}, device {device}") +``` + +## Common failures + +| Symptom | Cause | +| --- | --- | +| all ranks on a node run on GPU 0 | no device binding, or binding after `MPI_Init` | +| segfault in the first `Send`/`Recv` of a CuPy array | MPI is not CUDA-aware; `require_cuda_aware_mpi()` turns this into a clear error | +| received data is occasionally stale or zero | missing `synchronize_for_mpi()` before the call | +| `BufferError` or wrong values with array slices | non-contiguous buffer passed to MPI | diff --git a/docs/source/guides/portable-code.md b/docs/source/guides/portable-code.md new file mode 100644 index 0000000..abc66ca --- /dev/null +++ b/docs/source/guides/portable-code.md @@ -0,0 +1,171 @@ +# Writing backend-agnostic code + +The goal is one code base that runs unchanged on NumPy and CuPy. Most array +code already does; this page collects the rules that keep it that way. + +## Always go through the `xp` namespace + +```python +import cunumpy as xp + +grid = xp.linspace(0.0, 1.0, 1025) +field = xp.exp(-((grid - 0.5) ** 2) / 0.01) +energy = xp.sum(field**2) * (grid[1] - grid[0]) +``` + +Two habits break portability: + +* **Importing raw NumPy for array work.** `np.zeros(n)` always creates a host + array, even when the rest of the program runs on the GPU. Use `numpy` + directly only for things that must stay on the host: dtypes such as + `np.float64`, file I/O, and data you will explicitly convert. +* **Binding names from `cunumpy`.** `xp.zeros` is looked up each time it is + accessed, so it follows the active backend. `from cunumpy import zeros` + captures the function of the backend that was active at import time and + keeps calling it after `set_backend()`. Always write `xp.zeros(...)`. The + same applies to `from cunumpy import numpy_backend`; read `xp.numpy_backend` + where you need it. + +## Follow the input, not the global setting + +A function that *creates* arrays uses the active backend. A function that +*receives* arrays should compute on the library of its inputs, which may differ +from the active backend (a caller may have kept a NumPy reference while CuPy is +active). `get_array_module()` returns the matching compatibility module: + +```python +def gradient(f, dx): + array_xp = xp.get_array_module(f) + out = array_xp.empty_like(f) + out[1:-1] = (f[2:] - f[:-2]) / (2 * dx) + out[0] = (f[1] - f[0]) / dx + out[-1] = (f[-1] - f[-2]) / dx + return out +``` + +When a function combines several inputs, check that they agree. A mixed +NumPy/CuPy operation otherwise fails somewhere deep inside the library, or, +worse, silently converts: + +```python +def axpy(a, x, y): + xp.assert_same_backend(x, y) # TypeError naming both backends + return a * x + y +``` + +At the boundary of a library (the public constructor of a solver, for example), +normalize incoming arrays to the active backend once with `to_cunumpy()` and +keep them there: + +```python +class Solver: + def __init__(self, matrix, rhs): + self.matrix = xp.to_cunumpy(matrix) # copied once if needed + self.rhs = xp.to_cunumpy(rhs) +``` + +## Be explicit about dtypes + +NumPy and CuPy agree on dtype rules for array operations, but Python scalars and +lists are inferred, and code that mixes `float32` arrays with `float64` arrays +promotes silently. Pass a dtype when it matters, especially for arrays that will +be handed to CUDA kernels, which check dtypes strictly: + +```python +x = xp.zeros(n, dtype=xp.float64) +idx = xp.arange(n, dtype=xp.int32) +w = xp.asarray([0.25, 0.75], dtype=xp.default_float_dtype()) +``` + +`np.float64` and `xp.float64` are the same dtype object and both are accepted by +CuPy, so using NumPy dtypes for declarations is fine. + +## Random numbers + +`xp.get_rng(seed)` returns a `Generator` of the active backend, so random data is +generated where it is used: + +```python +rng = xp.get_rng(seed=42) +velocities = rng.normal(0.0, 1.0, size=(n_particles, 3)) +``` + +NumPy and CuPy generators do **not** produce the same sequence from the same +seed. When CPU and GPU runs must start from identical data (regression tests, +comparing results across machines), generate on the host and convert: + +```python +import numpy as np + +host = np.random.default_rng(42).normal(size=(n_particles, 3)) +velocities = xp.to_cunumpy(host) +``` + +## Watch for hidden synchronization + +On the GPU, operations are queued and run asynchronously. Anything that needs a +value on the host waits for the GPU and copies data back: + +* `float(x)`, `int(x)`, `bool(x)`, `x.item()` on a device scalar; +* `if xp.any(mask): ...`, `while residual > tol: ...` (a device comparison in a + Python condition); +* `print(array)`, `array.tolist()`, f-string formatting of a device value; +* iterating over a device array in a Python `for` loop. + +None of these are errors, and they are free on NumPy. In a hot loop they stall +the GPU. Check convergence every few iterations instead of every iteration, and +keep reductions on the device until a value is really needed. These copies +bypass CuNumpy and are not seen by `count_transfers()`; use a profiler (see +[Timing and profiling](profiling.md)). + +## Know where the libraries differ + +CuNumpy does not wrap operations; it forwards them. Coverage is therefore that +of the installed CuPy version. The common differences: + +* Some NumPy functions do not exist in CuPy, or exist with fewer options (for + example parts of `numpy.polynomial`, object and string dtypes, structured + arrays). Check the CuPy documentation for functions outside the core. +* Reductions such as `xp.sum(a)` return a NumPy scalar on NumPy but a 0-d + device array on CuPy. Keep it as an array, or call `float()` and accept the + synchronization. +* SciPy functions accept only NumPy arrays; CuPy has its own `cupyx.scipy`. + Convert with `to_numpy()` at that boundary, or branch on + `xp.get_array_backend(a)`. +* Plotting, HDF5 and most I/O libraries need host arrays: `to_numpy()` first. + +When a function genuinely needs a different implementation per backend, branch +on the input's backend in one place rather than throughout the code: + +```python +def solve_banded(ab, b): + if xp.is_gpu(b): + import cupyx.scipy.linalg as la + else: + import scipy.linalg as la + return la.solve_banded((1, 1), ab, b) +``` + +## Test on both backends + +`cunumpy.testing.BACKENDS` parametrizes a test over NumPy and CuPy, skipping the +CuPy case where there is no GPU, so the same test file runs on a laptop and in +GPU CI: + +```python +import pytest +from cunumpy.testing import BACKENDS + +import cunumpy as xp + + +@pytest.mark.parametrize("backend", BACKENDS) +def test_gradient_of_linear_function(backend): + with xp.use_backend(backend): + x = xp.linspace(0.0, 1.0, 11) + g = gradient(3.0 * x, x[1] - x[0]) + assert xp.get_array_backend(g) == backend + assert float(xp.max(xp.abs(g - 3.0))) < 1e-12 +``` + +See [Testing kernels](../kernels/testing.md) for the other test helpers. diff --git a/docs/source/guides/profiling.md b/docs/source/guides/profiling.md new file mode 100644 index 0000000..ae594e9 --- /dev/null +++ b/docs/source/guides/profiling.md @@ -0,0 +1,94 @@ +# Timing and profiling + +GPU code needs two adjustments to the usual timing habits: work is +asynchronous, so a timer must wait for the device, and the interesting events +(kernels, copies) happen outside Python, so a GPU profiler is needed to see +them. CuNumpy's helpers handle both and run unchanged on the NumPy backend. + +## Time a region: `timed_region` + +```python +with xp.timed_region("field solve") as timing: + solve(field) + +print(f"{timing.name}: {timing.elapsed * 1e3:.2f} ms (synced={timing.synced})") +``` + +On CuPy, `timed_region` synchronizes the device on entry (so earlier queued +work is not charged to this region) and on exit (so the region's own work is +included), and marks the region as an NVTX range. On NumPy it is a plain +`time.perf_counter()` timer and `synced` is `False`. With `sync=False` only the +host time, that is the launch overhead, is measured. + +Compare with the naive version, which on a GPU typically reports a few +microseconds regardless of the problem size: + +```python +start = time.perf_counter() +solve(field) # only queues the kernels +elapsed = time.perf_counter() - start # wrong on the GPU +``` + +Timing tips: + +* **Warm up first.** The first call of a CUDA kernel compiles it (or loads it + from CuPy's disk cache), and the first CuPy operations initialize the CUDA + context. Run one step before measuring, or compile at setup with + `catalog.compile_all()` (see [Pairing host and CUDA kernels](../kernels/dispatch.md)). +* **Time many iterations.** Wrap a loop of steps rather than one short step; + the synchronization itself costs a few microseconds. +* **Do not leave synchronizing timers in production hot loops.** Each one + stalls the pipeline; use NVTX ranges there and a profiler to read them. + +## Mark regions for Nsight: `nvtx_range` + +```python +with xp.nvtx_range("push markers"): + push(markers, dt, n_threads=n_markers) + + +@xp.nvtx_range("time step") +def step(state, dt): + ... +``` + +On CuPy this pushes an NVTX range, so the region appears as a labelled bar on +the Nsight Systems timeline directly above the kernels it launched. It does not +synchronize and costs well under a microsecond. On NumPy, or without NVTX +support, it does nothing. `color` selects an entry of NVTX's colour table. + +Record a profile with: + +```bash +nsys profile -t cuda,nvtx -o step_profile python simulate.py --gpu --steps 20 +nsys stats step_profile.nsys-rep # summary tables in the terminal +``` + +and open the `.nsys-rep` file in Nsight Systems to see the timeline. + +## Find unwanted transfers + +A host-device copy in the time loop is the most common reason a GPU port is +slower than expected. `count_transfers()` lists every copy made through CuNumpy +with its call site: + +```python +with xp.count_transfers() as counter: + for _ in range(10): + step(state, dt) +print(counter.report()) +``` + +Copies made outside CuNumpy (`float(x)` on a device value, `cupy.asarray`, +another library converting internally) show up in Nsight as `memcpy` rows on the +timeline; NVTX ranges tell you which part of the step issued them. See [Moving +data between host and device](data-movement.md). + +## A profiling checklist + +1. Run a few steps with `count_transfers()` and remove every transfer from the + step. +2. Wrap the phases of a step in `nvtx_range()` and record an `nsys` profile. +3. Look for gaps between kernels (host overhead, synchronization) and for the + longest kernels. +4. Use `timed_region()` around a loop of steps to measure the improvement. diff --git a/docs/source/index.md b/docs/source/index.md index 29b1c3f..54d3585 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -1,28 +1,86 @@ # CuNumpy -CuNumpy provides a NumPy-like interface that can create and operate on NumPy -arrays on the CPU or CuPy arrays on an NVIDIA GPU. Select the backend once -for a program or temporarily for a section of code; use explicit conversion -helpers when data needs to cross between host and device. +CuNumpy lets one Python code base run on NumPy (CPU) or CuPy (NVIDIA GPU). +Replace `import numpy as np` with `import cunumpy as xp`, choose a backend, and +the array code you already have runs on the selected device. For codes whose +time is spent in compiled loops, CuNumpy adds a kernel layer to port those +kernels to CUDA one at a time, with the CPU version kept as the tested +reference. ```python import cunumpy as xp -xp.set_backend("cupy") +xp.set_backend("cupy") # falls back to NumPy without a usable GPU values = xp.arange(1_000) -print(xp.get_backend()) # 'cupy' when CuPy/CUDA is functional +print(xp.sum(values**2), xp.get_backend()) ``` -Install with `python -m pip install cunumpy`. GPU use also requires a CuPy -installation that matches the system's CUDA setup. CuNumpy falls back to -NumPy if CuPy cannot be used. +Install with `python -m pip install cunumpy`, plus a CuPy wheel matching your +CUDA version for GPU use. Start with the [Quickstart](quickstart.md). + +## Where to go + +* **New to CuNumpy**: [Installation](installation.md), then the + [Quickstart](quickstart.md). +* **Running array code on CPU and GPU**: the *User guide* pages on backends, + backend-agnostic code, data movement, devices, MPI, and profiling. +* **Porting compiled kernels to CUDA**: start at [Porting kernels to the + GPU](kernels/overview.md). +* **Complete programs**: [Worked examples](examples/index.md). +* **Looking up a function**: the [API reference](api.md). +* **Using an AI assistant**: point it at [the guide for AI + assistants](ai-assistants.md), which ships with the package as + `cunumpy/LLM_GUIDE.md`. ```{toctree} :maxdepth: 2 -:caption: Guides and reference: +:caption: Getting started +installation quickstart +``` + +```{toctree} +:maxdepth: 2 +:caption: User guide + +guides/backends +guides/portable-code +guides/data-movement +guides/gpu-devices +guides/mpi +guides/profiling array-api-compat -api +``` + +```{toctree} +:maxdepth: 2 +:caption: Porting kernels + +kernels/overview +kernels/pyccel-kernel +kernels/cuda-kernel +kernels/dispatch +kernels/arguments +kernels/accumulation +kernels/debugging +kernels/testing +``` + +```{toctree} +:maxdepth: 2 +:caption: Examples + +examples/index +``` + +```{toctree} +:maxdepth: 1 +:caption: More + +best-practices +troubleshooting pyodide +ai-assistants +api ``` diff --git a/docs/source/installation.md b/docs/source/installation.md new file mode 100644 index 0000000..3d2a3fe --- /dev/null +++ b/docs/source/installation.md @@ -0,0 +1,82 @@ +# Installation + +## Install the package + +```bash +python -m pip install cunumpy +``` + +This installs CuNumpy with its two dependencies, NumPy and `array-api-compat`. +That is all that is needed for the CPU (NumPy) backend, and CuNumpy runs on any +platform NumPy runs on, including [Pyodide](pyodide.md). CuNumpy requires +Python 3.10 or newer. + +## Add GPU support + +The GPU backend uses [CuPy](https://cupy.dev). CuNumpy does not install CuPy or +CUDA itself, because the right CuPy package depends on the CUDA version of the +machine. Install the CuPy wheel that matches your CUDA toolkit or driver, for +example: + +```bash +python -m pip install cupy-cuda12x # CUDA 12.x +``` + +On HPC clusters, load the site's CUDA module first and check the CuPy +installation guide for the matching package name. Then verify that CuNumpy can +use it: + +```python +import cunumpy as xp + +print("CuPy usable:", xp.cupy_available()) +print("visible GPUs:", xp.device_count()) + +xp.set_backend("cupy") +print("active backend:", xp.get_backend()) # 'cupy' if the GPU works +``` + +If `cupy_available()` is `False`, CuNumpy quietly falls back to NumPy when CuPy +is requested. See [Troubleshooting](troubleshooting.md) for the usual causes. + +## Optional extras + +| Extra | Installs | Use it for | +| --- | --- | --- | +| `cunumpy[test]` | `pytest`, `coverage` | running the test suite, using `cunumpy.testing` | +| `cunumpy[test-compiled]` | the above plus `pyccel` | tests that compile host kernels with Pyccel | +| `cunumpy[docs]` | Sphinx, MyST, the book theme | building this documentation | +| `cunumpy[dev]` | all of the above plus formatters | developing CuNumpy itself | + +For MPI programs, install `mpi4py` against your MPI library as usual. For GPU +MPI the library must be CUDA-aware; see [Multi-GPU programs with +MPI](guides/mpi.md). + +## Install from source + +```bash +git clone https://github.com/max-models/cunumpy.git +cd cunumpy +python -m pip install -e ".[dev]" +python -m pytest tests/unit +``` + +Tests that need a GPU are skipped automatically where CuPy is not functional. + +## Environment variables + +| Variable | Effect | +| --- | --- | +| `ARRAY_BACKEND=cupy` | start with the CuPy backend instead of NumPy (read once, at import) | +| `CUNUMPY_CUDA_DEBUG=1` | enable [CUDA debug mode](kernels/debugging.md) for all kernels | + +MPI launchers also export node-local rank variables (`OMPI_COMM_WORLD_LOCAL_RANK`, +`SLURM_LOCALID`, ...), which `xp.local_rank()` reads to pick a GPU per process. + +## Build the documentation + +```bash +python -m pip install ".[docs]" +cd docs +make html # output in docs/build/html +``` diff --git a/docs/source/kernels/accumulation.md b/docs/source/kernels/accumulation.md new file mode 100644 index 0000000..4ccad9d --- /dev/null +++ b/docs/source/kernels/accumulation.md @@ -0,0 +1,115 @@ +# Accumulation kernels + +Particle-in-cell and finite-element codes often *scatter* contributions: every +particle adds its charge to the grid cells it overlaps. Two problems appear on +the GPU: + +1. **Many threads write the same cell.** A plain `+=` loses updates when + threads race. The writes must be atomic. +2. **The target buffer belongs to someone else.** It is often a NumPy array + owned by a host library (a distributed stencil vector, a solver's + right-hand side) that keeps using it on the host, for example for MPI + reductions. It cannot simply be replaced by a CuPy array. + +CuNumpy addresses the first with the header `cunumpy/atomic.cuh` and the second +with `DeviceMirror`. + +## Atomic adds: `cunumpy/atomic.cuh` + +```c +#include + +double cunumpy_atomic_add(double* p, double v); // *p += v atomically, returns old value +float cunumpy_atomic_add(float* p, float v); + +// indexed variants for C-contiguous arrays of shape (n0, n1) and (n0, n1, n2) +double cunumpy_atomic_add_2d(double* data, long long n1, long long i, long long j, double v); +double cunumpy_atomic_add_3d(double* data, long long n1, long long n2, + long long i, long long j, long long k, double v); +``` + +The header is found by every `CudaKernel` automatically. On compute capability +6.0 and newer, `double` atomics are a hardware instruction; older devices use a +compare-and-swap loop. + +## `DeviceMirror` + +```python +mirror = xp.DeviceMirror(host_array) +``` + +| Member | CuPy backend | NumPy backend | +| --- | --- | --- | +| `mirror.device` | CuPy array, allocated on first access as a copy of the host | the host array itself | +| `mirror.to_device()` | copies host into device array | no-op | +| `mirror.to_host()` | copies device into host array **in place** | no-op | +| `mirror.zero()` | zeroes the device array | zeroes the host array | +| `mirror.rebind(new_host)` | follows a reallocation by the owner | same | + +`to_host()` writes into the existing host array, so the owning library keeps +its reference and sees the new values. Every method returns the mirror, so +calls chain. If the owner reallocated the host array with a new shape or dtype +and `rebind()` was not called, the mirror raises `ValueError` instead of copying +into the wrong buffer. + +## A complete charge deposition + +```python +import numpy as np + +import cunumpy as xp + +DEPOSIT = r""" +#include +#include + +extern "C" __global__ +void deposit(const double* x, const double* w, double* rho, + long long n_particles, long long n_cells, double dx) { + CUNUMPY_THREAD_1D(ip, n_particles); + long long cell = (long long)(x[ip] / dx); + if (cell >= 0 && cell < n_cells) cunumpy_atomic_add(&rho[cell], w[ip] / dx); +} +""" + + +def deposit_host(x, w, rho, n_particles, n_cells, dx): + cells = (x[:n_particles] / dx).astype(np.int64) + inside = (cells >= 0) & (cells < n_cells) + np.add.at(rho, cells[inside], w[:n_particles][inside] / dx) + + +deposit = xp.Kernel(deposit_host, xp.CudaKernel(DEPOSIT, "deposit")) + +# owned by a host library, e.g. a distributed vector +rho_host = np.zeros(64) +rho = xp.DeviceMirror(rho_host) + +rng = np.random.default_rng(0) +x = xp.to_cunumpy(rng.uniform(0.0, 1.0, 10_000)) +w = xp.to_cunumpy(np.full(10_000, 1e-4)) + +rho.zero() +deposit(x, w, rho.device, x.size, rho_host.size, 1.0 / 64, n_threads=x.size) +rho.to_host() # rho_host now holds the charge density, on both backends + +# the host library continues with rho_host, e.g. an MPI Allreduce +``` + +The same lines run on the NumPy backend, where `rho.device is rho_host`, the +host kernel writes into it directly, and `to_host()` does nothing. + +## Guidelines + +* Create the mirror once and keep it with the owner of the buffer. Allocation + happens on the first `device` access. +* Make exactly one `to_host()` per accumulation, after all kernels have written. + The mirror's copies are explicit method calls, so they are easy to find in the + code; note that only the first allocation of `device` (through `to_cupy()`) is + recorded by `count_transfers()`, while `to_host()` and `to_device()` are not. +* `zero()` before accumulating unless the kernel should add to the previous + contents; then call `to_device()` first if the host side changed. +* Atomics on one hot cell serialize. If most particles hit few cells, consider + sorting particles by cell or accumulating per block in shared memory first. +* Floating-point atomics make the summation order non-deterministic. Results + differ between runs in the last bits; compare with a tolerance in tests. diff --git a/docs/source/kernels/arguments.md b/docs/source/kernels/arguments.md new file mode 100644 index 0000000..6dce885 --- /dev/null +++ b/docs/source/kernels/arguments.md @@ -0,0 +1,211 @@ +# Kernel arguments and structs + +Kernels of real simulations take many arguments: the marker array, its shape, +the grid spacing, spline degrees, knot vectors, a dozen parameters. Listing +them at every call site is error-prone, and every new field means editing every +kernel signature. CuNumpy offers three ways to pass a group of values as one +argument: + +| Tool | Host kernel sees | CUDA kernel sees | Use when | +| --- | --- | --- | --- | +| `CudaArguments` | (not used) | several parameters, flattened | grouping device arguments only | +| `KernelArguments` | one object (`__host_args__()`) | several parameters (`__cuda_args__()`) | the host kernel takes an argument class, the CUDA kernel flat parameters | +| `CudaStruct` | (not used) | one C struct parameter | many fields; one definition shared by all CUDA kernels | + +They combine: a `KernelArguments` object can return a `CudaStruct` value from +`__cuda_args__()`. + +## `CudaArguments`: flatten into several parameters + +Any object with a `__cuda_args__()` method returning a tuple is replaced by +that tuple when passed to a `CudaKernel`. Subclassing `CudaArguments` gives you +the method for free: + +```python +import numpy as np + +import cunumpy as xp + + +class DeviceParticles(xp.CudaArguments): + def __init__(self, positions, velocities): + self.positions = xp.as_device_array(positions, np.float64, ndim=2, name="positions") + self.velocities = xp.as_device_array(velocities, np.float64, ndim=2, name="velocities") + super().__init__(self.positions, self.velocities, self.positions.shape[0]) + + +PUSH = r""" +extern "C" __global__ +void push(double dt, double* x, const double* v, int n) { ... } +""" +push = xp.CudaKernel(PUSH, "push") +particles = DeviceParticles(x, v) +push(0.1, particles, n_threads=particles.positions.shape[0]) # -> push(0.1, x, v, n) +``` + +Build the object once and reuse it. `as_device_array()` references existing +device arrays of the right dtype and copies everything else exactly once (see +[Data movement](../guides/data-movement.md), section "Build device arguments once"). + +## `KernelArguments`: one object, two forms + +A host kernel compiled with Pyccel often receives a group of arrays as an +instance of an argument class, while the CUDA kernel receives the arrays as +separate parameters. `KernelArguments` lets one object stand for both, so the +call site never branches on the backend: + +```python +class MarkerArguments: # the host argument class, e.g. compiled with Pyccel + def __init__(self, markers: "float[:, :]", n_markers: int): + self.markers = markers + self.n_markers = n_markers + + +class ParticleArguments(xp.KernelArguments): + def __init__(self, particles): + self._particles = particles + self._host = None + self._cuda = None + + def __host_args__(self): + if self._host is None: + markers = self._particles.markers + self._host = MarkerArguments(markers, markers.shape[0]) + return self._host + + def __cuda_args__(self): + if self._cuda is None: + markers = self._particles.markers + self._cuda = (markers, markers.shape[0], markers.shape[1]) + return self._cuda + + +args = ParticleArguments(particles) +push(args, dt, n_threads=n_markers) # a Kernel: same call on both backends +``` + +* `Kernel` (on the NumPy backend) and `PyccelKernel` replace the argument with + `__host_args__()`; `CudaKernel` flattens `__cuda_args__()`. +* Only top-level arguments are resolved, not objects inside lists or dicts. +* Build both forms lazily, as above: a CPU run never builds device arguments, a + GPU run never builds the host object. +* **Invalidate the cache** when the underlying arrays are replaced (resizing, + `deepcopy`, unpickling): reset the stored forms or create a new arguments + object. Stale cached arguments point at the old arrays. +* `xp.resolve_host_args(args, kwargs)` applies the host replacement, for code + that calls host kernels without `Kernel`. + +## `CudaStruct`: one C struct, defined once + +A `CudaStruct` defines a C struct in Python: its C declaration for the kernel +source, its exact memory layout, and a packer for values. The kernel takes the +struct by value as one parameter: + +```python +Particles = xp.CudaStruct( + "Particles", + [("x", "double*"), ("v", "double*"), ("n", "long long"), ("charge", "double")], +) + +PUSH = Particles.declaration + r""" +extern "C" __global__ void push(Particles p, double dt) { + long long i = blockDim.x * (long long)blockIdx.x + threadIdx.x; + if (i < p.n) p.x[i] += dt * p.charge * p.v[i]; +} +""" +push = xp.CudaKernel(PUSH, "push", structs=[Particles]) + +value = Particles(x=x, v=v, n=x.size, charge=-1.0) +push(value, 0.1, n_threads=x.size) +``` + +Fields may be scalars, pointers to scalar types (or `void*`), and array views +`Array1D` to `Array3D`. 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. + +The packed value (`CudaStructValue`) holds only device *addresses*. It keeps +references to the arrays so they stay alive, but if you replace an array (for +example `particles.x = new_array`), pack a new value. `value["charge"]` reads a +field back; `Particles.dtype` is the NumPy structured dtype of the layout. + +`CudaKernel(..., structs=[Particles])` also checks that a struct definition in +the kernel source matches the Python definition, so a hand-edited header that +drifted raises a `ValueError` instead of reading fields at wrong offsets. + +### Generate the struct from the host argument class + +If the host kernels already use an annotated argument class (Pyccel style), the +struct can be derived from it, making the Python class the single definition of +the arguments on host and device: + +```python +class MarkerArguments: + def __init__(self, markers: "float[:, :]", n_markers: int, valid: "bool[:]"): + ... + + +MarkerArgs = xp.CudaStruct.from_signature(MarkerArguments.__init__, "MarkerArgs") +print(MarkerArgs.declaration) +``` + +```c +struct MarkerArgs { + Array2D markers; + long long n_markers; + Array1D valid; +}; +``` + +Mappings: `float` to `double`, `int` to `long long` (Pyccel integers are +64-bit; change with `int_type=`), `bool` to `bool`, NumPy scalar types to their +C types, `"float[:, :]"` to `Array2D` (1 to 3 dimensions). `Final[...]` +and `const` are ignored. `scalar_names={"float": "float"}` switches to single +precision. Parameters without a mappable annotation raise `ValueError`. + +Array view fields accept non-contiguous arrays, and the kernel indexes them like +the host kernel does: + +```c +#include "marker_args.cuh" +#include + +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); +} +``` + +### Write the struct to a header + +Kernels in `.cu` files include the struct from a header. Generate it from the +Python definition and commit it: + +```python +xp.write_cuda_header("kernels/marker_args.cuh", [MarkerArgs, DomainArgs]) +``` + +The header gets an include guard (`MARKER_ARGS_CUH`), the `array_view.cuh` +include when needed, and the definitions. `MarkerArgs.to_header(path)` does the +same for one struct. A test keeps the committed file in sync: + +```python +from pathlib import Path + + +def test_marker_args_header_is_up_to_date(tmp_path): + generated = xp.write_cuda_header(tmp_path / "marker_args.cuh", [MarkerArgs, DomainArgs]) + assert Path("kernels/marker_args.cuh").read_text() == generated +``` + +## Choosing + +* Start with plain arguments. Group when the same set of five or more values + appears in several kernels. +* Use `KernelArguments` when the host kernels already take argument objects; + it keeps the call sites identical. +* Use a `CudaStruct` when many CUDA kernels take the same group, or when + signatures become too long to read and keep in sync. +* Generate the struct with `from_signature` when a host argument class exists, + so the two cannot drift apart. diff --git a/docs/source/kernels/cuda-kernel.md b/docs/source/kernels/cuda-kernel.md new file mode 100644 index 0000000..195789a --- /dev/null +++ b/docs/source/kernels/cuda-kernel.md @@ -0,0 +1,266 @@ +# Writing CUDA kernels + +`CudaKernel` wraps a CUDA C kernel so it can be called from Python like the +host kernel it mirrors, plus a launch shape. It builds on CuPy's `RawKernel` +(compiled at run time with NVRTC, no `nvcc` or build step needed) and adds what +`RawKernel` lacks: signature checks, launch-shape arithmetic, header handling, +and debug support. + +## A first kernel + +```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]; +} +""" + +axpy = xp.CudaKernel(AXPY, "axpy") + +xp.set_backend("cupy") +x = xp.arange(10_000, dtype=xp.float64) +y = xp.zeros_like(x) +axpy(2.0, x, y, x.size, n_threads=x.size) +``` + +The essentials: + +* The kernel is declared `extern "C" __global__` so its name is not mangled + (templates are the exception, see below). +* Arguments are passed in the order of the C signature. `n_threads` is the + number of threads to launch; CuNumpy computes the grid as + `ceil(n_threads / block_size)`. Threads beyond `n` must return early, as + above, because the last block is usually only partially used. +* Compilation happens on the first call, or when `axpy.compile()` is called. + CuPy caches the binary on disk (`~/.cupy/kernel_cache`), so later runs load it. +* Creating a `CudaKernel` does not import CuPy. Kernels can be defined at module + level in code that also runs on machines without a GPU. + +## What the signature checks catch + +`RawKernel` reads each argument with the size declared in C and never checks +it: a Python `int` passed to a `double` parameter, a `float64` array passed as +`float*`, or a NumPy array instead of a CuPy array produce wrong numbers or +crashes without an error. `CudaKernel` parses the signature once and checks +every call: + +| Parameter type | Accepted | Rejected | +| --- | --- | --- | +| `double*`, `const int*`, ... | C-contiguous CuPy array of exactly that dtype | NumPy arrays (`TypeError`, never copied), other dtypes, non-contiguous views | +| `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` in range, `bool`, matching NumPy integers | out-of-range values (`OverflowError`), floats | +| `Array1D` ... `Array3D` | 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 +`long long` are `int64`, `float` is `float32`, `double` is `float64`. +`xp.ctype_of(np.float64)` returns `"double"`, useful when generating source. + +A wrong argument count, dtype or layout raises before anything is launched, +with the parameter name in the message: + +```python +axpy(2, x, y, x.size, n_threads=x.size) # fine: 2 is cast to double +axpy(2.0, xp.to_numpy(x), y, x.size, n_threads=x.size) +# TypeError: argument 1 (double* x) must be a CuPy array, got ndarray; arrays are never copied to the device +``` + +Non-contiguous views such as `a[:, 0]` are rejected for pointer parameters +because the kernel would read them as a flat buffer. Make them contiguous with +`xp.ascontiguousarray()` (a copy), or use an array view parameter (below), +which carries the strides. + +The checks cost about 0.3 µs per argument. For tiny kernels called millions of +times, `check_signature=False` disables them once the calls are known to be +correct; the kernel then behaves like a raw `RawKernel`. + +## Launch configuration + +```python +kernel(*args, n_threads=None, grid=None, block=None, shared_mem=0, stream=None) +``` + +* **1D**: `n_threads=n` with the default `block_size=128` (set per kernel with + `CudaKernel(..., block_size=256)`). +* **2D/3D**: `n_threads=(nx, ny)` and a tuple block, e.g. + `CudaKernel(src, "stencil", block_size=(16, 16))`. In the kernel, use + `blockIdx.y`/`threadIdx.y` for the second dimension. +* **Explicit grid**: `grid=(n_blocks,)` instead of `n_threads`, for kernels that + loop internally (grid-stride loops) or need a specific number of blocks. +* **Per-call block**: `block=(32, 8)` overrides `block_size` for one call. +* **Dynamic shared memory**: `shared_mem` bytes per block for + `extern __shared__` arrays. +* **Stream**: `stream=s` queues the launch on a stream, see [Devices, memory and + streams](../guides/gpu-devices.md). + +A zero-sized launch (`n_threads=0`) launches nothing. `launch_shape()` returns +the `(grid, block)` a call would use, which is how to size per-block outputs: + +```python +BLOCK_SUM = r""" +extern "C" __global__ void block_sum(const double* x, double* out, int n) { + extern __shared__ double buffer[]; + int i = blockDim.x * blockIdx.x + threadIdx.x; + buffer[threadIdx.x] = i < n ? x[i] : 0.0; + __syncthreads(); + for (int s = blockDim.x / 2; s > 0; s /= 2) { + if (threadIdx.x < s) buffer[threadIdx.x] += buffer[threadIdx.x + s]; + __syncthreads(); + } + if (threadIdx.x == 0) out[blockIdx.x] = buffer[0]; +} +""" +block_sum = xp.CudaKernel(BLOCK_SUM, "block_sum", block_size=128) +(n_blocks,), (threads,) = block_sum.launch_shape(x.size) +partial = xp.zeros(n_blocks) +block_sum(x, partial, x.size, n_threads=x.size, shared_mem=threads * 8) +total = float(partial.sum()) +``` + +## Kernels in files + +Keeping CUDA source in `.cu` files gives editor support and lets kernels share +headers: + +```python +push = xp.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` +suffix (override with `name=` or `suffix=`), and adds the file's directory to the +include path, so `#include "helpers.cuh"` next to it works. + +A file with several small kernels is loaded at once with `all_from_file`, which +returns a dict by name; the kernels share one compilation: + +```python +ops = xp.CudaKernel.all_from_file("kernels/vector_ops.cu", block_size=256) +ops["scale"](x, 2.0, x.size, n_threads=x.size) +ops["shift"](x, 1.0, x.size, n_threads=x.size) +``` + +`xp.cuda_kernel_names(source)` lists the `__global__` functions of a source. + +## Headers and the compile cache + +Pass extra include directories with `include_dirs=[...]` and NVRTC flags with +`options=("-std=c++17",)`. The headers shipped with CuNumpy are always found: + +| Header | Provides | +| --- | --- | +| `` | `CUNUMPY_THREAD_1D(i, n)`, `_2D`, `_3D`, `CUNUMPY_GRID_STRIDE_1D(i, n)` | +| `` | strided views `Array1D`, `Array2D`, `Array3D` | +| `` | `cunumpy_atomic_add` and indexed 2D/3D variants, see [Accumulation kernels](accumulation.md) | + +`xp.cuda_include_dir()` returns their directory for use with other compilers. + +CuPy's disk cache is keyed on the source string and the options only, so +editing an included header would normally *not* trigger a recompile. +`CudaKernel` resolves the `#include "..."` files of a source (recursively), +hashes their contents, and adds `-DCUNUMPY_INCLUDE_HASH=0x...` to the options, +so a changed header produces a new cache entry. `kernel.included_headers` lists +the resolved files and `kernel.compile_options()` the options actually used. +System includes in angle brackets are not hashed. + +## Index macros and array views + +Kernels ported from Python index multi-dimensional arrays like `markers[ip, 3]`. +With bare pointers, the CUDA version has to receive the shape and compute +`markers[ip * n_cols + 3]` by hand, and breaks for non-contiguous arrays. Array +views carry shape and strides instead: + +```python +SCALE_COLUMN = r""" +#include +#include + +extern "C" __global__ +void scale_column(Array2D a, long long column, double factor) { + CUNUMPY_THREAD_1D(i, a.shape[0]); // long long i; returns if i >= shape[0] + a(i, column) *= factor; +} +""" +scale_column = xp.CudaKernel(SCALE_COLUMN, "scale_column") + +markers = xp.zeros((1000, 7)) +view = markers[::2, 1:5] # non-contiguous view is fine +scale_column(view, 1, 10.0, n_threads=view.shape[0]) +``` + +An `ArrayND` parameter takes a CuPy array of dtype `T` and `N` dimensions; it +is passed by value as pointer, shape and strides (in elements). `a(i, j)` +returns a reference, `a.shape[k]` and `a.size()` give the extents. Compiling +with `-DCUNUMPY_BOUNDS_CHECK` (automatic in [debug mode](debugging.md)) checks +every index against the shape and traps with a message on violation. + +The index macros remove the boilerplate of thread index computation: + +```c +CUNUMPY_THREAD_1D(i, n); // long long i; return if i >= n +CUNUMPY_THREAD_2D(i, j, ni, nj); // 2D launch, n_threads=(ni, nj) +CUNUMPY_THREAD_3D(i, j, k, ni, nj, nk); +CUNUMPY_GRID_STRIDE_1D(i, n) { ... } // loop; launch with any grid +``` + +## Templates and generated kernels + +C++ function templates are instantiated with `template_args`; dtypes are +converted with `ctype_of`, and the instantiated signature is checked as usual: + +```python +import numpy as np + +SCALE = r""" +template +__global__ void scale(T* x, T factor, int n) { + int i = blockDim.x * blockIdx.x + threadIdx.x; + if (i < n) x[i] *= factor; +} +""" +scale_f64 = xp.CudaKernel(SCALE, "scale", template_args=(np.float64,)) +scale_f32 = xp.CudaKernel(SCALE, "scale", template_args=(np.float32,)) +``` + +When the source itself is generated per variant (unrolled loops per dimension, +per spline degree, per dtype), `CudaKernelVariants` creates one kernel per key +on first use and caches it: + +```python +def make_matvec(ndim, dtype): + return xp.CudaKernel(generate_source(ndim, xp.ctype_of(dtype)), "matvec") + + +matvec = xp.CudaKernelVariants(make_matvec) +matvec.get(3, np.float64)(mat, x, out, n_threads=out.size) +matvec.compile_all([(3, np.float64), (3, np.complex128)], jobs=4) # at setup +``` + +## Compile at setup + +The first call of each kernel compiles it, which takes from a fraction of a +second to several seconds. Call `kernel.compile()` (or +`catalog.compile_all(jobs=8)` for a whole catalog, in parallel threads) during +setup so the first time step is not slower than the rest and compilation errors +appear before the simulation starts. `kernel.is_compiled` tells whether it has +happened. Compiling requires a GPU (`RuntimeError` otherwise). + +## Tips for writing kernels + +* **Mirror the host kernel's argument list.** Same order, same names. Then the + call site is identical on both backends (see [Pairing host and CUDA + kernels](dispatch.md)) and parity tests can share argument builders. +* **Use `long long` for indices into large arrays.** `int` overflows above + 2^31 elements; the index macros already declare `long long`. +* **Guard the tail.** Every 1D kernel needs `if (i >= n) return;` (or + `CUNUMPY_THREAD_1D`). +* **Keep `const` on read-only pointers.** It documents intent and lets the + compiler use the read-only cache. +* **Many threads writing one location need atomics.** Use `cunumpy_atomic_add` + ([Accumulation kernels](accumulation.md)). +* **When something crashes, enable debug mode** before anything else ([Debugging + CUDA kernels](debugging.md)). diff --git a/docs/source/kernels/debugging.md b/docs/source/kernels/debugging.md new file mode 100644 index 0000000..e952e69 --- /dev/null +++ b/docs/source/kernels/debugging.md @@ -0,0 +1,95 @@ +# Debugging CUDA kernels + +Kernel launches are asynchronous. When a kernel reads out of bounds, the error +is reported by the *next* operation that waits for the GPU (a `.get()`, an MPI +call, the next allocation), often far from the kernel that caused it, with a +message like `cudaErrorIllegalAddress: an illegal memory access was +encountered`. Out-of-bounds accesses that stay inside allocated memory produce +no error at all, only wrong numbers. + +## Step 1: debug mode + +Enable debug mode in any of these ways: + +```bash +CUNUMPY_CUDA_DEBUG=1 python simulate.py # whole process +``` + +```python +xp.set_cuda_debug(True) # globally, from now on +with xp.cuda_debug(): # for a block + ... +xp.CudaKernel(src, "push", debug=True) # one kernel, regardless of the global setting +``` + +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`/`Array2D`/`Array3D` 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 + the CUDA error chained. + +```python +with xp.cuda_debug(): + push = xp.CudaKernel(SOURCE, "push") + push(markers, dt, n, n_threads=n) +# RuntimeError: CUDA error after launching kernel 'push' with grid (79,) and block (128,): ... +``` + +Compile options are fixed when a kernel is compiled. A kernel that was already +compiled before debug mode was enabled keeps its options (the synchronization +still applies). Enable debug mode before the kernels are first called, or set +`CUNUMPY_CUDA_DEBUG=1` in the environment. `kernel.debug_active()` and +`kernel.compile_options()` show what applies to a kernel. + +Use `-DCUNUMPY_BOUNDS_CHECK` in your own code too: + +```c +#ifdef CUNUMPY_BOUNDS_CHECK + if (cell < 0 || cell >= n_cells) { printf("bad cell %lld\n", cell); __trap(); } +#endif +``` + +`-G` (full device debug symbols) is not available with NVRTC. + +## Step 2: compute-sanitizer + +Debug mode tells you *which* kernel failed. NVIDIA's memory checker tells you +*where*, thanks to `-lineinfo`, and also finds out-of-bounds accesses that do +not crash: + +```bash +CUNUMPY_CUDA_DEBUG=1 compute-sanitizer python -m pytest tests/test_push.py +CUNUMPY_CUDA_DEBUG=1 compute-sanitizer --tool racecheck python check_deposit.py +``` + +The `memcheck` tool (default) finds invalid reads and writes, `racecheck` +shared-memory races, `initcheck` reads of uninitialized device memory. Run them +on a small problem; they slow execution down by a large factor. + +## Step 3: compare with the host kernel + +Most porting bugs are not crashes but differences: an index off by one, a +wrong stride, a missing `if`. The host kernel is the reference. +`assert_kernels_agree()` runs both on the same inputs and names the argument +that differs (see [Testing kernels](testing.md)). For `__device__` helpers, +`device_function_kernel()` exposes a single function to Python so it can be +compared value by value with its host version. + +## Common causes + +| Symptom | Likely cause | +| --- | --- | +| `TypeError` at the call, before launch | argument does not match the C signature: host array, wrong dtype, non-contiguous view, wrong count. Read the message; it names the parameter. | +| illegal memory access | index beyond the array: missing `if (i >= n) return;`, wrong `n`, wrong stride arithmetic. Use array views with bounds checks. | +| correct on small inputs, wrong on large ones | `int` overflow in index computations; use `long long`. | +| results differ slightly between runs | floating-point atomics in a different order; compare with a tolerance. | +| results occasionally garbage after MPI | missing `synchronize_for_mpi()` before the MPI call. | +| a header change has no effect | the header is included with angle brackets (not hashed) or the kernel object was compiled before the change; see "Headers and the compile cache" in [Writing CUDA kernels](cuda-kernel.md). | + +After an illegal memory access the CUDA context of the process is unusable: all +later CUDA calls fail. Restart the process (or the pytest run) after fixing the +kernel. diff --git a/docs/source/kernels/dispatch.md b/docs/source/kernels/dispatch.md new file mode 100644 index 0000000..43d541f --- /dev/null +++ b/docs/source/kernels/dispatch.md @@ -0,0 +1,160 @@ +# Pairing host and CUDA kernels + +`Kernel` holds a host kernel and, once it is written, its CUDA counterpart, and +calls the one matching the active backend. `KernelCatalog` collects all kernels +of a package. Together they let a code base port its kernels incrementally while +every call site stays the same. + +## `Kernel` + +```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_host(a, x, y, n): + for i in range(n): + y[i] += a * x[i] + + +axpy = xp.Kernel(axpy_host, xp.CudaKernel(AXPY, "axpy"), name="axpy") + +axpy(2.0, x, y, x.size, n_threads=x.size) +``` + +* On the **NumPy backend** the host kernel is called with the positional + arguments; `n_threads`, `grid`, `block`, `shared_mem` and `stream` are + ignored. +* On the **CuPy backend** the CUDA kernel is launched; `n_threads` (or `grid`) + is required. +* A plain host function is wrapped in a [`PyccelKernel`](pyccel-kernel.md); + pass `host_options={"outputs": (2,)}` to configure that wrapper, or pass a + `PyccelKernel` you built yourself. +* Kernels take positional arguments only. [`KernelArguments`](arguments.md) + objects are resolved per backend. +* `kernel.compile()` compiles the CUDA kernel now; `kernel.has_cuda` tells + whether there is one. + +### Kernels without CUDA version + +`Kernel(host, None)` is a kernel not yet ported. On the CuPy backend: + +* `missing_cuda="raise"` (default) raises `NotImplementedError`, naming + `cuda_path` if given, the file where the CUDA kernel is expected. +* `missing_cuda="fallback"` calls the host kernel through `PyccelKernel`, + copying the arrays to the host and back at every call. A `RuntimeWarning` is + emitted once per kernel, and `count_transfers()` records a `fallback` event + per call. + +For the fallback to find device arrays inside your own objects, the wrapper +needs to know them: `host_options={"object_modules": ("my_package.",)}`. + +`kernel.get_kernel()` returns the kernel for the active backend. Calling it +once at setup surfaces a missing CUDA kernel immediately instead of in the +middle of a run. + +## `KernelCatalog`: one folder per kernel + +`KernelCatalog.from_package()` finds kernels by convention: + +```text +my_sim/kernels/ +├── __init__.py # catalog = xp.KernelCatalog.from_package(__name__) +├── push/ +│ ├── push_kernels.py # def push(...): ... host kernel +│ └── push_cuda.cu # __global__ void push(...) CUDA kernel +├── deposit/ +│ ├── deposit_kernels.py +│ └── deposit_cuda.cu +└── sort/ + └── sort_kernels.py # not ported yet +``` + +```python +# my_sim/kernels/__init__.py +import cunumpy as xp + +catalog = xp.KernelCatalog.from_package(__name__, missing_cuda="fallback") +``` + +```python +# anywhere in the code +from my_sim.kernels import catalog + +catalog["push"](markers, dt, n_markers, n_threads=n_markers) +``` + +For every subfolder `` containing `_kernels.py`, the function +`` in that module is the host kernel, and `_cuda.cu` in the same +folder, if present, provides the CUDA kernel `__global__ void (...)`. +Other `__global__` functions in that file are ignored by the catalog (load them +with `CudaKernel.all_from_file`). The suffixes are configurable with +`host_suffix` and `cuda_suffix`. + +Options: + +* `missing_cuda`: passed to every `Kernel`. +* `host_options`: `PyccelKernel` options for all host kernels, or a function of + the kernel name for per-kernel output declarations: + + ```python + OUTPUTS = {"push": (0,), "deposit": (2,)} + catalog = xp.KernelCatalog.from_package( + __name__, + host_options=lambda name: {"outputs": OUTPUTS.get(name)}, + ) + ``` + +* `include_dirs`: extra include directories for the CUDA kernels. By default + the source root of the top-level package is one, so a kernel can write + `#include "my_sim/kernels/common.cuh"`. Each kernel's own folder is always + searched. +* Other keyword arguments (`block_size`, `structs`, `options`, ...) go to every + `CudaKernel.from_file()`. + +A catalog is a read-only mapping: `catalog["push"]`, `"push" in catalog`, +`len(catalog)`, iteration over names. `KernelCatalog()` plus +`catalog.register(kernel)` builds one by hand. + +## Porting status and setup + +```python +print(catalog.summary()) +# CUDA kernels: 2 of 3 (missing: sort) + +catalog.without_cuda # ['sort'] +catalog.with_cuda # ['deposit', 'push'] +``` + +`summary()` is handy for a `--status` command line flag or the start-up log. + +At setup, on the GPU backend, compile everything at once: + +```python +if xp.cupy_backend: + catalog.compile_all(jobs=8) # threads; jobs=None uses all CPUs +``` + +All kernels are compiled even if one fails; the first error is raised +afterwards. On later runs CuPy loads the binaries from its disk cache, so this +is fast. + +## Recommended project layout + +* One folder per kernel, the host kernel and its CUDA port side by side, so a + reviewer sees both at once. +* Shared device helpers (`__device__` functions) in `.cuh` headers next to the + kernels or in a `common/` folder, included with quotes. Test them with + `device_function_kernel` ([Testing kernels](testing.md)). +* Generated struct headers (`write_cuda_header`) committed next to the kernels, + with a test that they are up to date. +* One parametrised parity test over `catalog.parity_cases()`. +* `missing_cuda="fallback"` while porting, `"raise"` once the time loop is fully + ported. diff --git a/docs/source/kernels/overview.md b/docs/source/kernels/overview.md new file mode 100644 index 0000000..3190b9b --- /dev/null +++ b/docs/source/kernels/overview.md @@ -0,0 +1,78 @@ +# Porting kernels to the GPU + +Array-level code (`xp.sum`, `xp.fft`, slicing, broadcasting) runs on the GPU as +soon as the backend is CuPy. Many scientific codes, however, spend their time in +*kernels*: loops over particles or grid cells written in Python and compiled +for the CPU, for example with [Pyccel](https://github.com/pyccel/pyccel) or +Numba. Those kernels take NumPy arrays and cannot run on device memory. + +CuNumpy's kernel layer lets such a code base move to the GPU **one kernel at a +time**, with the CPU version kept as the reference and the call sites +unchanged. + +## The building blocks + +| Class | What it does | Guide | +| --- | --- | --- | +| `PyccelKernel` | calls a host kernel with CuPy arrays by copying them to the host and back | [Host kernels with GPU data](pyccel-kernel.md) | +| `CudaKernel` | wraps a CUDA C kernel, checks every call against its signature | [Writing CUDA kernels](cuda-kernel.md) | +| `Kernel` | a host kernel plus its CUDA kernel; calls the one matching the backend | [Pairing host and CUDA kernels](dispatch.md) | +| `KernelCatalog` | all `Kernel`s of a package, found by folder convention | [Pairing host and CUDA kernels](dispatch.md) | +| `CudaArguments`, `KernelArguments`, `CudaStruct` | pass a group of arrays and scalars as one argument | [Kernel arguments and structs](arguments.md) | +| `DeviceMirror`, `cunumpy/atomic.cuh` | scatter-add into a buffer owned by a host library | [Accumulation kernels](accumulation.md) | +| `cunumpy.testing` | check that host and CUDA kernels compute the same | [Testing kernels](testing.md) | + +## Which one do I need? + +* **"I just want my existing code to run with `ARRAY_BACKEND=cupy`."** Wrap the + host kernels in `PyccelKernel`. Everything works, but every kernel call + copies its arrays to the host and back. This is a correct starting point, not + a fast one. +* **"I have one hot kernel and want it fast on the GPU."** Write a CUDA version + and launch it with `CudaKernel`. +* **"I have dozens of kernels and will port them over months."** Put each + host/CUDA pair in a `Kernel`, collect them in a `KernelCatalog`, and test each + pair with `assert_kernels_agree`. Unported kernels either raise or fall back + to the host version. +* **"My kernels take ten arrays each."** Group them with `KernelArguments` (host + form and device form of the same data) or a `CudaStruct`. + +## A typical porting workflow + +1. **Run everything on the host.** The code base uses `import cunumpy as xp` + and the NumPy backend; kernels are plain host functions (Pyccel, Numba, or + pure Python). +2. **Make the GPU backend work, slowly.** Collect the kernels in a + `KernelCatalog` with `missing_cuda="fallback"`. On CuPy every kernel now runs + through `PyccelKernel`, with host copies. Results must match the CPU run. + `count_transfers()` shows how many copies each step makes; `catalog.summary()` + shows how many kernels are left (`CUDA kernels: 0 of 42`). +3. **Port the most expensive kernel.** Write `name/name_cuda.cu` next to the + host kernel, mirroring its argument list. The catalog picks it up + automatically. A parity test compares it with the host version. +4. **Repeat**, ordered by a profile (see [Timing and + profiling](../guides/profiling.md)), until the kernels in the time loop are + ported. `assert_no_transfers()` around a time step verifies that no fallback + or copy is left. +5. **Switch to `missing_cuda="raise"`** so a kernel added later without CUDA + version fails loudly instead of silently copying. + +The [particle pusher example](../examples/particle-pusher.md) walks through +these steps on a small code. + +## Design rules the kernel layer follows + +Knowing these makes the error messages predictable: + +* **CUDA kernels never copy.** `CudaKernel` accepts only device arrays for + pointer parameters. A NumPy array raises `TypeError` instead of being copied + to the device behind your back. +* **Signatures are checked.** CuPy's `RawKernel` reads arguments with the size + declared in C and never checks them, so a Python `int` passed as `double` + arrives as garbage. `CudaKernel` parses the signature and casts or rejects + every argument. +* **Copies of the fallback path are visible.** `PyccelKernel` conversions and + `Kernel` fallbacks are recorded by `count_transfers()` and warned about once. +* **The CUDA kernel mirrors the host kernel.** Same name, same argument order, + plus a launch shape (`n_threads`). The call site does not branch on the + backend. diff --git a/docs/source/kernels/pyccel-kernel.md b/docs/source/kernels/pyccel-kernel.md new file mode 100644 index 0000000..8a9f319 --- /dev/null +++ b/docs/source/kernels/pyccel-kernel.md @@ -0,0 +1,117 @@ +# Host kernels with GPU data + +`PyccelKernel` adapts a function that expects NumPy arrays so it can be called +with CuPy arrays. It is named after Pyccel, whose compiled functions accept +only NumPy arrays, but it wraps any callable: a Numba function, a C extension, +or plain Python. It does not compile anything. + +## What it does + +```python +import cunumpy as xp + + +def smooth(field, out): + out[1:-1] = 0.25 * field[:-2] + 0.5 * field[1:-1] + 0.25 * field[2:] + out[0], out[-1] = field[0], field[-1] + + +smooth_kernel = xp.PyccelKernel(smooth, outputs=(1,)) +``` + +On every call it decides whether conversion is needed (the active backend is +CuPy, or a CuPy array is among the arguments): + +* **No conversion** (NumPy backend, NumPy arrays): the function is called + directly. Identity, mutation, return values and exceptions behave as without + the wrapper; the overhead is a type check per argument. +* **Conversion**: every CuPy array in the arguments is copied to the host, the + function runs on the host copies, the arrays it may have written are copied + back into the original CuPy arrays, and returned NumPy arrays are converted to + CuPy. + +```python +with xp.use_backend("cupy"): + field = xp.random.random(1000) + out = xp.empty_like(field) + smooth_kernel(field, out) # out is updated on the device +``` + +## Declare the outputs + +The wrapper cannot know which arguments the function writes. By default it +copies back *every* converted array, which is correct but doubles the +transfers. `outputs` lists the arguments that may be written: + +```python +xp.PyccelKernel(smooth, outputs=(1,)) # positional argument 1 +xp.PyccelKernel(update, outputs=(0, -1)) # first and last argument +xp.PyccelKernel(solve, outputs=("out",)) # solve(a, b, out=out) +xp.PyccelKernel(norm, outputs=()) # writes nothing +``` + +* Positional arguments are declared by index (negative indices count from the + end), keyword arguments by name. The two forms are not interchangeable, + because compiled functions usually do not expose a Python signature that + would map names to positions. +* A container or object declared as output has all its arrays copied back. +* **A missing declaration is a silent bug**: if the function writes an + argument that is not declared, the device array keeps its old values. When in + doubt, leave `outputs=None`. + +## Arrays inside containers and objects + +Lists, tuples and dicts are traversed recursively, so `kernel([x, y], {"v": v})` +works. Instances of your own classes are traversed only if their class's module +starts with one of the `object_modules` prefixes: + +```python +kernel = xp.PyccelKernel(push, object_modules=("my_simulation.",), outputs=(0,)) +kernel(particles, dt) # particles.positions etc. are converted +``` + +Such objects are shallow-copied and their array attributes replaced on the +copy, so the caller's object keeps referencing the device arrays. Objects of +other modules are passed through unchanged. + +Other details: + +* **Aliasing is preserved.** If the same device array appears twice, the host + function sees one host array twice. +* **Return values**: NumPy arrays (also inside returned tuples and lists) are + converted to CuPy; `is_array` customizes which return values count as arrays + (for NumPy subclasses). +* **`use_cupy`** forces conversion on (`True`) or off (`False`); the default + `None` decides per call. +* **Argument objects with two forms** (`KernelArguments`) are replaced by their + `__host_args__()` first; see [Kernel arguments and structs](arguments.md). + +## Costs and when to move on + +Each converted call copies every input array to the host and every output back. +For a kernel called every time step with large arrays, this dominates the run +time. `count_transfers()` records one `kernel_conversion` event per converted +call, naming the kernel and the number of arrays: + +```python +with xp.count_transfers() as counter: + step(state, dt) +for event in counter.kernel_conversion_calls: + print(event.where, event.description) +``` + +`PyccelKernel` is the right tool for: + +* getting a code base to run on the GPU backend before any kernel is ported; +* kernels outside the time loop (setup, I/O, rare diagnostics); +* the fallback of a `Kernel` that has no CUDA version yet. + +For kernels inside the time loop, write a CUDA version ([Writing CUDA +kernels](cuda-kernel.md)) and pair it with the host kernel ([Pairing host and +CUDA kernels](dispatch.md)). + +## Pyodide + +In Pyodide only the no-conversion path exists. Plain Python functions can still +be wrapped, so code that uses `PyccelKernel` runs in the browser; see +[Pyodide](../pyodide.md). diff --git a/docs/source/kernels/testing.md b/docs/source/kernels/testing.md new file mode 100644 index 0000000..145a168 --- /dev/null +++ b/docs/source/kernels/testing.md @@ -0,0 +1,204 @@ +# Testing kernels + +A GPU port is only as trustworthy as its comparison with the CPU version. +`cunumpy.testing` provides pytest helpers for exactly that, designed so that the +same test suite runs on a laptop without a GPU (GPU cases are skipped) and on a +GPU runner (everything runs). + +```python +from cunumpy.testing import ( + BACKENDS, + assert_kernels_agree, + backend, + device_function_kernel, + requires_cupy, +) +``` + +`cunumpy.testing` is not imported by `import cunumpy`, and it imports pytest only +when one of its pytest objects is used. + +## Run a test on both backends + +`BACKENDS` is `["numpy", pytest.param("cupy", marks=requires_cupy)]`: + +```python +import pytest + +import cunumpy as xp +from cunumpy.testing import BACKENDS + + +@pytest.mark.parametrize("backend", BACKENDS) +def test_norm(backend): + with xp.use_backend(backend): + assert float(xp.linalg.norm(xp.ones(4))) == 2.0 +``` + +The `backend` fixture does the same and also activates the backend for the +whole test. Import it into `conftest.py` to make it available everywhere: + +```python +# conftest.py +from cunumpy.testing import backend # noqa: F401 +``` + +```python +def test_energy_is_conserved(backend): + state = make_state() # arrays land on the active backend + e0 = energy(state) + for _ in range(100): + step(state, 1e-3) + assert abs(float(energy(state)) - float(e0)) < 1e-10 +``` + +`requires_cupy` is a plain skip marker for GPU-only tests: + +```python +from cunumpy.testing import requires_cupy + + +@requires_cupy +def test_kernel_compiles(): + with xp.use_backend("cupy"): + assert catalog["push"].compile() +``` + +## Compare host and CUDA kernels: `assert_kernels_agree` + +```python +import numpy as np + +import cunumpy as xp +from cunumpy.testing import assert_kernels_agree + + +def make_axpy_args(backend, seed): + rng = np.random.default_rng(seed) + x = xp.to_cunumpy(rng.random(1000)) + y = xp.to_cunumpy(rng.random(1000)) + return (2.0, x, y, 1000) + + +def test_axpy_parity(): + assert_kernels_agree(axpy, make_axpy_args, n_threads=1000) +``` + +For each backend, NumPy then CuPy, it activates the backend, builds the +arguments with `make_args(backend, seed)`, calls the kernel (`n_calls` times), +then copies the CUDA results to the host and compares them with +`numpy.testing.assert_allclose(rtol=1e-12, atol=0)`. A failure names the +argument that differs. Without a GPU the test is skipped. + +Things to know: + +* **Build random data on the host.** NumPy and CuPy generators produce + different sequences from the same seed, so use `numpy.random.default_rng` and + convert with `to_cunumpy()`, as above. +* **Which arguments are compared**: `outputs=(2,)` selects them by index; + otherwise the host kernel's declared `outputs` are used, and if there are + none, every array argument. Arrays held by argument objects (one level deep, + e.g. a `CudaArguments` object or a list) are compared too. +* **Tolerances**: the default `rtol=1e-12` suits deterministic kernels. Kernels + with atomics or a different summation order need looser tolerances, for + example `rtol=1e-10, atol=1e-14`. +* **`n_calls=10`** runs the kernel repeatedly on the same arguments, which + catches state that is not reset between calls. +* It returns the host results by argument name for additional assertions. + +### One test for the whole catalog + +```python +import pytest + +from my_sim.kernels import catalog +from my_sim.kernels.test_args import MAKE_ARGS, N_THREADS + + +@pytest.mark.parametrize("name, kernel", catalog.parity_cases()) +def test_parity(name, kernel): + assert_kernels_agree(kernel, MAKE_ARGS[name], n_threads=N_THREADS[name]) +``` + +`parity_cases()` yields the kernels that have a CUDA version, so every newly +ported kernel is tested as soon as its `.cu` file exists (it needs an entry in +`MAKE_ARGS`, which fails loudly with a `KeyError` if forgotten). + +## Test `__device__` helpers: `device_function_kernel` + +Helpers such as B-spline evaluation or coordinate maps are `__device__` +functions in headers. Testing them through a whole kernel is indirect. +`device_function_kernel(header_source, signature)` generates an elementwise +kernel that calls the function once per thread: + +```python +from pathlib import Path + +import numpy as np + +import cunumpy as xp +from cunumpy.testing import device_function_kernel, requires_cupy + +BSPLINES = Path("my_sim/kernels/common/bsplines.cuh").read_text() + + +@requires_cupy +def test_find_span_matches_host(): + find_span = device_function_kernel( + BSPLINES, "int find_span(const double* t, int p, double eta)" + ) + t = np.linspace(0.0, 1.0, 17) + eta = np.random.default_rng(0).random(1000) + expected = np.array([find_span_host(t, 3, e) for e in eta]) + + with xp.use_backend("cupy"): + out = xp.empty(eta.size, dtype=xp.int32) + find_span( + xp.to_cupy(t), + xp.full(eta.size, 3, dtype=xp.int32), # scalar parameter -> one per thread + xp.to_cupy(eta), + out, + eta.size, + n_threads=eta.size, + ) + np.testing.assert_array_equal(xp.to_numpy(out), expected) +``` + +In the generated kernel, pointer parameters are passed unchanged to every +thread (shared data), scalar parameters become per-thread arrays, the return +value of thread `i` goes to `out[i]`, and `n` is the number of elements. Extra +keyword arguments (`include_dirs`, `includes`, `block_size`) go to the +`CudaKernel`. + +## Test that a step stays on the device + +```python +@requires_cupy +def test_time_step_has_no_transfers(): + with xp.use_backend("cupy"): + state = make_state() + step(state, 1e-3) # warm-up: compilation, allocations + with xp.assert_no_transfers(): + step(state, 1e-3) +``` + +This catches `to_numpy()` calls, `PyccelKernel` conversions and `Kernel` +fallbacks that crept into the step. It does not see copies made outside +CuNumpy (see [Data movement](../guides/data-movement.md)). + +## Test generated headers + +When struct headers are generated with `write_cuda_header()` and committed, +add a test that regenerates them and compares (see [Kernel arguments and +structs](arguments.md), "Write the struct to a header"). A forgotten regeneration +then fails CI instead of producing a kernel that reads fields at wrong +offsets. + +## CI setup + +* Run the suite on a normal CPU runner: everything on NumPy runs, GPU cases + are reported as skipped. +* Run the same suite on a GPU runner, optionally with `CUNUMPY_CUDA_DEBUG=1` so + kernels are bounds-checked and errors are attributed to the right launch. +* Every so often, run the GPU suite under `compute-sanitizer` (see [Debugging + CUDA kernels](debugging.md)). diff --git a/docs/source/quickstart.md b/docs/source/quickstart.md index e8976c4..3cd2e99 100644 --- a/docs/source/quickstart.md +++ b/docs/source/quickstart.md @@ -1,284 +1,119 @@ -# User guide +# Quickstart -CuNumpy gives CPU and GPU programs a common NumPy-like entry point. This guide -starts with the default CPU workflow and then shows backend selection, -array-aware dispatch, data transfers, and GPU utilities. +This page is a ten-minute tour. Each section ends with a link to the guide that +covers the topic in depth. -## Install and import - -Install the package with pip: - -```bash -python -m pip install cunumpy -``` - -CuNumpy depends on NumPy and `array-api-compat`. To use NVIDIA GPUs, install a -CuPy distribution that matches your CUDA environment separately. CUDA itself -is not installed by CuNumpy. - -`array-api-compat` is a small adapter that gives NumPy and CuPy a more -consistent interface for shared array operations. You do not need to import -it directly when using CuNumpy. See [Why CuNumpy uses -`array-api-compat`](array-api-compat.md) for a beginner-friendly explanation -and examples. - -Import CuNumpy using the familiar alias `xp`: +## 1. Replace the NumPy import ```python import cunumpy as xp -a = xp.array([1.0, 2.0, 3.0]) -b = xp.asarray([4.0, 5.0, 6.0]) -print(xp.dot(a, b)) +a = xp.linspace(0.0, 1.0, 5) +b = xp.sin(a) ** 2 + xp.cos(a) ** 2 +print(b, xp.get_backend()) # [1. 1. 1. 1. 1.] numpy ``` -Most common array operations are available through the active NumPy-like -namespace, including array creation, indexing, arithmetic, reductions, -linear algebra, random operations via `xp.random`, and FFTs via `xp.fft`. -CuNumpy aims to keep the interface familiar; NumPy and CuPy are separate -libraries, so niche functions and edge-case behavior can differ. Consult the -upstream library documentation for operation-specific details. +`xp` behaves like NumPy: array creation, arithmetic, reductions, `xp.linalg`, +`xp.fft` and `xp.random` are all available. Behind the scenes every attribute +is forwarded to the `array-api-compat` module of the selected library, so the +results are ordinary NumPy (or CuPy) arrays. -## Configure the active backend +## 2. Run the same code on a GPU -The default is NumPy. Select CuPy using an environment variable before -starting Python: +Select CuPy before the arrays are created, either from the shell: ```bash -ARRAY_BACKEND=cupy python my_program.py -``` - -Or change the active backend within a running program: - -```python -xp.set_backend("cupy") -print(xp.get_backend()) # active backend name - -data = xp.arange(1000) # created on the GPU while CuPy is active +ARRAY_BACKEND=cupy python my_script.py ``` -Only `"numpy"` and `"cupy"` are valid names. If CuPy is requested but is -missing or fails its availability check, initialization falls back to NumPy. -Inspect `xp.get_backend()` after selection if the effective backend matters. -The boolean properties `xp.numpy_backend` and `xp.cupy_backend` are also -available for conditional code. - -Use a context manager when only a section should use a particular backend: +or in the program: ```python xp.set_backend("cupy") - -with xp.use_backend("numpy"): - reference = xp.linspace(0, 1, 100) - print(xp.get_backend()) # 'numpy' - -print(xp.get_backend()) # 'cupy' again +a = xp.linspace(0.0, 1.0, 5) # now a cupy.ndarray on the GPU +print(xp.get_backend()) # 'cupy', or 'numpy' if no usable GPU ``` -The old backend is restored even if the block raises an exception. Backend -selection is a process-wide setting and is not thread-safe; concurrent code -must not switch it independently from different threads or async tasks. - -## The active backend and an array's backend +If CuPy cannot be used, CuNumpy falls back to NumPy, so the same script runs on +a laptop and on a GPU node. More in [Choosing a backend](guides/backends.md). -These are deliberately separate concepts: +## 3. Write functions that follow their input -* `xp.get_backend()` returns the globally selected backend for new CuNumpy - operations. -* `xp.get_array_backend(array)` returns `"numpy"` or `"cupy"` according to - where a particular array is stored. -* `xp.get_array_module(array)` returns the matching `array_api_compat` module. - -Changing the global selection does not migrate arrays already created. For -example, the active backend can be NumPy while a CuPy array is still alive: - -```python -xp.set_backend("cupy") -gpu_values = xp.arange(4) - -xp.set_backend("numpy") -print(xp.get_backend()) # 'numpy' -print(xp.get_array_backend(gpu_values)) # 'cupy' -``` - -For a function that should follow its input array rather than global state, -use `get_array_module()`: +The active backend decides where *new* arrays are created. A reusable function +should instead follow the arrays it receives: ```python def normalize(values): - array_xp = xp.get_array_module(values) - length = array_xp.sqrt(array_xp.sum(values * values)) - return values / length -``` - -This is often the right approach for reusable functions called with arrays -created by other libraries. `is_cpu(values)` and `is_gpu(values)` are concise -boolean checks. `same_backend(a, b)` checks whether all given arrays share a -backend, and `assert_same_backend(a, b)` raises a `TypeError` with the detected -backends if they do not. - -## Convert arrays explicitly - -Use explicit conversion at CPU/GPU boundaries: - -```python -host_array = xp.to_numpy(gpu_array) -device_array = xp.to_cupy(host_array) -active_array = xp.to_cunumpy(host_array) + array_xp = xp.get_array_module(values) # numpy- or cupy-compat module + return values / array_xp.linalg.norm(values) ``` -`to_numpy()` always returns a host-side NumPy array. For a NumPy input it -uses `numpy.asarray`, so an existing array or view can be returned without a -copy. For a CuPy input it transfers data from device to host. `to_cupy()` -converts an array-like object to a CuPy array, and requires working CuPy/CUDA. -`to_cunumpy()` converts to whichever backend is currently active. - -Conversions do not mutate the source array. To avoid repeated transfer costs, -keep data on one device for a whole computational phase and transfer results -once at a boundary: - -```python -with xp.use_backend("cupy"): - signal_gpu = xp.asarray(signal_host) - spectrum_gpu = xp.fft.rfft(signal_gpu) - spectrum_host = xp.to_numpy(spectrum_gpu) -``` - -Mixing NumPy and CuPy arrays in the same operation is usually an error. Use -`assert_same_backend` to fail early in APIs that require matched arrays, or -convert the inputs to a common backend first. - -## Random generators and floating-point dtypes - -`xp.get_rng(seed)` chooses a random generator for the active backend: - -```python -rng = xp.get_rng(seed=1234) -noise = rng.normal(loc=0.0, scale=1.0, size=10_000) -``` +More in [Writing backend-agnostic code](guides/portable-code.md). -The interface is similar across NumPy and CuPy, but generated values are not -expected to be bitwise identical across backends. When reproducibility -matters, keep the backend and library versions fixed as well as the seed. +## 4. Move data explicitly -`xp.default_float_dtype()` returns the active module's `float64` dtype object. -It can be used to request an explicit portable precision: +Existing arrays never move by themselves. Convert at boundaries such as file +output or plotting: ```python -weights = xp.asarray([0.25, 0.75], dtype=xp.default_float_dtype()) +host = xp.to_numpy(a) # always a NumPy array on the host +device = xp.to_cupy(host) # a CuPy array (needs a GPU) +active = xp.to_cunumpy(host) # whatever backend is active ``` -Explicit dtypes are especially useful when code must not rely on backend- or -platform-specific inference for Python integers and floats. - -## Work with CUDA devices - -`xp.device_count()` reports the number of visible CUDA devices, or zero if -CuPy/CUDA is unavailable. It checks hardware visibility regardless of the -active array backend. +Keep the arrays on the GPU for the whole computation and transfer once. A test +can check that a time step makes no transfer at all: ```python -count = xp.device_count() -if count: - xp.set_backend("cupy") - xp.set_device(0) - values = xp.arange(100) +with xp.assert_no_transfers(): + step(state, dt) ``` -`set_device(device_id)` selects a CUDA device when the active backend is -CuPy and does nothing on NumPy. For a common one-rank-per-GPU MPI layout, -`set_device_for_rank(rank)` selects `rank % device_count()` and returns the -selected ID. Its default assumes ranks are arranged in contiguous device -blocks on each node; supply `devices_per_node` or call `set_device()` yourself -if your process mapping differs. - -`memory_info()` returns `(free_bytes, total_bytes)` for the active CUDA device, -or `None` on NumPy. This queries CUDA's runtime memory accounting, not only -CuPy's allocator. CuPy keeps freed blocks in pools for reuse, so cached memory -may remain visible as allocated. `free_memory()` releases free cached blocks -from CuPy's device and pinned-host pools; live arrays remain allocated. +More in [Moving data between host and device](guides/data-movement.md). -## Streams and synchronization +## 5. Port a compute kernel -CUDA operations are generally asynchronous with respect to the host. A stream -orders queued work; synchronization waits until that work has completed. -`xp.stream()` creates a non-blocking CuPy stream and yields it. On NumPy it is -a no-op context manager that yields `None`: +Compiled host kernels (for example with [Pyccel](https://github.com/pyccel/pyccel)) +take NumPy arrays. CuNumpy lets you pair each with a CUDA kernel and calls the +one that matches the backend: ```python -with xp.stream() as stream: - device_values = xp.to_cupy(host_values) - result = xp.exp(device_values) +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]; +} +""" -xp.synchronize() # wait before host code consumes the result -host_result = xp.to_numpy(result) -``` - -For finer control, synchronize the yielded stream with `stream.synchronize()`. -Use `xp.synchronize()` when you need to wait for all work on the current -device. Synchronization is a no-op on NumPy. -`pin_memory(host_array)` creates a pinned (page-locked) host copy. Pinned -memory can improve host/device transfer throughput, especially for workloads -that overlap transfers and computation, but should be used when transfer -profiling indicates it is useful: +def axpy(a, x, y, n): # the host version + y[:n] += a * x[:n] -```python -pinned = xp.pin_memory(host_array) -``` - -## Adapt NumPy-only or Pyccel kernels - -`xp.PyccelKernel` wraps a callable that expects NumPy arrays. It is useful for -compiled Pyccel functions and for other host-only kernels. It is an adapter, -not a compiler: compilation and imports remain the application's -responsibility. - -With the NumPy backend and NumPy arguments, the wrapped callable runs -directly. When conversion is needed (active CuPy backend or CuPy arrays among -the arguments), CuNumpy makes host copies of device arrays, invokes the -callable, copies declared in-place outputs back to the device, and converts -returned NumPy arrays to CuPy. - -```python -def axpy(alpha, x, y, out): - out[:] = alpha * x + y - return out -axpy_cpu = xp.PyccelKernel(axpy, outputs=(3,)) +kernel = xp.Kernel(axpy, xp.CudaKernel(AXPY, "axpy")) -with xp.use_backend("cupy"): - x = xp.arange(8, dtype=xp.float64) - y = xp.ones_like(x) - out = xp.empty_like(x) - result = axpy_cpu(2.0, x, y, out) +x = xp.arange(1000, dtype=xp.float64) +y = xp.zeros(1000) +kernel(2.0, x, y, 1000, n_threads=1000) # host on NumPy, CUDA on CuPy ``` -By default (`outputs=None`), all converted arrays are copied back after a -converted call. Specify only arguments the kernel may mutate to avoid copying -read-only inputs back. Positional arguments are declared by index; keyword -arguments by name. An argument passed by keyword must be declared by name, -because compiled builtins do not expose a signature for mapping it to a -position. A wrong output declaration can silently discard an in-place update -on the host copy, so include every argument the kernel writes. - -Lists, tuples, and dictionaries are traversed recursively. For instances from -your own modules, pass their module prefixes in `object_modules`, for example -`object_modules=("my_simulation.",)`. The wrapper shallow-copies those objects -and traverses their attributes, preserving the caller's original references. -`is_array` customizes which host result types are converted back to CuPy; its -default recognizes `numpy.ndarray`. - -Aliased device arrays are converted only once per call, so if the same array -appears in multiple arguments the kernel sees the same host array. Reference -cycles in supported containers are handled. `use_cupy=True` forces conversion -on and `use_cupy=False` forces it off; leave it as `None` for per-call -automatic selection. - -## Pyodide and WebAssembly - -CuNumpy's NumPy backend works in Pyodide. CuPy/CUDA and native Pyccel -compilation are not supported in that environment. Python-source functions -can still be wrapped with `PyccelKernel`; keep `use_cupy` unset or false. -Follow the [Pyodide guide](pyodide.md) for an install-and-run example and -runtime constraints. +More in [Porting kernels to the GPU](kernels/overview.md). + +## Which guide do I need? + +| I want to ... | Read | +| --- | --- | +| run array code on CPU or GPU with one code base | [Choosing a backend](guides/backends.md), [Backend-agnostic code](guides/portable-code.md) | +| understand when data is copied and avoid slow transfers | [Data movement](guides/data-movement.md) | +| control GPUs, memory pools and streams | [Devices, memory and streams](guides/gpu-devices.md) | +| run one MPI rank per GPU | [Multi-GPU programs with MPI](guides/mpi.md) | +| time GPU code or see it in Nsight | [Timing and profiling](guides/profiling.md) | +| call existing NumPy/Pyccel kernels with GPU arrays | [Host kernels with GPU data](kernels/pyccel-kernel.md) | +| write and launch CUDA kernels from Python | [Writing CUDA kernels](kernels/cuda-kernel.md) | +| port a code base kernel by kernel | [Pairing host and CUDA kernels](kernels/dispatch.md) | +| pass particle or grid data to kernels as one object | [Kernel arguments and structs](kernels/arguments.md) | +| scatter-add into a buffer another library owns | [Accumulation kernels](kernels/accumulation.md) | +| find an illegal memory access | [Debugging CUDA kernels](kernels/debugging.md) | +| test that CPU and GPU kernels agree | [Testing kernels](kernels/testing.md) | +| see everything working together | [Worked examples](examples/index.md) | diff --git a/docs/source/troubleshooting.md b/docs/source/troubleshooting.md new file mode 100644 index 0000000..245ae20 --- /dev/null +++ b/docs/source/troubleshooting.md @@ -0,0 +1,94 @@ +# Troubleshooting + +## Backend and installation + +**`set_backend("cupy")` but `get_backend()` returns `"numpy"`.** +CuPy is missing or not functional, and CuNumpy fell back to NumPy. Check +`xp.cupy_available()`, then try `python -c "import cupy; cupy.arange(3)"` to see +the real error. Typical causes: no CuPy installed, a CuPy wheel for a different +CUDA major version, no GPU visible (`CUDA_VISIBLE_DEVICES` empty, a login node +without GPUs), or a driver too old for the CUDA runtime. + +**`ARRAY_BACKEND=cupy` has no effect.** The variable is read once, when +CuNumpy is first imported. Setting `os.environ["ARRAY_BACKEND"]` after the +import does nothing; use `xp.set_backend()`. + +**Code still runs on the CPU after `set_backend("cupy")`.** Look for +`import numpy as np` used for array creation, and for `from cunumpy import +zeros`-style imports, which bind the function of the backend active at import +time. See [Writing backend-agnostic code](guides/portable-code.md). + +## Arrays and transfers + +**`TypeError` mixing NumPy and CuPy arrays.** An operation got one array of +each kind. Find where the host array came from (often `np.` instead of `xp.`, +or data loaded from disk without `to_cunumpy()`), or check inputs with +`xp.assert_same_backend()` at function entry. + +**The GPU version is slower than the CPU version.** Usually transfers in the +loop or host synchronization. Run a step inside `xp.count_transfers()` and read +the report; look for `float()`, `.item()`, `print()` or `if` on device values; +profile with `nsys` ([Timing and profiling](guides/profiling.md)). Also check +that the problem is large enough: GPUs need many thousands of elements per +operation to pay off. + +**GPU memory looks full although arrays were deleted.** CuPy's memory pool +keeps freed blocks for reuse. `xp.free_memory()` returns them to the driver. +Memory that remains in use is referenced by live arrays. + +## Kernels + +**`TypeError: argument ... must be a CuPy array ... never copied to the +device`.** A NumPy array was passed to a `CudaKernel` pointer parameter. +CUDA kernels never copy implicitly; create the array on the device or convert it +once with `xp.to_cupy()` / `xp.as_device_array()`. + +**`TypeError` about a dtype, or a non-contiguous array.** The array's dtype +must equal the C pointer type (`double*` needs `float64`, `int*` needs `int32`). +Views such as `a[:, 0]` are rejected for pointer parameters; use +`xp.ascontiguousarray()` or an `Array2D` parameter. + +**`OverflowError` for a scalar argument.** The Python integer does not fit the +declared C type, for example a size above 2^31 passed as `int`. Use `long long` +in the kernel. + +**`NotImplementedError: No CUDA version of kernel ...`.** A `Kernel` without +CUDA version was called on the CuPy backend with `missing_cuda="raise"`. Port +the kernel, or use `missing_cuda="fallback"` while porting. + +**`RuntimeWarning: No CUDA version of kernel ...: calling the host kernel`.** +The fallback is running, with host copies at every call. `catalog.summary()` +lists the kernels still to port. + +**After a `PyccelKernel` call, device arrays have old values.** The kernel +wrote an argument that is not in `outputs`. Add it, or remove `outputs` to copy +everything back. + +**`illegal memory access` somewhere unrelated.** An earlier kernel failed +asynchronously. Re-run with `CUNUMPY_CUDA_DEBUG=1` to raise at the guilty +launch, then use `compute-sanitizer` ([Debugging CUDA +kernels](kernels/debugging.md)). Restart the process afterwards; the CUDA +context is unusable. + +**A change in a `.cuh` header is ignored.** Headers included with quotes are +hashed into the compile options and trigger recompilation; headers included +with angle brackets are not. Also, a kernel object compiles once per process: +restart the process after editing. + +**The first time step is much slower.** Kernels are compiled on first call. +Compile at setup with `catalog.compile_all()` or `kernel.compile()`; later runs +load from CuPy's disk cache. + +## MPI + +**All ranks use GPU 0.** Call `xp.bind_local_device()` before `from mpi4py +import MPI`. + +**Segfault in the first MPI call with a CuPy array.** The MPI library is not +CUDA-aware. `xp.require_cuda_aware_mpi()` at start-up gives a clear message; +load a CUDA-aware MPI module or rebuild `mpi4py` against one. + +**Occasionally wrong data after an exchange.** A kernel was still writing the +send buffer. Call `xp.synchronize_for_mpi(send, recv)` before the MPI call. + +See [Multi-GPU programs with MPI](guides/mpi.md). diff --git a/pyproject.toml b/pyproject.toml index 6991fae..41f70a1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -51,7 +51,7 @@ urls."Source" = "https://github.com/max-models/cunumpy" where = [ "src" ] [tool.setuptools.package-data] -cunumpy = [ "py.typed", "*.pyi", "cuda/include/cunumpy/*.cuh" ] +cunumpy = [ "py.typed", "*.pyi", "LLM_GUIDE.md", "cuda/include/cunumpy/*.cuh" ] [tool.isort] profile = "black" diff --git a/src/cunumpy/LLM_GUIDE.md b/src/cunumpy/LLM_GUIDE.md new file mode 100644 index 0000000..f1fe3ba --- /dev/null +++ b/src/cunumpy/LLM_GUIDE.md @@ -0,0 +1,372 @@ +# cunumpy: guide for AI coding assistants + +This file tells AI coding assistants (and people in a hurry) how to write +correct code with `cunumpy`. It ships inside the installed package +(`/cunumpy/LLM_GUIDE.md`; locate it with +`python -c "import cunumpy, pathlib; print(pathlib.Path(cunumpy.__file__).parent / 'LLM_GUIDE.md')"`). +It is self-contained. The full documentation is at +https://max-models.github.io/cunumpy/ and in `docs/source/` of the repository. + +## What cunumpy is + +* A drop-in NumPy-like namespace: `import cunumpy as xp`, then `xp.zeros`, + `xp.fft.rfft`, `xp.linalg.norm`, ... Every attribute is forwarded at access + time to `array_api_compat.numpy` or `array_api_compat.cupy`, depending on the + **active backend** (`"numpy"` = CPU, `"cupy"` = NVIDIA GPU). Results are plain + `numpy.ndarray` / `cupy.ndarray`; there is no wrapper array type. +* Helpers for explicit host/device transfers, device/stream/memory control, MPI + with one rank per GPU, profiling. +* A kernel layer for porting compiled CPU kernels (Pyccel, Numba, Python loops) + to CUDA one at a time: `PyccelKernel`, `CudaKernel`, `Kernel`, + `KernelCatalog`, `KernelArguments`, `CudaStruct`, `DeviceMirror`, and test + helpers in `cunumpy.testing`. + +## Hard rules + +1. **Always `import cunumpy as xp` and call `xp.(...)` at the point of use.** + Never `from cunumpy import zeros, arange, ...`: that binds the function of the + backend active at import time. Same for `numpy_backend`/`cupy_backend`: read + `xp.cupy_backend` each time. +2. **Do not use `np.` for data that should follow the backend.** + `np.zeros` always allocates on the host. Using `numpy` for dtypes + (`np.float64`), host-only I/O, and host-side random test data is fine. +3. **Select the backend once, at the program entry point** (`ARRAY_BACKEND=cupy` + env var before import, or `xp.set_backend("cupy")` early). Library code must + not call `set_backend()`. Use `with xp.use_backend(...)` for scoped switches + (tests, CPU references). Backend state is process-global and not thread-safe. +4. **Requesting CuPy can silently fall back to NumPy** (no CuPy, no GPU, CUDA + mismatch). Check `xp.get_backend()` after `set_backend("cupy")` if it matters. +5. **Changing the backend never moves existing arrays.** Transfers are always + explicit: `xp.to_numpy(a)`, `xp.to_cupy(a)`, `xp.to_cunumpy(a)`. +6. **No transfers or host syncs inside time loops.** Avoid `to_numpy`, `float(x)`, + `x.item()`, `print(x)`, `if device_value:` in hot loops. Verify with + `xp.assert_no_transfers()` in tests. +7. **`CudaKernel` never copies.** Pointer parameters require C-contiguous CuPy + arrays of the exact dtype; NumPy arrays raise `TypeError`. Do not "fix" that + error by disabling checks; convert once with `xp.to_cupy` / + `xp.as_device_array` where the data is created. +8. **Keep `check_signature=True`** (the default). Only disable it for tiny, + already-tested kernels in hot loops. +9. **A CUDA kernel mirrors its host kernel**: same name, same positional + argument order. `Kernel` calls take positional arguments only, plus the + launch keywords `n_threads=`/`grid=`/`block=`/`shared_mem=`/`stream=`. +10. **Declare `outputs` correctly on `PyccelKernel`** (and via `host_options`). + An argument the kernel writes but that is not declared leaves stale device + data, silently. If unsure, leave `outputs=None` (copies everything back). + +## Decision guide + +| Task | Use | +| --- | --- | +| array math that should run on CPU or GPU | `xp.*` functions | +| function receiving arrays from elsewhere | `array_xp = xp.get_array_module(a)`; `xp.assert_same_backend(a, b)` | +| normalize inputs at an API boundary | `xp.to_cunumpy(a)` once | +| hand data to SciPy/matplotlib/h5py | `xp.to_numpy(a)` | +| call an existing NumPy-only kernel with GPU arrays (slow, correct) | `xp.PyccelKernel(fn, outputs=(...))` | +| launch a hand-written CUDA C kernel | `xp.CudaKernel(source, "name")` / `CudaKernel.from_file(path)` | +| host kernel + CUDA port, chosen by backend | `xp.Kernel(host_fn, cuda_kernel_or_None)` | +| many kernels in a package, ported incrementally | `xp.KernelCatalog.from_package(__name__, missing_cuda="fallback")` | +| group arrays/scalars into one kernel argument | `xp.CudaArguments` (device only), `xp.KernelArguments` (host object + device tuple), `xp.CudaStruct` (C struct) | +| CUDA struct from a Pyccel argument class | `xp.CudaStruct.from_signature(Cls.__init__, "Name")`, `xp.write_cuda_header(...)` | +| kernel writes into a host buffer owned by another library | `xp.DeviceMirror(host_array)` + `` | +| N-D indexing in CUDA, non-contiguous arrays | `Array1D`..`Array3D` 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.timed_region("name") as t:` → `t.elapsed` | +| profiler markers | `xp.nvtx_range("name")` (context manager or decorator) | +| find transfers | `with xp.count_transfers() as c: ...; print(c.report())` | +| debug an illegal memory access | `CUNUMPY_CUDA_DEBUG=1`, then `compute-sanitizer` | +| test on both backends | `cunumpy.testing.BACKENDS`, `backend` fixture, `requires_cupy` | +| test CUDA vs host kernel | `cunumpy.testing.assert_kernels_agree(kernel, make_args, n_threads=...)` | + +## API cheat sheet + +Backend and inspection: + +```python +xp.set_backend("numpy" | "cupy") # global; falls back to numpy if cupy unusable +xp.get_backend() -> "numpy" | "cupy" # active backend +with xp.use_backend("numpy"): ... # temporary, exception-safe +xp.numpy_backend, xp.cupy_backend # bools for the active backend +xp.cupy_available() -> bool # CuPy importable and functional (cached) +xp.get_array_backend(a) -> "numpy" | "cupy" # where this array lives +xp.get_array_module(a) # array_api_compat module for a's backend +xp.is_cpu(a), xp.is_gpu(a) +xp.same_backend(*arrays) -> bool +xp.assert_same_backend(*arrays) # TypeError if mixed +xp.__version__ +``` + +Transfers: + +```python +xp.to_numpy(a) # -> numpy.ndarray; copies only if a is CuPy +xp.to_cupy(a) # -> cupy.ndarray; copies only if a is not CuPy; ImportError w/o CuPy +xp.to_cunumpy(a) # -> array of the active backend +xp.as_device_array(value, dtype=None, ndim=None, *, name=None) + # contiguous CuPy array of dtype: returned as is; else one copy; + # RuntimeError on the NumPy backend +with xp.count_transfers() as c: ... # c.total, c.to_host, c.to_device, + # c.kernel_conversions, c.fallbacks, c.events, c.report() +with xp.assert_no_transfers(): ... # AssertionError with report if anything copied +``` + +Only transfers through cunumpy are counted (not raw `cupy.asarray`, `.get()`, +`float(device_scalar)`, or `DeviceMirror.to_host()/to_device()`). + +Random numbers and dtypes: + +```python +rng = xp.get_rng(seed=None) # numpy or cupy Generator for the active backend +xp.default_float_dtype() # float64 of the active backend +``` + +NumPy and CuPy generators give different sequences for the same seed. For +identical data on both backends: `xp.to_cunumpy(np.random.default_rng(s).random(n))`. + +Devices, memory, streams (all safe on NumPy: no-ops / neutral values): + +```python +xp.device_count() -> int # visible GPUs, 0 without CuPy +xp.set_device(i) +xp.memory_info() -> (free, total) | None +xp.free_memory() # release CuPy pool cached blocks +xp.synchronize() +with xp.stream() as s: ... # s is None on NumPy; prefer xp.synchronize() +xp.pin_memory(host_array) # pinned copy; needs CuPy +``` + +MPI: + +```python +xp.set_backend("cupy") +xp.bind_local_device() # BEFORE `from mpi4py import MPI`; uses local_rank() +from mpi4py import MPI +xp.require_cuda_aware_mpi() # collective; RuntimeError if MPI can't take device buffers +xp.mpi_is_cuda_aware(comm=None) -> bool +xp.local_rank() -> int # node-local rank from launcher env vars +xp.synchronize_for_mpi(*buffers) # before every MPI call that touches device buffers +``` + +Profiling: + +```python +with xp.timed_region("name", sync=True) as t: ... # t.name, t.elapsed (s), t.synced +with xp.nvtx_range("name", color=None): ... # also usable as @decorator +``` + +`PyccelKernel`: + +```python +k = xp.PyccelKernel(fn, use_cupy=None, object_modules=(), is_array=None, outputs=None) +k(*args, **kwargs) +``` + +NumPy path: calls `fn` directly. CuPy path: copies CuPy arrays (also inside +lists/tuples/dicts and objects whose class module starts with an +`object_modules` prefix) to the host, calls `fn`, copies `outputs` (positional +indices or keyword names; default: all) back, converts returned NumPy arrays to +CuPy. Does not compile anything. + +`CudaKernel`: + +```python +k = xp.CudaKernel(source, name, *, block_size=128, options=(), include_dirs=(), + source_dir=None, structs=(), template_args=None, + check_signature=True, debug=None) +k = xp.CudaKernel.from_file("push/push_cuda.cu") # name "push" +ks = xp.CudaKernel.all_from_file("ops.cu") # dict name -> kernel +k(*args, n_threads=None, grid=None, block=None, shared_mem=0, stream=None) +k.compile(); k.is_compiled; k.launch_shape(n_threads) -> (grid, block) +k.included_headers; k.compile_options(); k.debug_active() +xp.CudaKernelVariants(factory).get(*key); .compile_all(keys, jobs=1) +xp.ctype_of(np.float64) == "double" +xp.cuda_kernel_names(source); xp.parse_cuda_signature(source, name) +xp.cuda_include_dir() +``` + +* Source must declare `extern "C" __global__` (templates: plain `__global__` + plus `template_args=(np.float64, 3)`). +* Type mapping (LP64): `int`=int32, `long`/`long long`=int64, `float`=float32, + `double`=float64, `complex`=complex128, `bool`=bool; `int64_t`, + `size_t` etc. supported. +* Scalars: Python `int`/`float`/`bool` are cast with range checks + (`OverflowError`); NumPy scalars must match or cast safely. +* Pointer params: C-contiguous CuPy arrays, exact dtype (`void*` any). + `ArrayND` params: CuPy arrays of dtype T and ndim N, any strides. +* Compiled lazily with NVRTC on first call; cached on disk by CuPy; quoted + `#include "..."` headers are hashed into the options so edits recompile. +* Creating a `CudaKernel` does not import CuPy; compiling needs a GPU. + +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..Array3D: data, shape[], strides[] (elements), a(i, j), size() +#include // cunumpy_atomic_add(double*|float*, v), _2d(data, n1, i, j, v), _3d(...) +``` + +`-DCUNUMPY_BOUNDS_CHECK` (on in debug mode) bounds-checks array view indexing. + +`Kernel` and `KernelCatalog`: + +```python +k = xp.Kernel(host_kernel, cuda_kernel=None, *, name=None, + missing_cuda="raise" | "fallback", cuda_path=None, host_options=None) +k(*args, n_threads=..., grid=None, block=None, shared_mem=0, stream=None) +k.get_kernel(); k.compile(); k.has_cuda + +catalog = xp.KernelCatalog.from_package(__name__, *, host_suffix="_kernels", + cuda_suffix="_cuda.cu", missing_cuda="raise", host_options=None, + include_dirs=None, **cuda_options) +catalog["push"]; catalog.summary(); catalog.with_cuda; catalog.without_cuda +catalog.compile_all(jobs=1); catalog.parity_cases(); catalog.register(kernel) +``` + +Package layout for `from_package`: `pkg//_kernels.py` defines +function `` (host); optional `pkg//_cuda.cu` defines +`__global__ void (...)`. Kernel folders must be importable packages. +`host_options` may be a function of the kernel name, e.g. +`lambda name: {"outputs": OUTPUTS[name]}`. + +Argument objects: + +```python +class Dev(xp.CudaArguments): # flattened into several CUDA params + def __init__(self, x, n): super().__init__(x, n) + +class Args(xp.KernelArguments): # one object, host form + device form + def __host_args__(self): return host_object # host kernel gets this + def __cuda_args__(self): return (arr, n, ...) # CUDA kernel gets these, flattened + +S = xp.CudaStruct("S", [("x", "double*"), ("n", "long long"), ("a", "Array2D")]) +S.declaration; S.dtype; S.to_header(path); value = S(x=..., n=..., a=...) +S = xp.CudaStruct.from_signature(Cls.__init__, "S", int_type="long long") +xp.write_cuda_header("args.cuh", [S1, S2]) +xp.CudaKernel(S.declaration + src, "k", structs=[S]) +xp.resolve_host_args(args, kwargs) +``` + +Only top-level arguments are resolved. Cache both forms lazily and invalidate +them when the underlying arrays are replaced. A packed struct holds device +addresses: re-pack after replacing an array. + +`DeviceMirror`: + +```python +m = xp.DeviceMirror(host_numpy_array) # TypeError if not numpy.ndarray +m.device # CuPy copy (lazy) on CuPy; the host array itself on NumPy +m.zero(); m.to_device(); m.to_host() # to_host copies in place; no-ops on NumPy +m.rebind(new_host_array) # after the owner reallocates +``` + +Debugging: + +```python +xp.set_cuda_debug(True); xp.get_cuda_debug(); with xp.cuda_debug(): ... +xp.CudaKernel(..., debug=True) # env: CUNUMPY_CUDA_DEBUG=1 +``` + +Debug mode adds `-lineinfo -DCUNUMPY_BOUNDS_CHECK` at compile time and +synchronizes after each launch (errors become `RuntimeError` naming the kernel). +Enable it before kernels compile. After an illegal memory access, the process +must be restarted. + +Testing (`import cunumpy.testing`; not imported by `import cunumpy`): + +```python +from cunumpy.testing import BACKENDS, backend, requires_cupy, \ + assert_kernels_agree, device_function_kernel + +@pytest.mark.parametrize("backend", BACKENDS) # "numpy" always, "cupy" if GPU +def test_x(backend): + with xp.use_backend(backend): ... + +assert_kernels_agree(kernel, make_args, *, n_threads=None, grid=None, block=None, + rtol=1e-12, atol=0.0, n_calls=1, outputs=None, seed=0) +# make_args(backend, seed) -> tuple of positional args, built with the backend active + +k = device_function_kernel(header_source, "int f(const double* t, int p, double x)") +k(t, p_array, x_array, out, n, n_threads=n) # scalars become per-thread arrays +``` + +## Canonical patterns + +Backend-agnostic function: + +```python +def normalize(v): + array_xp = xp.get_array_module(v) + return v / array_xp.linalg.norm(v) +``` + +Script entry point: + +```python +import cunumpy as xp + +def main(use_gpu: bool): + xp.set_backend("cupy" if use_gpu else "numpy") + print("backend:", xp.get_backend()) + data = xp.to_cunumpy(load_host_data()) # one transfer in + for _ in range(n_steps): + data = update(data) # no transfers here + save(xp.to_numpy(data)) # one transfer out +``` + +Kernel pair: + +```python +SRC = r""" +#include +extern "C" __global__ +void scale(double* x, double a, long long n) { + CUNUMPY_THREAD_1D(i, n); + x[i] *= a; +} +""" + +def scale_host(x: "float[:]", a: float, n: int): + for i in range(n): + x[i] *= a + +scale = xp.Kernel(scale_host, xp.CudaKernel(SRC, "scale"), + host_options={"outputs": (0,)}) +scale(x, 2.0, x.size, n_threads=x.size) +``` + +Parity test: + +```python +def make_args(backend, seed): + x = xp.to_cunumpy(np.random.default_rng(seed).random(1000)) + return (x, 2.0, x.size) + +def test_scale(): + assert_kernels_agree(scale, make_args, n_threads=1000) +``` + +## Common mistakes to avoid + +* Writing `if xp.get_backend() == "cupy": cupy.foo(a) else: numpy.foo(a)` for + functions that `xp.foo` already dispatches. Branch only when the libraries + genuinely differ (SciPy vs `cupyx.scipy`), and branch on the *array* + (`xp.is_gpu(a)`), not the global setting. +* Calling `xp.to_cupy()` inside a loop or per kernel call. +* Passing NumPy arrays or non-contiguous views (`a[:, 0]`) to `CudaKernel` + pointer parameters. Use `xp.ascontiguousarray()` or an `Array2D` param. +* `int` loop indices in CUDA for arrays that may exceed 2^31 elements; use + `long long` (the index macros do). +* Forgetting `if (i >= n) return;` (or `CUNUMPY_THREAD_1D`) in a kernel. +* Plain `+=` from many threads into one cell; use `cunumpy_atomic_add`. +* Timing GPU code with `time.perf_counter()` without synchronizing; use + `xp.timed_region()`. +* Importing `mpi4py.MPI` before `xp.bind_local_device()`; skipping + `xp.synchronize_for_mpi()` before MPI calls on device buffers. +* Comparing CPU and GPU results with exact equality where the GPU uses atomics + or a different reduction order; use a tolerance. +* Assuming `xp.random.seed(s)` gives the same numbers on both backends. +* Expecting `DeviceMirror` copies or raw CuPy conversions to show up in + `count_transfers()`. +* Writing code that requires CuPy at import time. `import cunumpy` and creating + `CudaKernel` objects work without CuPy; import `cupy` lazily, only on the GPU + path, if you need it at all.