diff --git a/CHANGELOG.md b/CHANGELOG.md index 4b30ca5..1b26e2b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [Unreleased] +### Fixed +- Preserve Fortran order when copying CuPy arrays to the host through + `to_numpy`, `host_call`, `evaluate_on_host`, and `PyccelKernel` conversions. + ### Added - `kernel_testing.emulated_launches()` runs every `CudaKernel` launch in a block on the CPU, on the host buffers of the fake CuPy arrays, so code that launches diff --git a/docs/source/api.md b/docs/source/api.md index 195a2a1..38dab05 100644 --- a/docs/source/api.md +++ b/docs/source/api.md @@ -206,6 +206,10 @@ Converts to a host-side NumPy array. CuPy arrays are copied from device to host. Other array-like inputs are passed through `numpy.asarray`; NumPy arrays may therefore be returned as-is rather than copied. +F-contiguous CuPy inputs keep F order on the host; other CuPy inputs become +C-contiguous. Arbitrary device strides are not preserved. See +[Array ordering and strides](guides/array-ordering.md). + ### `to_cupy(array)` Converts an array-like input to a CuPy array. Raises `ImportError` if CuPy or diff --git a/docs/source/guides/array-ordering.md b/docs/source/guides/array-ordering.md new file mode 100644 index 0000000..be572f1 --- /dev/null +++ b/docs/source/guides/array-ordering.md @@ -0,0 +1,139 @@ +# Array ordering and strides + +CuNumpy supports C-ordered (row-major) and F-ordered (column-major) arrays on +both NumPy and CuPy. Choose the layout when allocating an array, and keep it +consistent with the compiled kernels that consume it. + +For a two-dimensional array `a[row, column]`, C order stores each row together; +F order stores each column together. Indexing and shape stay the same. F order +can help code that processes columns of a particle buffer, but measure the +whole workload, including kernels and copies, before choosing a layout. + +## Allocate and inspect + +```python +import cunumpy as xp + +markers = xp.zeros((100, 6), dtype=xp.float64, order="F") +assert markers.flags.f_contiguous + +positions = markers[:, 1:4] +assert positions.flags.f_contiguous +assert positions.shape == (100, 3) +``` + +A block of complete columns of an F-contiguous array remains F-contiguous. +A block of complete rows of a C-contiguous array remains C-contiguous. +Stepped slices such as `markers[::2, :]` generally have neither layout. +Check `a.flags.c_contiguous`, `a.flags.f_contiguous`, and `a.strides` when +debugging a kernel boundary. Strides are measured in bytes. Arrays with +singleton dimensions or no elements can satisfy both contiguity flags. + +Use `xp.asfortranarray(a)` to require F order or `xp.ascontiguousarray(a)` to +require C order. These may allocate a copy; a transpose changes the shape and +indexing meaning and is not a substitute for an order conversion. + +## Host and device conversions + +| Operation | Layout behavior | +| --- | --- | +| `xp.to_numpy(cupy_array)` | F-contiguous input becomes F-contiguous host storage; other input becomes C-contiguous host storage. | +| `xp.to_numpy(numpy_array)` | Uses `numpy.asarray`; an ordinary NumPy array keeps its layout and strides without copying. | +| `xp.to_cupy(numpy_array)` | Uses `cupy.asarray`, preserving C or F order for contiguous input. | +| `xp.to_cunumpy(a)` | Uses the conversion for the active backend. | +| `xp.host_call`, `xp.evaluate_on_host` | Device arguments follow `to_numpy`; returned NumPy arrays follow `to_cupy`. | +| `xp.kernels.PyccelKernel` | Device arguments use the same host conversion as `to_numpy`; declared outputs are copied back into the original device arrays. | + +The device-to-host conversion uses `array.get(order="A")`. This preserves +F-contiguous layout, including complete column blocks, but does not preserve +arbitrary strides or shared storage between distinct views. Calling CuPy's +`array.get()` directly defaults to C order; see the +[CuPy array documentation](https://docs.cupy.dev/en/stable/reference/generated/cupy.ndarray.html). + +The host function itself controls the layout of any new arrays it creates. +Returning a C-ordered result from `host_call` does not make that result F-ordered +just because an input was F-ordered. + +```{note} +Preserving F order in CuPy-to-host conversions is an unreleased fix. Earlier +releases used `array.get()` and produced C-ordered host copies. Applications +depending on F order at a host-kernel boundary need a release containing this fix. +``` + +## Copies, masks, and assignment + +`a.copy()` defaults to C order on both NumPy and CuPy. Use `a.copy(order="K")` +to retain the layout of a C- or F-contiguous array, or `a.copy(order="F")` to +require F order. `order="K"` is not a promise to preserve arbitrary strides. + +Boolean row masks and fancy indexing create new arrays. Do not rely on those +results retaining F order. For example, a boolean row selection from an +F-ordered marker buffer produces a C-ordered result on CuPy. Normalize an +independent result before passing it to a kernel requiring F order: + +```python +alive = markers[:, 0] >= 0 +selected = xp.asfortranarray(markers[alive]) +assert selected.flags.f_contiguous +``` + +Assignment into an existing buffer keeps the destination's layout: + +```python +markers[...] = markers.copy() # C-ordered temporary, F-ordered destination +assert markers.flags.f_contiguous +``` + +Likewise, `xp.algorithms.compact_by_mask(alive, markers)` gathers selected rows +and writes them into the front of the original buffer, preserving that +buffer's layout. Only the returned count of leading rows is valid. However, +`markers[:n_kept]` generally loses F contiguity when it trims rows: the columns +still have the original buffer's stride. Pass the full buffer plus a valid-row +count to a kernel designed for that interface, or make an F-contiguous copy +of the prefix for a kernel that requires one. A stride-aware CUDA array view +can consume the prefix directly. + +Check or normalize arrays created by concatenation, file readers, and other +libraries before passing them across a layout-sensitive boundary. + +## Compiled kernel boundaries + +Pyccel array signatures specify an expected order. A two-dimensional +`"float[:, :]"` argument expects C order; use `"float[:, :](order=F)"` for an +F-ordered argument. Every caller must supply the declared layout, including +temporary arrays and array attributes in argument objects. +`PyccelKernel` preserves F-contiguous device arguments when moving them to +the host; it does not inspect signatures or fix mismatches for you. The same +requirement applies to a `Kernel` using its host fallback without a CUDA version. +See [Host kernels with GPU data](../kernels/pyccel-kernel.md). + +CUDA `Array2D` and other `ArrayND` view arguments carry shape and strides, +so indexing through the view supports C order, F order, and strided slices. +The same applies to array view fields in `CudaStruct` and to emulated launches. +Raw pointer arguments and `CArrayND` views require C-contiguous storage. +Hand-written indexing such as `data[row * n_columns + column]` assumes C order; +use the array view's indexing to handle both layouts. See +[Index macros and array views](../kernels/cuda-kernel.md). + +Some convenience helpers deliberately require C order: + +| Helper | Contract | +| --- | --- | +| `xp.as_device_array` | Returns a C-contiguous device array, copying if needed. | +| `xp.kernels.as_kernel_array` | Returns a C-contiguous array on the side of `like`. | +| `xp.kernels.kernel_output` | Yields a C-contiguous working buffer and copies updates back into the original output when needed. | + +For an F-order kernel, allocate or normalize its arguments explicitly instead +of using these helpers to prepare F-ordered inputs. + +## Check values and layout separately + +Numerical equality does not check memory order. Tests for layout-sensitive +code should assert the host argument's contiguity inside the called function, +as well as checking values. Use a two-dimensional shape with both dimensions +greater than one so C and F order are distinguishable. + +The fake CuPy supports these conversion checks without a GPU. Its +`kernel_testing.host_buffer(a)` exposes the existing NumPy storage without a +copy, preserving layout and strides. It does not exercise the `to_numpy` +transfer path; test that path separately. See [Testing kernels](../kernels/testing.md). diff --git a/docs/source/guides/data-movement.md b/docs/source/guides/data-movement.md index 72d5980..21c7603 100644 --- a/docs/source/guides/data-movement.md +++ b/docs/source/guides/data-movement.md @@ -21,6 +21,10 @@ implicitly: every transfer is a visible function call. `to_numpy()`, on CuPy like `to_cupy()`. * None of them modify the source array or change the active backend. +CuPy-to-host conversions preserve F-contiguous layout; other device inputs +become C-contiguous host arrays. See [Array ordering and strides](array-ordering.md) +for conversion rules, copies, and compiled kernel requirements. + ## Transfer at boundaries, not in loops Load or generate data, move it to the device once, run the whole computation diff --git a/docs/source/index.md b/docs/source/index.md index aed8c8f..49f1fd2 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -47,6 +47,7 @@ quickstart guides/backends guides/portable-code guides/data-movement +guides/array-ordering guides/gpu-devices guides/execution-helpers guides/mpi diff --git a/docs/source/kernels/pyccel-kernel.md b/docs/source/kernels/pyccel-kernel.md index 5be0abd..05315d6 100644 --- a/docs/source/kernels/pyccel-kernel.md +++ b/docs/source/kernels/pyccel-kernel.md @@ -37,6 +37,11 @@ with xp.use_backend("cupy"): smooth_kernel(field, out) # out is updated on the device ``` +Device-to-host conversion preserves F-contiguous arrays. The host kernel must +accept that order: the wrapper does not adapt arrays to a compiled signature. +See [Array ordering and strides](../guides/array-ordering.md) for Pyccel signatures +and the operations that can change layout. + ## Declare the outputs The wrapper cannot know which arguments the function writes. By default it diff --git a/pyproject.toml b/pyproject.toml index 8346503..0b053c7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ requires = [ "setuptools", "wheel" ] [project] name = "cunumpy" -version = "0.6.1" +version = "0.6.2" description = "Simple wrapper for numpy and cupy. Replace `import numpy as np` with `import cunumpy as xp`." readme = "README.md" keywords = [ "python" ] diff --git a/src/cunumpy/xp.py b/src/cunumpy/xp.py index d667280..bb94eb5 100644 --- a/src/cunumpy/xp.py +++ b/src/cunumpy/xp.py @@ -269,7 +269,7 @@ def synchronize() -> None: def _to_numpy(array: Any) -> np.ndarray: """`to_numpy` without transfer counting, for internal use.""" if get_array_backend(array) == "cupy": - return array.get() + return array.get(order="A") return np.asarray(array) @@ -289,6 +289,7 @@ def to_numpy(array: Any) -> np.ndarray: A CuPy array is copied to the host, which `count_transfers()` counts as a ``to_host`` transfer; anything else is passed through `numpy.asarray`. + Fortran-contiguous CuPy arrays keep F order; other CuPy arrays use C order. """ result = _to_numpy(array) if _COUNTERS and get_array_backend(array) == "cupy": diff --git a/tests/unit/test_cupy.py b/tests/unit/test_cupy.py index a8e804d..f03baaf 100644 --- a/tests/unit/test_cupy.py +++ b/tests/unit/test_cupy.py @@ -2,6 +2,50 @@ import pytest import cunumpy as xp +from cunumpy.kernels import PyccelKernel + + +@pytest.mark.parametrize("order", ["C", "F"]) +@pytest.mark.parametrize("column_block", [False, True]) +@pytest.mark.parametrize( + "conversion", ["to_numpy", "host_call", "evaluate_on_host", "pyccel_kernel"] +) +def test_host_conversion_preserves_order(order, column_block, conversion): + if not xp.cupy_available(): + pytest.skip("CuPy not installed or not functional") + + expected = np.array(np.arange(30.0).reshape(5, 6), order=order) + with xp.use_backend("cupy"): + device = xp.to_cupy(expected) + if column_block: + device = device[:, 1:4] + expected = expected[:, 1:4] + + def check_host(array): + assert isinstance(array, np.ndarray) + np.testing.assert_array_equal(array, expected) + assert array.flags.f_contiguous == (order == "F") + assert array.flags.c_contiguous == (order == "C") + return array + + if conversion == "to_numpy": + check_host(xp.to_numpy(device)) + else: + if conversion == "host_call": + result = xp.host_call(check_host, device) + elif conversion == "evaluate_on_host": + + class Model: + @xp.evaluate_on_host + def evaluate(self, array): + return check_host(array) + + result = Model().evaluate(device) + else: + result = PyccelKernel(check_host)(device) + + assert xp.is_gpu(result) + check_host(xp.to_numpy(result)) def test_to_cupy_available(): diff --git a/tests/unit/test_execution_transfers.py b/tests/unit/test_execution_transfers.py index da5ba0a..aa9806e 100644 --- a/tests/unit/test_execution_transfers.py +++ b/tests/unit/test_execution_transfers.py @@ -31,10 +31,10 @@ def set(self, host): cp.copies.append("set") np.copyto(self.array, host) - def get(self, out=None, stream=None): + def get(self, out=None, stream=None, order="C"): cp.copies.append(("get", stream)) if out is None: - return self.array.copy() + return self.array.copy(order=order) np.copyto(out, self.array) return out diff --git a/tests/unit/test_transfers.py b/tests/unit/test_transfers.py index 90201bf..f72493e 100644 --- a/tests/unit/test_transfers.py +++ b/tests/unit/test_transfers.py @@ -42,8 +42,8 @@ def shape(self): def dtype(self): return self.data.dtype - def get(self): - return self.data.copy() + def get(self, order="C"): + return self.data.copy(order=order) def __setitem__(self, key, value): self.data[key] = value.data if isinstance(value, _FakeDeviceArray) else value