diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 5f88caa..0bf30df 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -1,19 +1,55 @@ +name: Continuous integration + on: push: - branches: [ main ] - pull_request: + branches: [main, f/improvements] + workflow_dispatch: -name: Continuous integration +permissions: + contents: read jobs: ci: - runs-on: ${{ matrix.runner }} - strategy: - matrix: - runner: [a40, orin] + name: A40 CUDA and Python tests + runs-on: [self-hosted, linux, x64, a40] + timeout-minutes: 40 + + env: + JOBS: "4" + SANITIZER_REPEATS: "3" + COMPUTE_SANITIZER: /usr/local/cuda/bin/compute-sanitizer + steps: - - name: checkout code - uses: actions/checkout@v4 - - name: run test script - run: bash ./ci/script.sh + - name: Check out code + uses: actions/checkout@v7 + with: + persist-credentials: false + + - name: Check CUDA + run: | + nvidia-smi + /usr/local/cuda/bin/nvcc --version + /usr/local/cuda/bin/compute-sanitizer --version + + - name: Prepare Python environment + run: | + python3 -m venv build/test-venv + build/test-venv/bin/python -m pip install --upgrade pip + build/test-venv/bin/python -m pip install ./python pytest + + - name: Build, run unit tests, and check CUDA memory + run: | + bash scripts/check_memory.sh -DSM_ARCH=86 -DCPPVERSION=20 + + - name: Run Python tests with pytest + run: bash scripts/run_python_tests.sh + - name: Save diagnostic logs + if: always() + uses: actions/upload-artifact@v7 + with: + name: a40-test-logs + path: | + build/test/Testing/Temporary/ + build/sanitizer/ + if-no-files-found: ignore diff --git a/.gitignore b/.gitignore index 032fffb..0fbce20 100644 --- a/.gitignore +++ b/.gitignore @@ -269,4 +269,23 @@ __pycache__/ */build/ -cmake-* \ No newline at end of file +cmake-* + +# CMake/CTest artefacts (including accidental in-source builds) +CMakeCache.txt +CMakeFiles/ +CTestTestfile.cmake +Testing/ +Makefile +cmake_install.cmake +install_manifest.txt +**/_deps/ + +# Build products from this project +/gputils_main +/lib/ +/test/gputils_test* + +# Temporary files produced by the C++ file-I/O tests +/myTest.bt +/myTest.dtensor diff --git a/CHANGELOG.md b/CHANGELOG.md index 5d853fe..da196a8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,68 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). + +## v2.0.0 - 04-11-2026 + +### Added + +- `scripts/run_python_tests.sh` runs the Python pytest suite with configurable interpreter/build paths and optional C++ exchange tests. +- File-I/O tests cover shapes, precision, binary layout, non-contiguous inputs, malformed files, filesystem errors, and actual C++/Python exchanges using temporary fixtures. +- `LeastSquaresBatchedWorkspace` preallocates reusable status storage and one-entry matrix-pointer arrays only for batches of one. +- `Svd::rankRelative(rtol, atol = 0)` adds GPU rank calculation using per-matrix thresholds `max(atol, rtol * sigmaMax)` and validated tolerances. +- Build/test and CUDA memory-check scripts run C++/Python tests, repeated leak checks, and save sanitizer logs. +- Size-safe `Session::recordAllocation(size_t)`, `recordDeallocation(size_t)`, and counted `cudaRelease` methods. +- `DTensor::getRowsInto(destination, rowsFrom, rowsTo, matIdx)` validates and copies rows into preallocated storage on the source stream. +- `DTensor::leastSquaresBatched(b, infoArray)` supports reusable device status storage without per-call allocation. +- `Session::numStreams()` exposes the initialized stream count. +- Configurable GPU benchmarks compare batched GEMM and row-extraction methods, with correctness checks, host/GPU timings, and CSV results. +- Regression tests cover validation, files, GEMM coefficients, stream ordering, ownership, rank tolerances, and reusable buffers. + +### Changed + +- **Breaking:** replace `incrementAllocatedBytes(int)` with separate size-safe allocation/deallocation methods. +- **Breaking:** copy assignment deep-copies; move assignment transfers storage and empties the source. Use slices for aliases and `deviceCopyTo` for preallocated copies. +- Reshape retains and reuses owned matrix-pointer capacity where possible. +- Column-major uploads use the supplied host vector directly; only row-major uploads allocate a conversion buffer. +- `setStreamIdx` returns a reference for lvalues and moves rvalues, avoiding implicit allocation and deep copying. +- `getRows` uses one asynchronous 2D copy and inherits the source stream; `getRowsInto` preserves the destination's stream index. +- Tensor constructors validate dimensions, allocation overflow, vector length, and storage mode before allocating. +- Validate reshape sizes, slices, random bounds, row extraction, debug access, empty factoriser inputs, and QR output shapes. +- Reject zero streams and reconfiguration after initialization; validate stream indices against the initialized count. +- Document numerical preconditions, buffer lifetimes, stream dependencies, reusable workspaces, and binary-format assumptions. +- Ignore generated build files, dependencies, binaries, and temporary test artifacts, including in-source builds. + +### Fixed + +- Nullspace projection uses constructor-allocated scratch and stream-ordered copy-back to avoid GEMM aliasing; tests cover dense batches, repeated calls, streams, and allocation reuse. +- Givens annihilation uses the identity rotation for an already-zero pair, avoiding NaNs and preserving both rows; covered for float and double. +- Python file I/O rejects incomplete headers, payload-size mismatches, zero dimensions, empty/scalar arrays, and object/zero-size dtypes before loading or overwriting data. +- Python binary loading preserves column-major element order for single matrices, with a rectangular round-trip regression test. +- Batched Cholesky and least-squares support single matrices using constructor-allocated pointer storage without allocations during reusable solves. +- Size QR workspace for the largest requirement across factorisation, least-squares, and Q extraction; add large square/tall regressions for both precisions. +- Share stream configuration across translation units using an inline static `Session` member, with a two-file regression test. +- Prevent SVD constructor exception leaks using a `unique_ptr` owner and raw observer; verify cleanup with allocation-failure tests. +- Tensor constructors release GPU storage on subsequent host-allocation failures; fault-injection tests verify cleanup for float, double, and int. +- Session shutdown completes pending work, releases library handles, and destroys owned CUDA streams. +- Single and batched Cholesky queries, factorisations, solves, and status tensors use the matrix-selected stream. +- Host transfers, scalar reads, initialization, and pointer uploads complete on their selected stream; transpose results inherit the source stream. +- Device copies wait for prior destination work and complete on the source stream before returning. +- QR debugging-helper uploads complete on the selected solver stream before returning, fixing intermittent stream-ordering failures. +- Track owned bytes and pointer capacity independently of shape, count successful allocations/frees, and preserve accounting across reshape, views, and moves. +- Reset `Svd::rank()` counters on each call; SVD factorisation, rank kernels, and results use the input stream. +- Assignment releases previous storage and safely handles metadata, self-assignment, and empty tensors; deep copies complete on the source stream. +- Enable testing at the project root so `ctest --test-dir build` discovers the test suite. +- Make the Python-format binary-reader test independent of generated Python files and the working directory. +- Reject malformed text and incomplete or incorrectly sized binary files; report write failures and handle short filenames and empty data safely. +- Avoid null-pointer arithmetic when slicing a single-matrix tensor along the matrix axis. +- QR operations use the matrix-selected cuSolver/cuBLAS stream; `getQR` synchronizes before reading R on the host. +- Givens coefficient kernels use the selected cuBLAS stream, check launch errors, and restore the previous pointer mode. +- Make test-stream setup independent of execution order and add nonblocking-stream regressions for QR, Givens, and row extraction. +- Correct the row-slicing test to select a column before slicing its rows. + + @@ -53,8 +115,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added - The project now uses the GNU General Public License (GPL) v3; license file added -- Introduces new Python package for storing and loading numpy arrays; in can be installed with `pip install gputils-api`; - unit tests and documentation +- Python package for storing and loading NumPy arrays, with tests and documentation; install with `pip install gputils-api`. ### Fixed @@ -70,8 +131,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -- Method `saveToFile` becomes `DTensor::saveToFile(std::string pathToFile, Serialisation ser)`, i.e., the user can - choose whether to save the file as a text (ASCII) file, or a binary one +- `DTensor::saveToFile(std::string pathToFile, Serialisation ser)` supports text (ASCII) and binary output. ### Changed @@ -115,12 +175,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed -- (Breaking change) The methods `CholeskyFactoriser::factorise`, `CholeskyFactoriser::solve`, `QRFactoriser::factorise`, - `QRFactoriser::leastSquares`, and `QRFactoriser::getQR` are now `void` and do not - return a status code. Instead, a status code is returned by called `statusCode`. - This change leads to a reduction in data being downloaded from the GPU. -- In `Svd` a status code (`bool`) is returned from `Svd::factorise` only if the - `#GPUTILS_DEBUG_MODE` is defined, otherwise, the method returns always `true`. +- **Breaking:** Cholesky factorisation/solve and QR factorisation/least-squares/extraction return `void`; query status separately to reduce GPU downloads. +- `Svd::factorise` returns a status only with `GPUTILS_DEBUG_MODE`; otherwise it always returns `true`. - New base class `IStatus` used for a universal implementation of `info()` @@ -142,7 +198,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed -- Memory management improvements: we got rid of `pointerToMatrices`, which would unnecessarily allocate memory and `addAB` does not allocate any new memory internally. +- Remove unnecessary `pointerToMatrices` allocations; `addAB` no longer allocates memory internally. y (device memory to device memory) ``` -The copy constructor has also been implemented; to hard-copy a `DTensor` just -do `DTensor myCopy(existingTensor)`. +Copy construction and copy assignment create independent owning storage: + +```c++ +DTensor myCopy(existingTensor); +myCopy = existingTensor; // deep copy, replacing myCopy's allocation +myCopy = std::move(existingTensor); // transfer storage; existingTensor becomes empty +``` + +**Compatibility change:** copy assignment previously aliased the source. It now +matches the copy constructor's deep-copy semantics. Copies preserve the source's +shape and stream index, and complete on the source's stream before returning. +Self-copy and self-move assignment are safe no-ops; empty and moved-from tensors +can be assigned again. Moving an owning tensor transfers its allocation without +copying device data. Replacing an owning destination releases its old allocation +and invalidates views referring to it. + +For intentional aliasing, use the slice constructor explicitly. Slices remain +non-owning: their backing allocation must outlive them. Copying a slice produces +independent owning storage; moving a slice preserves its non-owning nature and +does not extend the backing allocation's lifetime. In particular, do not move a +slice into the tensor that owns its backing allocation; copy it instead. +Assigning into a slice replaces that view rather than writing through it. + +For repeated operations requiring no allocation, preallocate the destination +and use `existingTensor.deviceCopyTo(myCopy)` instead of copy assignment. +Ensure previous asynchronous uses of a destination have completed before +replacing its storage, and establish dependencies for source writes from other +streams. Move assignment can free old storage; it is not a guarantee of bounded +real-time latency. Lastly, a not so efficient method that should only be used for debugging, if at all, is the `()` operator (e.g., `x(i, j, k)`), which fetches @@ -239,6 +284,48 @@ Then, we can solve the system by A.leastSquaresBatched(B); ``` +For repeated solves, construct a workspace once during setup. It supports any +batch size, including one, and subsequent calls allocate no GPUtils memory: + +```c++ +LeastSquaresBatchedWorkspace workspace(A); +A.leastSquaresBatched(B, workspace); +Session::getInstance().synchronizeStream(A.streamIdx()); +std::vector status; +workspace.info().download(status); +``` + +The workspace owns one status integer per matrix and allocates two one-entry +device pointer arrays only for a batch of one. Reuse it with tensors of the same +shape; their storage addresses may differ. Keep it alive until the solve stream +completes and use separate workspaces for concurrent streams. + +For batches with existing device pointer arrays, the status-buffer overload +also remains available: + +```c++ +DTensor info(A.numMats()); // setup: one integer per matrix +A.leastSquaresBatched(B, info); // no GPUtils workspace allocation in this call +Session::getInstance().synchronizeStream(A.streamIdx()); +std::vector status; +info.download(status); +``` + +The original one-argument overload remains available and constructs its own +temporary workspace. The status-only overload rejects single matrices that lack +device pointer arrays; use an explicit workspace for allocation-free calls in +that case. The supplied status buffer may be larger than the batch; +only its first `A.numMats()` entries are written. Both overloads overwrite A +and B; restore the original inputs before solving a new problem. Keep the +inputs and supplied buffer alive until A's stream completes. Reuse on that +stream is ordered; use separate buffers or synchronization for concurrent +streams. The buffer's own stream index does not select the solve stream. + +`CholeskyBatchFactoriser` also supports a batch of one. Its constructor +allocates two one-entry device pointer arrays only for that case; factorisation +and repeated solves allocate no GPUtils memory. Ordinary tensors retain their +existing allocation behaviour, and larger batches reuse tensor pointer arrays. + The `DTensor` `B` will be overwritten with the solution. > [!IMPORTANT] @@ -264,6 +351,11 @@ If the file extension is `.bt` (binary tensor), the data will be stored in binar The structure of the binary encoding is similar to that of the text encoding: the first three `uint64_t`-sized positions correspond to the number of rows, columns and matrices, followed by the elements of the tensor. +The dimensions and payload are stored in the host's native byte order and the +payload uses the native representation of `T`; the file does not record its +element type. The Python API uses little-endian dimensions, so interoperability +requires a little-endian host and a matching NumPy dtype. Readers reject missing, +malformed, truncated, or extra data. To load a tensor from a file, the static function `DTensor::parseFromFile(filename)` can be used. For example: @@ -352,6 +444,12 @@ the given matrix $B$. If you want the factoriser to keep your matrix, you need to set the third argument of the above constructor to `false`. +With `destroyMatrix=false`, SVD owns an independent input copy through a +`std::unique_ptr`; that copy is released automatically on normal destruction or +if later constructor setup throws. With `destroyMatrix=true`, SVD only borrows +the caller's tensor and never frees it. The caller must keep borrowed input +alive and retain its shape/stream binding while the solver uses it. + After you have factorised the matrix, you can access $S$, $V'$ and, perhaps, $U$. You can do: @@ -371,6 +469,35 @@ auto U = svdEngine.leftSingularVectors(); if (U) std::cout << "U = " << U.value(); ``` +### Numerical rank tolerances + +After a successful factorisation, `svdEngine.rank(epsilon)` retains its existing +absolute-cutoff behavior. To use a scale-aware cutoff explicitly: + +```c++ +const auto& ranks = svdEngine.rankRelative(1e-6f); // atol defaults to zero +const auto& noisyRanks = svdEngine.rankRelative(1e-6f, 1e-8f); +``` + +`rankRelative(rtol, atol)` counts singular values **strictly greater than** +`max(atol, rtol * sigmaMax)`. Each matrix in a batch uses its own largest +singular value. Both tolerances must be finite and nonnegative. Zero matrices +have rank zero, including when both tolerances are zero. Computation and counter +reset stay on the selected stream without downloading singular values or +allocating new workspace. + +With `atol = 0`, rescaling a matrix rescales its threshold, preserving numerical +rank within floating-point limits. A positive `atol` intentionally imposes an +absolute cutoff, useful when data have a known noise scale. `rtol` is required; +one conventional choice is `max(rows, cols) * std::numeric_limits::epsilon()`, +but applications with noisy data may need a larger tolerance. + +Both rank methods return the same reusable internal device buffer. Download or +deep-copy a result before making another rank call if you need a snapshot, and +order consumers on other streams explicitly. For example, the second call above +overwrites the result referenced by `ranks`. `Nullspace` continues to use its +existing absolute rank policy; this new method does not change its behavior. + ## 4. Projection onto a nullspace The nullspace of a matrix is computed by SVD. @@ -390,11 +517,23 @@ nullspace projection operators `NN' = (N1N1', ..., NkNk')`. This allows the user to project-in-place onto the nullspace. ```c++ -DTensor vectors(m, 1, k); +DTensor vectors(n, 1, k); N.project(vectors); std::cout << vectors << "\n"; ``` +`Nullspace` allocates an `n × 1 × k` scratch tensor in its constructor. Projection +computes into that separate storage and copies back into the existing vector +allocation, avoiding overlapping GEMM input/output and per-call allocations. +The constructor completes projector setup before returning. Each `project` call +queues GEMM and copy-back on the vector tensor's selected stream and returns +without synchronizing it; the vector's data pointer and stream index stay unchanged. +Keep the vectors and `Nullspace` object alive until that stream completes. +Calls on one stream may be queued consecutively. Before reusing the same object +on another stream, complete its prior work or establish a stream dependency; +concurrent projections need separate `Nullspace` objects. Inputs must contain +one `n`-element vector per original matrix; incompatible shapes are rejected. + ## 5. Other We can get the total allocated bytes (on the GPU) with @@ -403,12 +542,30 @@ size_t allocatedBytes = Session::getInstance().totalAllocatedBytes(); ``` +This counts live allocations managed by GPUtils, including retained matrix-pointer +capacity after reshaping, rather than only the bytes implied by the current shape. +Shrinking preserves pointer capacity; growing within that capacity reuses it. +Slices do not count borrowed element storage a second time. A reshaped slice may +own a separate pointer array, which is counted and released independently. +Moves and swaps transfer allocation sizes along with ownership. + +For manual accounting, `recordAllocation(size_t)` and `recordDeallocation(size_t)` +replace `incrementAllocatedBytes(int)`. They check overflow/underflow without +narrowing byte counts. Prefer `Session::cudaAllocate` and +`Session::cudaRelease(pointer, allocatedBytes)` so accounting follows successful +CUDA calls automatically. Do not also manually account for those wrapper calls. + GPUtils supports multiple streams. By default a single stream is created, but you can set the number of streams you need with ```c++ /* This needs to be the first line in your code */ Session::setStreams(4); // create 4 strems ``` +Call `setStreams` with a positive count before the first `Session` or tensor +allocation initializes the session. Changing the count afterward is rejected. +The configured count is shared across translation units through C++17 inline +static storage: configuration in one source file also applies when another source +file first initializes the singleton. Then, you can use `setStreamIdx` to select a stream to go with your instance of `DTensor` ```c++ auto a = DTensor::createRandomTensor(3, 6, 4, -1, 1).setStreamIdx(0); @@ -417,4 +574,63 @@ auto b = DTensor::createRandomTensor(3, 6, 4, -1, 1).setStreamIdx(1); Session::getInstance().synchronizeAllStreams(); ``` +Stream ordering rules: + +- Single and batched Cholesky, QR, and SVD execute on their bound matrix's selected + stream. Select that stream before constructing the solver. Keep solver workspaces + and operands alive until completion; operands produced on another stream require + an explicit dependency before the solve. Solver `info()` tensors follow the + matrix stream, as do SVD outputs and transpose results. +- Existing `upload`, `download`, debug element access, and deep copies wait for + completion on the tensor/source stream. They do not depend on legacy default-stream + synchronization. The row-major upload conversion buffer stays alive until transfer + completion. No new public asynchronous host-transfer API is introduced. +- `source.deviceCopyTo(destination)` waits for prior work on the destination's + selected stream, then completes the copy on the source stream before returning. + Neither stream index changes, so subsequent destination-stream consumers are safe. + This is a blocking convenience operation, not a transfer-overlap API. +- Arithmetic involving several tensors does not automatically connect all operand + streams. Callers must synchronize or establish CUDA-event dependencies between + different producer and consumer streams. A selected stream records an execution + convention, not a record of every external access to a tensor's raw pointer. + +## 6. Build, test, and check CUDA memory + +The scripts in `scripts/` can be launched from any working directory. They require +CMake, a supported CUDA toolkit/GPU, a C++ compiler, and Python with NumPy. For an +isolated Python environment (from the repository root): + +```bash +python3 -m venv build/test-venv +build/test-venv/bin/python -m pip install numpy +``` + +Run all C++ and Python unit tests: + +```bash +./scripts/build_and_test.sh -DSM_ARCH=120 +``` + +Use `SM_ARCH=86` for A40 or `SM_ARCH=87 -DCPPVERSION=17` for Orin instead of `120`. +Additional arguments are forwarded to CMake; the default build type is +`RelWithDebInfo`. `GPUTILS_BUILD_DIR`, `JOBS`, and `PYTHON` override the build +directory, build parallelism (default 4), and Python interpreter. Without `PYTHON`, +the script uses `test-venv/bin/python` in the build directory if available, then +`python3`. Python fixtures stay in the build tree. + +Build/test and run the complete CUDA suite under Compute Sanitizer memcheck: + +```bash +./scripts/check_memory.sh -DSM_ARCH=120 +``` + +This checks invalid device memory accesses and device/heap leaks, runs the full +suite three times in shuffled order within one process, and checks `gputils_main` +separately. It enables CUDA line information in the default build configuration, +preserves failure exit codes, and requires clean error/leak summaries. Complete +logs and generated fixtures are stored under `build/sanitizer/`. Set +`COMPUTE_SANITIZER` to override its executable path and `SANITIZER_REPEATS` to +change the repeat count. Compute Sanitizer checks CUDA memory, not arbitrary +C++ host-memory leaks. Do not use sanitizer runs for performance measurements. + ## Happy number crunching! diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt new file mode 100644 index 0000000..790b4b5 --- /dev/null +++ b/benchmark/CMakeLists.txt @@ -0,0 +1,13 @@ +cmake_minimum_required(VERSION 3.20) +project(GPUtilsBenchmarks LANGUAGES CXX CUDA) + +find_package(CUDAToolkit REQUIRED) +if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + set(CMAKE_BUILD_TYPE Release CACHE STRING "Build type" FORCE) +endif() +add_executable(gputils_benchmark benchmark.cu) +target_compile_features(gputils_benchmark PRIVATE cxx_std_17 cuda_std_17) +target_include_directories(gputils_benchmark PRIVATE ../include) +target_link_libraries(gputils_benchmark PRIVATE CUDA::cublas CUDA::cusolver) +target_compile_definitions(gputils_benchmark PRIVATE + GPUTILS_BENCH_CONFIG="$") diff --git a/benchmark/README.md b/benchmark/README.md new file mode 100644 index 0000000..70a373a --- /dev/null +++ b/benchmark/README.md @@ -0,0 +1,73 @@ +# GPU benchmarks + +Build this standalone project from the repository root. It does not modify the +library's algorithms or require GoogleTest. Choose an architecture supported by +your GPU and CUDA toolkit (the current remote Blackwell GPU uses `120`). + +```bash +cmake -S benchmark -B benchmark/build \ + -DCMAKE_CUDA_COMPILER=/usr/local/cuda/bin/nvcc \ + -DCMAKE_BUILD_TYPE=Release -DCMAKE_CUDA_ARCHITECTURES=120 +cmake --build benchmark/build -j4 +benchmark/build/gputils_benchmark > benchmark/build/results.csv \ + 2> benchmark/build/metadata.txt +``` + +Defaults: 30 warm-up iterations and 200 measured samples per variant, for both +`float` and `double`. The suite covers four GEMM and four extraction shapes. +Use representative application dimensions when available: + +```bash +benchmark/build/gputils_benchmark --samples 1000 --warmup 100 \ + --gemm 64 32 48 128 --rows 512 256 64 > benchmark/build/custom.csv +``` + +`--gemm M N K BATCH` selects `(M x K) * (K x N)` and replaces the default GEMM +cases; `--rows ROWS COLS SELECTED` replaces the default extraction cases. Both +families always run. Extraction uses the middle contiguous range of rows from a +single column-major matrix. All options require positive integers. + +## Comparisons + +- GEMM `preallocated`: current `DTensor::addAB` versus strided batched cuBLAS, + with existing tensors. This isolates execution; both share the already-built + tensor pointer arrays. For a batch of one the current implementation uses + ordinary GEMM. +- GEMM `allocate`: construction of A, B and C, device-to-device initialization + of A and B, GEMM, completion, and destruction. The pointer variant uses actual + DTensors; the strided variant uses raw device buffers without matrix-pointer + arrays. This estimates the benefit of a future pointer-free storage path; + changing `addAB` alone will not remove the current DTensor pointer allocations. +- Rows `preallocated`: original per-row cuBLAS copy loop versus one asynchronous + 2D device copy versus one extraction kernel, using an existing destination. +- Rows `allocate`: actual `getRows` is reported separately, alongside the old + per-row loop and alternatives allocating an identically shaped DTensor. + Completion and destruction are included. Source upload is not timed. + +All candidates are checked against expected values before timing. GEMM checks +nonzero inputs and non-default alpha/beta, then times alpha=1.25, beta=0 so the +output does not accumulate across repetitions. Row checks compare every output +element to the source range. Check failures cause a nonzero exit. + +## Interpreting output + +CSV dimensions are M,N,K for GEMM and ROWS,COLS,SELECTED for extraction. Times +are in microseconds; percentiles use nearest ranks. Host timings measure one +completed operation using `steady_clock`, including the final stream wait. +Allocation-inclusive timings also include allocation, initialization and frees. +CUDA-event timings are collected in a separate pass for preallocated cases; +they include device work and any gaps between host submissions. Event timings +are not emitted for allocation cases, where host work is central to the result. + +This measures isolated, synchronized latency, not pipelined throughput. The +stream is the session's cuBLAS stream. Library defaults for math mode are kept. +GPU name, compute capability, driver/runtime versions, build type and sample +counts are written to stderr. Run on an otherwise idle GPU, retain the metadata, +record the source revision, compiler version, GPU clocks/power/temperature and +persistence state, and repeat the run to check stability. Persistence avoids +driver teardown; it does not lock clocks or eliminate thermal and scheduling +variation. p99 and maximum are observations, not real-time bounds. CUDA event +resolution and synchronization overhead can dominate very short operations. + +The default shapes are exploratory; use application workloads before deciding +which production implementation to adopt. No performance threshold is enforced. diff --git a/benchmark/benchmark.cu b/benchmark/benchmark.cu new file mode 100644 index 0000000..0539622 --- /dev/null +++ b/benchmark/benchmark.cu @@ -0,0 +1,238 @@ +#include "tensor.cuh" +#include +#include +#include +#include + +// Standalone comparisons: no production implementations are replaced here. +struct Buffer { + void *data = nullptr; + explicit Buffer(size_t bytes) { gpuErrChk(cudaMalloc(&data, bytes)); } + ~Buffer() { gpuErrChk(cudaFree(data)); } + Buffer(const Buffer &) = delete; + Buffer &operator=(const Buffer &) = delete; + template T *as() const { return static_cast(data); } +}; + +struct Events { + cudaEvent_t start, stop; + Events() { gpuErrChk(cudaEventCreate(&start)); gpuErrChk(cudaEventCreate(&stop)); } + ~Events() { gpuErrChk(cudaEventDestroy(start)); gpuErrChk(cudaEventDestroy(stop)); } +}; + +struct GemmShape { int m, n, k, batch; }; +struct RowShape { int rows, cols, selected; }; +struct Options { + int samples = 200, warmup = 30; + std::vector gemms{{8,8,8,128}, {32,32,32,128}, {64,32,48,64}, {128,128,128,16}}; + std::vector rows{{32,32,1}, {128,128,8}, {512,256,64}, {1024,512,256}}; +}; + +template +void stridedGemm(T *a, T *b, T *c, GemmShape s, T beta = 0) { + const T alpha = T(1.25); + auto handle = Session::getInstance().cuBlasHandle(); + if constexpr (std::is_same_v) { + gpuErrChk(cublasSgemmStridedBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, + s.m, s.n, s.k, &alpha, a, s.m, 1LL*s.m*s.k, b, s.k, 1LL*s.k*s.n, + &beta, c, s.m, 1LL*s.m*s.n, s.batch)); + } else { + gpuErrChk(cublasDgemmStridedBatched(handle, CUBLAS_OP_N, CUBLAS_OP_N, + s.m, s.n, s.k, &alpha, a, s.m, 1LL*s.m*s.k, b, s.k, 1LL*s.k*s.n, + &beta, c, s.m, 1LL*s.m*s.n, s.batch)); + } +} + +template +__global__ void extractRows(const T *src, T *dst, int rows, int cols, int first, int count) { + size_t idx = blockIdx.x * blockDim.x + threadIdx.x; + if (idx < size_t(count)*cols) dst[idx] = src[first + idx%count + (idx/count)*rows]; +} + +template +void copyRows(const T *src, T *dst, RowShape s, int first, int method, cudaStream_t stream) { + if (method == 0) { + for (int row = 0; row < s.selected; ++row) { + auto handle = Session::getInstance().cuBlasHandle(); + if constexpr (std::is_same_v) { + gpuErrChk(cublasScopy(handle, s.cols, src+first+row, s.rows, dst+row, s.selected)); + } else { + gpuErrChk(cublasDcopy(handle, s.cols, src+first+row, s.rows, dst+row, s.selected)); + } + } + } else if (method == 1) { + gpuErrChk(cudaMemcpy2DAsync(dst, s.selected*sizeof(T), src+first, s.rows*sizeof(T), + s.selected*sizeof(T), s.cols, cudaMemcpyDeviceToDevice, stream)); + } else { + extractRows<<>>( + src, dst, s.rows, s.cols, first, s.selected); + gpuErrChk(cudaGetLastError()); + } +} + +double percentile(std::vector values, double fraction) { + std::sort(values.begin(), values.end()); + return values[static_cast(std::ceil(fraction*values.size()))-1]; +} + +void measure(const Options &options, const std::string &key, cudaStream_t stream, + bool deviceTiming, const std::function &operation) { + for (int i = 0; i < options.warmup; ++i) { operation(); gpuErrChk(cudaStreamSynchronize(stream)); } + std::vector host, device; + Events events; + for (int i = 0; i < options.samples; ++i) { + // No other work is in flight in this standalone process. + auto start = std::chrono::steady_clock::now(); + operation(); + gpuErrChk(cudaStreamSynchronize(stream)); + auto stop = std::chrono::steady_clock::now(); + host.push_back(std::chrono::duration(stop-start).count()); + } + // Separate pass: event instrumentation must not inflate host latency. + if (deviceTiming) { + for (int i = 0; i < options.samples; ++i) { + gpuErrChk(cudaEventRecord(events.start, stream)); + operation(); + gpuErrChk(cudaEventRecord(events.stop, stream)); + gpuErrChk(cudaEventSynchronize(events.stop)); + float ms; + gpuErrChk(cudaEventElapsedTime(&ms, events.start, events.stop)); + device.push_back(ms*1000.0); + } + } + auto print = [&](const char *clock, const std::vector &values) { + std::cout << key << ',' << clock << ',' << options.samples << ',' + << percentile(values, .5) << ',' << percentile(values, .95) << ',' + << percentile(values, .99) << ',' << *std::max_element(values.begin(), values.end()) << '\n'; + }; + print("host", host); + if (deviceTiming) print("cuda_event", device); +} + +template +void checkValues(const T *data, const std::vector &expected) { + std::vector actual(expected.size()); + gpuErrChk(cudaMemcpy(actual.data(), data, actual.size()*sizeof(T), cudaMemcpyDeviceToHost)); + for (size_t i = 0; i < actual.size(); ++i) { + if (!std::isfinite(actual[i]) || std::abs(actual[i]-expected[i]) > T(1e-5)*std::max(T(1), std::abs(expected[i]))) + throw std::runtime_error("benchmark correctness check failed at element " + std::to_string(i)); + } +} + +template +void benchmarkGemm(const Options &options, GemmShape s, cudaStream_t stream) { + std::vector a(size_t(s.m)*s.k*s.batch, T(1)), b(size_t(s.k)*s.n*s.batch, T(2)); + std::vector c(size_t(s.m)*s.n*s.batch, T(3)); + DTensor A(a,s.m,s.k,s.batch), B(b,s.k,s.n,s.batch), C(c,s.m,s.n,s.batch); + C.addAB(A,B,T(1.25),T(.5)); + gpuErrChk(cudaStreamSynchronize(stream)); + checkValues(C.raw(), std::vector(c.size(), T(2.5*s.k+1.5))); + C.upload(c); + stridedGemm(A.raw(),B.raw(),C.raw(),s,T(.5)); + gpuErrChk(cudaStreamSynchronize(stream)); + checkValues(C.raw(), std::vector(c.size(), T(2.5*s.k+1.5))); + std::string key = std::string("gemm,") + (sizeof(T)==4?"float":"double") + ',' + + std::to_string(s.m)+','+std::to_string(s.n)+','+std::to_string(s.k)+','+std::to_string(s.batch); + measure(options,key+",preallocated,pointer",stream,true,[&] { C.addAB(A,B,T(1.25),T(0)); }); + measure(options,key+",preallocated,strided",stream,true,[&] { stridedGemm(A.raw(),B.raw(),C.raw(),s); }); + // Includes all three allocations, pointer-array initialization, D2D input + // initialization, GEMM completion, and destruction. No host uploads here. + measure(options,key+",allocate,pointer",stream,false,[&] { + DTensor x(s.m,s.k,s.batch), y(s.k,s.n,s.batch), z(s.m,s.n,s.batch); + gpuErrChk(cudaMemcpyAsync(x.raw(),A.raw(),a.size()*sizeof(T),cudaMemcpyDeviceToDevice,stream)); + gpuErrChk(cudaMemcpyAsync(y.raw(),B.raw(),b.size()*sizeof(T),cudaMemcpyDeviceToDevice,stream)); + z.addAB(x,y,T(1.25),T(0)); + gpuErrChk(cudaStreamSynchronize(stream)); + }); + measure(options,key+",allocate,strided",stream,false,[&] { + Buffer x(a.size()*sizeof(T)), y(b.size()*sizeof(T)), z(c.size()*sizeof(T)); + gpuErrChk(cudaMemcpyAsync(x.data,A.raw(),a.size()*sizeof(T),cudaMemcpyDeviceToDevice,stream)); + gpuErrChk(cudaMemcpyAsync(y.data,B.raw(),b.size()*sizeof(T),cudaMemcpyDeviceToDevice,stream)); + stridedGemm(x.as(),y.as(),z.as(),s); + gpuErrChk(cudaStreamSynchronize(stream)); + }); +} + +template +void benchmarkRows(const Options &options, RowShape s, cudaStream_t stream) { + std::vector source(size_t(s.rows)*s.cols); + for (size_t i=0; i A(source,s.rows,s.cols), result(s.selected,s.cols); + int first = (s.rows-s.selected)/2; + std::vector expected(size_t(s.selected)*s.cols); + for (int c=0;c output(s.selected,s.cols); + copyRows(A.raw(),output.raw(),s,first,method,stream); + gpuErrChk(cudaStreamSynchronize(stream)); + }); + } +} + +int main(int argc, char **argv) { + try { + Options options; + auto integer = [&](int &i) { + if (++i>=argc) throw std::invalid_argument("missing option value"); + std::string value(argv[i]); + size_t used=0; + long n=std::stol(value,&used); + if (used!=value.size() || n<=0 || n>INT_MAX) throw std::invalid_argument("expected positive integer"); + return static_cast(n); + }; + for (int i=1;irows) throw std::invalid_argument("selected rows exceed source rows"); + options.rows={{rows,cols,selected}}; + } else throw std::invalid_argument("options: --samples N --warmup N --gemm M N K BATCH --rows ROWS COLS SELECTED"); + } + if (std::string(GPUTILS_BENCH_CONFIG)!="Release" && std::string(GPUTILS_BENCH_CONFIG)!="RelWithDebInfo") + throw std::runtime_error("use a Release or RelWithDebInfo build"); + cudaDeviceProp props{}; + int device=0, driver=0, runtime=0; + gpuErrChk(cudaGetDevice(&device)); gpuErrChk(cudaGetDeviceProperties(&props,device)); + gpuErrChk(cudaDriverGetVersion(&driver)); gpuErrChk(cudaRuntimeGetVersion(&runtime)); + std::cerr << "GPU: " << props.name << "; compute capability: " << props.major << '.' << props.minor + << "; driver API: " << driver << "; runtime: " << runtime << "; build: " << GPUTILS_BENCH_CONFIG + << "; warmup: " << options.warmup << "; samples: " << options.samples << '\n'; + cudaStream_t stream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(),&stream)); + std::cout << "operation,dtype,dim1,dim2,dim3,batch,scope,variant,clock,samples,median_us,p95_us,p99_us,max_us\n"; + std::cout << std::fixed << std::setprecision(3); + for (auto shape:options.gemms) { + checkedTensorSize(shape.m,shape.k,shape.batch); + checkedTensorSize(shape.k,shape.n,shape.batch); + checkedTensorSize(shape.m,shape.n,shape.batch); + benchmarkGemm(options,shape,stream); benchmarkGemm(options,shape,stream); + } + for (auto shape:options.rows) { + checkedTensorSize(shape.rows,shape.cols,1); + benchmarkRows(options,shape,stream); benchmarkRows(options,shape,stream); + } + gpuErrChk(cudaDeviceSynchronize()); + std::cerr << "All benchmark correctness checks passed.\n"; + } catch (const std::exception &error) { std::cerr << error.what() << '\n'; return 1; } +} diff --git a/ci/script.sh b/ci/script.sh index ca13ddf..322bcc4 100644 --- a/ci/script.sh +++ b/ci/script.sh @@ -45,6 +45,11 @@ tests() { # -- run tests ctest --test-dir ./build/test --output-on-failure + # The exchange tests require the C++ helper, so run them after building it. + PYTHONPATH="$PWD/python" \ + GPUTILS_FILE_EXCHANGE="$PWD/build/test/gputils_file_exchange" \ + "$PWD/python/venv/bin/python" "$PWD/python/test/test.py" CppFileExchangeTest -v + if [ -z "${hwInfoOrin}" ]; then # -- run compute sanitizer diff --git a/include/tensor.cuh b/include/tensor.cuh index b73e280..f3b0da6 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -15,6 +15,12 @@ #include #include #include +#include +#include +#include +#include +#include +#include /** @@ -116,11 +122,24 @@ inline void gpuAssert(T code, const char *file, int line, bool abort = true) { /* ================================================================================================ * SESSION * ================================================================================================ */ -/** - * Total number of allocated streams - * Can be changed with Session::setStreams() - */ -static size_t s_numStreams = 1; +namespace gputils_detail { +// Independent of CUDA allocation, allowing large-count arithmetic to be tested +// without requiring multi-gigabyte GPU allocations. +class AllocationCounter { + size_t m_bytes = 0; +public: + size_t bytes() const { return m_bytes; } + void recordAllocation(size_t bytes) { + if (bytes > std::numeric_limits::max() - m_bytes) + throw std::overflow_error("allocation counter overflow"); + m_bytes += bytes; + } + void recordDeallocation(size_t bytes) { + if (bytes > m_bytes) throw std::logic_error("allocation counter underflow"); + m_bytes -= bytes; + } +}; +} /** * Singleton for Cuda library handles. @@ -131,12 +150,17 @@ static size_t s_numStreams = 1; * The cuSolver handle can be accessed anywhere by `Session::getInstance().cuSolverHandle()` */ class Session { + // C++17 inline storage: one configured count across all translation units. + inline static size_t s_numStreams = 1; + static bool &initialized() { static bool value = false; return value; } public: /** * Sets the total number of available streams * @param numStreams number of streams (default: 1) */ static void setStreams(size_t numStreams) { + if (numStreams == 0) throw std::invalid_argument("stream count must be positive"); + if (initialized()) throw std::logic_error("streams must be configured before Session initialization"); s_numStreams = numStreams; } @@ -163,24 +187,30 @@ private: gpuErrChk(cusolverDnCreate(&m_cusolverHandles[i])); gpuErrChk(cusolverDnSetStream(m_cusolverHandles[i], m_cublasStreams[i])); } + initialized() = true; } ~Session() { + // Finish all owned streams before releasing handles or stream resources. + synchronizeAllStreams(); for (size_t i = 0; i < m_numCublasHandlesStreams; i++) { gpuErrChk(cublasDestroy(m_cublasHandles[i])); gpuErrChk(cusolverDnDestroy(m_cusolverHandles[i])); + gpuErrChk(cudaStreamDestroy(m_cublasStreams[i])); } } std::vector m_cublasHandles; std::vector m_cublasStreams; std::vector m_cusolverHandles; - size_t m_bytesAllocated = 0; + gputils_detail::AllocationCounter m_allocations; size_t m_numCublasHandlesStreams = 1; public: Session(Session const &) = delete; + size_t numStreams() const { return m_numCublasHandlesStreams; } + void operator=(Session const &) = delete; /** @@ -207,9 +237,10 @@ public: * @return CUDA error */ cudaError_t cudaAllocate(void **d, size_t s) { + if (s > std::numeric_limits::max() - totalAllocatedBytes()) return cudaErrorInvalidValue; cudaError_t err = cudaMalloc(d, s); if (err == cudaSuccess) { - m_bytesAllocated += s; + recordAllocation(s); } return err; } @@ -217,13 +248,22 @@ public: /** * @return Total allocated bytes */ - size_t totalAllocatedBytes() const { return m_bytesAllocated; } + size_t totalAllocatedBytes() const { return m_allocations.bytes(); } /** - * Increment counter of allocated bytes - * @param s allocated bytes (can be negative) + * Record externally managed allocations/deallocations without narrowing to + * int or unsigned negative arithmetic. Prefer cudaAllocate/cudaRelease for + * device storage; their counters update only after successful CUDA calls. */ - void incrementAllocatedBytes(int s) { m_bytesAllocated += s; } + void recordAllocation(size_t bytes) { m_allocations.recordAllocation(bytes); } + void recordDeallocation(size_t bytes) { m_allocations.recordDeallocation(bytes); } + + cudaError_t cudaRelease(void *data, size_t bytes) { + if ((!data && bytes != 0) || bytes > totalAllocatedBytes()) return cudaErrorInvalidValue; + cudaError_t err = cudaFree(data); + if (err == cudaSuccess) recordDeallocation(bytes); + return err; + } /** * Synchronize stream @@ -260,6 +300,24 @@ enum StorageMode { defaultMajor = columnMajor }; +// Explicitly sized tensors are nonempty. Default construction remains available +// for an uninitialized tensor. CUDA BLAS and solver APIs use int dimensions. +template +size_t checkedTensorSize(size_t rows, size_t cols, size_t mats) { + if (!rows || !cols || !mats || rows > INT_MAX || cols > INT_MAX || mats > INT_MAX) + throw std::invalid_argument("tensor dimensions must be positive and fit in int"); + const size_t limit = std::numeric_limits::max(); + if (rows > limit / cols || rows * cols > limit / mats || + rows * cols * mats > limit / sizeof(T) || mats > limit / sizeof(T *)) + throw std::invalid_argument("tensor dimensions overflow allocation size"); + return rows * cols * mats; +} + +inline void validateStorageMode(StorageMode mode) { + if (mode != columnMajor && mode != rowMajor) + throw std::invalid_argument("invalid storage mode"); +} + /** * This library uses tensors to store and manipulate data on a GPU device. @@ -268,9 +326,12 @@ enum StorageMode { * Tensors can be used to do a batched operation on many similar-sized matrices or vectors in parallel. * @tparam T type of data stored in tensor */ +template class LeastSquaresBatchedWorkspace; + TEMPLATE_WITH_TYPE_T class DTensor { private: + void leastSquaresBatchedImpl(DTensor &b, DTensor &infoArray, T **matrixPointers, T **rhsPointers); T *m_d_data = nullptr; ///< Pointer to device data T **m_d_ptrMatrices = nullptr; ///< Pointer to matrices in tensor size_t m_numRows = 0; ///< Number of rows @@ -278,19 +339,35 @@ private: size_t m_numMats = 0; ///< Number of matrices bool m_doDestroyData = false; ///< Whether to destroy memory bool m_doDestroyPtrMatrices = false; ///< Whether to destroy memory + size_t m_dataBytes = 0; ///< Actual size of owned element allocation + size_t m_matrixPointerCapacity = 0; ///< Capacity of owned matrix-pointer allocation size_t m_idxStream = 0; ///< Stream index (defaults to 0) + cudaStream_t selectedStream() const { + cudaStream_t stream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(m_idxStream), &stream)); + return stream; + } + + T downloadElement(size_t offset) const { + T value; + auto stream = selectedStream(); + gpuErrChk(cudaMemcpyAsync(&value, m_d_data + offset, sizeof(T), cudaMemcpyDeviceToHost, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); + return value; + } + void destroy() { if (m_doDestroyData) { - if (m_d_data) gpuErrChk(cudaFree(m_d_data)); - m_d_data = nullptr; - Session::getInstance().incrementAllocatedBytes(-m_numRows * m_numCols * m_numMats * sizeof(T)); + gpuErrChk(Session::getInstance().cudaRelease(m_d_data, m_dataBytes)); } if (m_doDestroyPtrMatrices) { - if (m_d_ptrMatrices) gpuErrChk(cudaFree(m_d_ptrMatrices)); - m_d_ptrMatrices = nullptr; - Session::getInstance().incrementAllocatedBytes(-m_numMats * sizeof(T *)); + gpuErrChk(Session::getInstance().cudaRelease(m_d_ptrMatrices, m_matrixPointerCapacity * sizeof(T *))); } + m_d_data = nullptr; + m_d_ptrMatrices = nullptr; + m_doDestroyData = m_doDestroyPtrMatrices = false; + m_dataBytes = m_matrixPointerCapacity = 0; } /** @@ -333,9 +410,11 @@ private: public: /** - * Set the stream ID + * Set the stream ID without copying device storage. Named tensors return + * a reference; temporary tensors transfer ownership into the returned value. */ - DTensor setStreamIdx(size_t); + DTensor &setStreamIdx(size_t) &; + DTensor setStreamIdx(size_t) &&; size_t streamIdx() const { return m_idxStream; } @@ -408,7 +487,22 @@ public: * Move constructor. * @param other tensor to move to newly constructed tensor */ - DTensor(DTensor &&other); + DTensor(DTensor &&other) noexcept; + + /** Exchange storage, ownership, shape and stream metadata without copying. */ + void swap(DTensor &other) noexcept { + using std::swap; + swap(m_d_data, other.m_d_data); + swap(m_d_ptrMatrices, other.m_d_ptrMatrices); + swap(m_numRows, other.m_numRows); + swap(m_numCols, other.m_numCols); + swap(m_numMats, other.m_numMats); + swap(m_doDestroyData, other.m_doDestroyData); + swap(m_doDestroyPtrMatrices, other.m_doDestroyPtrMatrices); + swap(m_dataBytes, other.m_dataBytes); + swap(m_matrixPointerCapacity, other.m_matrixPointerCapacity); + swap(m_idxStream, other.m_idxStream); + } /** * Slice constructor. @@ -465,12 +559,17 @@ public: /** * Device-to-device copy. + * Completes on the source stream before returning and waits for prior work + * on the destination's selected stream. Neither tensor's stream ID changes. * @param other target tensor */ void deviceCopyTo(DTensor &other) const; /** * Slices rows from specified matrix. + * The copy is enqueued on this tensor's stream; the result inherits that + * stream. Keep the source storage alive until the copy completes, and + * synchronize before consuming the result on another stream. * @param rowsFrom index to slice rows from (zero-indexed) * @param rowsTo index to slice rows to (inclusive) * @param matIdx index of matrix to slice rows from (zero-indexed) @@ -478,6 +577,16 @@ public: */ DTensor getRows(size_t rowsFrom, size_t rowsTo, size_t matIdx) const; + /** + * Copy rows into existing storage without allocating or resizing. + * destination must have shape (rowsTo-rowsFrom+1, numCols(), 1) and must + * not overlap this tensor's storage. Indices and shape are checked. + * The copy runs on this tensor's stream; destination.streamIdx() is not + * changed. Keep both buffers alive until completion and explicitly order + * accesses from any other stream, including prior uses of destination. + */ + void getRowsInto(DTensor &destination, size_t rowsFrom, size_t rowsTo, size_t matIdx) const; + /** * Transposes each (m,n)-matrix of this tensor. * Each transposed matrix is stored at same k-index in new tensor. @@ -549,7 +658,7 @@ public: * @param j second column index * @param c cos θ * @param minus_s minus sin θ - * @throws std::invalid_argument if either i or j is greater or equal ncols + * @pre i and j are less than numCols(); indices are not checked. */ void applyRightGivensRotation(size_t i, size_t j, const T *c, const T *minus_s); @@ -560,7 +669,7 @@ public: * @param j second row index * @param c cos θ * @param minus_s minus sin θ - * @throws std::invalid_argument if either i or j is greater or equal nrows + * @pre i and j are less than numRows(); indices are not checked. */ void applyLeftGivensRotation(size_t i, size_t j, const T *c, const T *minus_s); @@ -575,6 +684,21 @@ public: */ void leastSquaresBatched(DTensor &b); + /** + * Allocation-free overload using a caller-owned device status buffer. + * infoArray must have at least numMats() integers; the first numMats() + * entries receive cuBLAS per-matrix status codes. No resizing is performed. + * A is overwritten by its factorisation and b by the solutions/residuals. + * Work is submitted on A's stream. Keep A, b and infoArray alive until that + * stream completes. Reuse on the same stream is ordered; sharing the buffer + * across streams requires explicit synchronization. Its streamIdx is not used. + * Single matrices without pointer arrays require the workspace overload instead. + */ + void leastSquaresBatched(DTensor &b, DTensor &infoArray); + + /** Reuse constructor-allocated status/pointer storage; supports batches of one. */ + void leastSquaresBatched(DTensor &b, LeastSquaresBatchedWorkspace &workspace); + /** * Batched `C <- bC + a*A*B`. * Performs the operation `Ci <- bCi + a*Ai*Bi` for each k-index `i`. @@ -590,9 +714,10 @@ public: /** * Reshapes the tensor * - * If the new number of tensors is larger than the current one, - * this method will allocate a device array of type T* and length - * equal to the new number of matrices. + * The owned matrix-pointer capacity is retained when shrinking. This method + * allocates a device array of type T* only when the required count exceeds + * that capacity. Reshaped views allocate their own matrix-pointer arrays + * when needed, without taking ownership of the element storage. * * No new memory is allocated if newNumMats = 1 * @@ -618,8 +743,12 @@ public: /* ------------- OPERATORS ------------- */ + /** Deep copy, replacing this tensor's storage (also when this is a view). */ DTensor &operator=(const DTensor &other); + /** Transfer ownership/view metadata, leaving other empty. Invalidates old views. */ + DTensor &operator=(DTensor &&other) noexcept; + T operator()(size_t i, size_t j = 0, size_t k = 0) const; DTensor &operator*=(T scalar); @@ -660,15 +789,50 @@ public: } }; /* END OF DTENSOR */ +// Kernel arguments capture pointer values without a temporary host upload buffer. +template +__global__ void k_setSingleMatrixPointers(T **pointers, T *matrix, T *rhs) { + pointers[0] = matrix; + pointers[1] = rhs; +} + +/** Reusable least-squares storage for a fixed shape; all allocation is in construction. + * Keep it alive until the solve stream completes. Concurrent streams need separate workspaces. + */ +template +class LeastSquaresBatchedWorkspace { + friend class DTensor; + size_t m_rows, m_cols, m_batch; + DTensor m_info; + std::unique_ptr> m_singleMatrixPointers; +public: + explicit LeastSquaresBatchedWorkspace(const DTensor &A) + : m_rows(A.numRows()), m_cols(A.numCols()), m_batch(A.numMats()), m_info(A.numMats()) { + if (!m_rows || !m_cols || !m_batch || m_rows < m_cols) + throw std::invalid_argument("[Least squares workspace] A must be nonempty and tall or square"); + m_info.setStreamIdx(A.streamIdx()); + if (m_batch == 1) m_singleMatrixPointers = std::make_unique>(2); + } + + DTensor &info() { return m_info; } +}; + template -DTensor DTensor::setStreamIdx(size_t idx) { - if (idx >= s_numStreams) { +DTensor &DTensor::setStreamIdx(size_t idx) & { + if (idx >= Session::getInstance().numStreams()) { throw std::invalid_argument("Invalid stream index; it exceeds the max allocated streams"); } m_idxStream = idx; return *this; } +template +DTensor DTensor::setStreamIdx(size_t idx) && { + // Validate before moving so a rejected stream leaves the source intact. + setStreamIdx(idx); + return std::move(*this); +} + template void DTensor::initialisePointersToMatricesData() { /* Make sure m_d_ptrMatrices has been allocated */ @@ -684,17 +848,22 @@ void DTensor::initialisePointersToMatricesData() { } /* Upload data to m_d_ptrMatrices */ size_t buffer_size = m_numMats * sizeof(T *); - gpuErrChk(cudaMemcpy(m_d_ptrMatrices, h_pointers.data(), buffer_size, cudaMemcpyHostToDevice)); + auto stream = selectedStream(); + gpuErrChk(cudaMemcpyAsync(m_d_ptrMatrices, h_pointers.data(), buffer_size, cudaMemcpyHostToDevice, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); } template DTensor DTensor::createRandomTensor(size_t numRows, size_t numCols, size_t numMats, T low, T hi) { + const size_t size = checkedTensorSize(numRows, numCols, numMats); + if (!(low <= hi)) throw std::invalid_argument("invalid random range"); if constexpr (std::is_floating_point::value) { - auto randVec = generateRealRandomVector(numRows * numCols * numMats, low, hi); + if (!std::isfinite(low) || !std::isfinite(hi)) throw std::invalid_argument("random bounds must be finite"); + auto randVec = generateRealRandomVector(size, low, hi); DTensor a(randVec, numRows, numCols, numMats); return a; } else if constexpr (std::is_same_v) { - auto randVec = generateIntRandomVector(numRows * numCols * numMats, low, hi); + auto randVec = generateIntRandomVector(size, low, hi); DTensor a(randVec, numRows, numCols, numMats); return a; } @@ -717,45 +886,36 @@ data_t vectorFromTextFile(std::string path_to_file) { file.open(path_to_file, std::ios::in); if (!file.is_open()) { throw std::invalid_argument("[vectorFromTextFile] the file does not exist"); } - std::string line; - getline(file, line); - dataStruct.numRows = atoi(line.c_str()); - getline(file, line); - dataStruct.numCols = atoi(line.c_str()); - getline(file, line); - dataStruct.numMats = atoi(line.c_str()); - - size_t numElements = dataStruct.numRows * dataStruct.numCols * dataStruct.numMats; - std::vector vecDataFromFile(numElements); - - size_t i = 0; - while (getline(file, line)) { - if constexpr (std::is_same_v) { - vecDataFromFile[i] = atoi(line.c_str()); - } else if constexpr (std::is_same_v) { - vecDataFromFile[i] = std::stod(line.c_str()); - } else if constexpr (std::is_same_v) { - vecDataFromFile[i] = std::stof(line.c_str()); - } else if constexpr (std::is_same_v) { - vecDataFromFile[i] = std::stold(line.c_str()); - } else if constexpr (std::is_same_v) { - vecDataFromFile[i] = std::stol(line.c_str()); - } else if constexpr (std::is_same_v) { - vecDataFromFile[i] = std::stoll(line.c_str()); - } else if constexpr (std::is_same_v) { - vecDataFromFile[i] = std::stoul(line.c_str()); - } else if constexpr (std::is_same_v) { - vecDataFromFile[i] = std::stoull(line.c_str()); - } else if constexpr (std::is_same_v) { - sscanf(line.c_str(), "%zu", &vecDataFromFile[i]); - } else { - throw std::invalid_argument("data type not supported"); + auto readDimension = [&file]() -> size_t { + std::string line; + if (!std::getline(file, line)) throw std::invalid_argument("missing tensor dimension"); + if (!line.empty() && line.back() == '\r') line.pop_back(); + if (line.empty() || line.find_first_not_of("0123456789") != std::string::npos) + throw std::invalid_argument("invalid tensor dimension in text file"); + try { + unsigned long long value = std::stoull(line); + if (value > std::numeric_limits::max()) throw std::invalid_argument("dimension too large"); + return static_cast(value); + } catch (const std::out_of_range &) { throw std::invalid_argument("dimension too large"); } + }; + dataStruct.numRows = readDimension(); + dataStruct.numCols = readDimension(); + dataStruct.numMats = readDimension(); + size_t numElements = checkedTensorSize(dataStruct.numRows, dataStruct.numCols, dataStruct.numMats); + if constexpr (!std::is_arithmetic_v || std::is_same_v) { + throw std::invalid_argument("data type not supported"); + } else { + dataStruct.data.resize(numElements); + for (T &value : dataStruct.data) { + std::string line; + if (!std::getline(file, line)) throw std::invalid_argument("truncated tensor text file"); + std::istringstream input(line); + if (!(input >> value) || !(input >> std::ws).eof()) + throw std::invalid_argument("invalid tensor text data"); } - - if (++i == numElements) break; + std::string extra; + if (file >> extra) throw std::invalid_argument("extra tensor text data"); } - dataStruct.data = vecDataFromFile; - file.close(); return dataStruct; } @@ -766,27 +926,38 @@ data_t vectorFromBinaryFile(std::string path_to_file) { std::ifstream inFile; inFile.open(path_to_file, std::ios::binary); if (!inFile.is_open()) { throw std::invalid_argument("[vectorFromBinaryFile] the file does not exist"); } - inFile.read(reinterpret_cast(&(dataStruct.numRows)), sizeof(uint64_t)); - inFile.read(reinterpret_cast(&(dataStruct.numCols)), sizeof(uint64_t)); - inFile.read(reinterpret_cast(&(dataStruct.numMats)), sizeof(uint64_t)); - uint64_t numElements = dataStruct.numRows * dataStruct.numCols * dataStruct.numMats; - std::vector vecDataFromFile(numElements); - for (size_t i = 0; i < numElements; i++) { - T el; - inFile.read(reinterpret_cast(&el), sizeof(T)); - vecDataFromFile[i] = el; - } - inFile.close(); - dataStruct.data = vecDataFromFile; + uint64_t dimensions[3]{}; + if (!inFile.read(reinterpret_cast(dimensions), sizeof(dimensions))) + throw std::invalid_argument("truncated tensor binary header"); + for (auto dimension : dimensions) + if (dimension > std::numeric_limits::max()) throw std::invalid_argument("dimension too large"); + dataStruct.numRows = static_cast(dimensions[0]); + dataStruct.numCols = static_cast(dimensions[1]); + dataStruct.numMats = static_cast(dimensions[2]); + size_t numElements = checkedTensorSize(dataStruct.numRows, dataStruct.numCols, dataStruct.numMats); + size_t bytes = numElements * sizeof(T); + if (bytes > static_cast(std::numeric_limits::max())) + throw std::invalid_argument("tensor binary payload too large"); + auto payloadStart = inFile.tellg(); + inFile.seekg(0, std::ios::end); + auto payloadEnd = inFile.tellg(); + if (payloadStart < 0 || payloadEnd < payloadStart || + static_cast(payloadEnd - payloadStart) != bytes) + throw std::invalid_argument("invalid tensor binary payload size"); + inFile.seekg(payloadStart); + dataStruct.data.resize(numElements); + if (!inFile.read(reinterpret_cast(dataStruct.data.data()), static_cast(bytes))) + throw std::invalid_argument("truncated tensor binary payload"); return dataStruct; } template DTensor DTensor::parseFromFile(std::string path_to_file, StorageMode mode) { + validateStorageMode(mode); // Figure out file extension size_t pathToFileLength = path_to_file.length(); - std::string fileNameExtension = path_to_file.substr(pathToFileLength - 3); + std::string fileNameExtension = pathToFileLength >= 3 ? path_to_file.substr(pathToFileLength - 3) : ""; typedef data_t (*PARSER)(std::string); PARSER parser = (fileNameExtension == ".bt") ? vectorFromBinaryFile : vectorFromTextFile; auto parsedData = parser(path_to_file); @@ -801,7 +972,7 @@ void DTensor::saveToFile(std::string pathToFile) { // Figure out file extension size_t pathToFileLength = pathToFile.length(); - std::string fileNameExtension = pathToFile.substr(pathToFileLength - 3); + std::string fileNameExtension = pathToFileLength >= 3 ? pathToFile.substr(pathToFileLength - 3) : ""; // If the extension is .bt... if (fileNameExtension == ".bt") { uint64_t nr = (uint64_t) numRows(), @@ -809,18 +980,23 @@ void DTensor::saveToFile(std::string pathToFile) { nm = (uint64_t) numMats(); std::ofstream outFile; outFile.open(pathToFile, std::ios::binary); + if (!outFile) throw std::runtime_error("cannot open tensor output file"); outFile.write(reinterpret_cast(&nr), sizeof(uint64_t)); outFile.write(reinterpret_cast(&nc), sizeof(uint64_t)); outFile.write(reinterpret_cast(&nm), sizeof(uint64_t)); for (const T &el: myData) outFile.write(reinterpret_cast(&el), sizeof(T)); outFile.close(); + if (!outFile) throw std::runtime_error("failed writing tensor output file"); } else { std::ofstream file(pathToFile); + if (!file) throw std::runtime_error("cannot open tensor output file"); file << numRows() << std::endl << numCols() << std::endl << numMats() << std::endl; if constexpr (std::is_floating_point::value) { file << std::setprecision(std::numeric_limits::max_digits10); } for (const T &el: myData) file << el << std::endl; + file.close(); + if (!file) throw std::runtime_error("failed writing tensor output file"); } } @@ -828,7 +1004,7 @@ void DTensor::saveToFile(std::string pathToFile) { template void DTensor::reshape(size_t newNumRows, size_t newNumCols, size_t newNumMats) { if (m_numRows == newNumRows && m_numCols == newNumCols && m_numMats == newNumMats) return; - size_t newNumElements = newNumRows * newNumCols * newNumMats; + size_t newNumElements = checkedTensorSize(newNumRows, newNumCols, newNumMats); /* Check whether dimensions are compatible */ if (numEl() != newNumElements) { char errMessage[256]; @@ -838,20 +1014,18 @@ void DTensor::reshape(size_t newNumRows, size_t newNumCols, size_t newNumMats throw std::invalid_argument(errMessage); } - /* Only free/reallocate if newNumMats > m_numMats - * otherwise, reuse the already allocated memory space */ - if (newNumMats > m_numMats) { - /* Free the memory for m_d_ptrMatrices */ - if (m_d_ptrMatrices && m_doDestroyPtrMatrices) { - gpuErrChk(cudaFree(m_d_ptrMatrices)); - m_d_ptrMatrices = nullptr; - m_doDestroyPtrMatrices = false; - } - /* Reallocate memory for m_d_ptrMatrices, if necessary */ - if (newNumMats > 1) { - gpuErrChk(Session::getInstance().cudaAllocate((void**)&m_d_ptrMatrices, newNumMats * sizeof(T *))); - m_doDestroyPtrMatrices = true; - } + // Capacity belongs to the allocation, not the current shape. A view gets + // its own pointer array when reshaped; its borrowed element data stays borrowed. + if (newNumMats > 1 && newNumMats > m_matrixPointerCapacity) { + T **newPointers = nullptr; + gpuErrChk(Session::getInstance().cudaAllocate((void**)&newPointers, newNumMats * sizeof(T *))); + if (m_doDestroyPtrMatrices) + gpuErrChk(Session::getInstance().cudaRelease(m_d_ptrMatrices, m_matrixPointerCapacity * sizeof(T *))); + m_d_ptrMatrices = newPointers; + m_matrixPointerCapacity = newNumMats; + m_doDestroyPtrMatrices = true; + } else if (newNumMats == 1 && !m_doDestroyPtrMatrices) { + m_d_ptrMatrices = nullptr; } m_numRows = newNumRows; @@ -862,25 +1036,35 @@ void DTensor::reshape(size_t newNumRows, size_t newNumCols, size_t newNumMats template DTensor::DTensor(size_t m, size_t n, size_t k, bool zero) { + size_t size = checkedTensorSize(m, n, k); m_numRows = m; m_numCols = n; m_numMats = k; - size_t size = m * n * k; - allocateOnDevice(size, zero); - /* Initialise m_d_ptrMatrices */ - initialisePointersToMatricesData(); + try { + allocateOnDevice(size, zero); + initialisePointersToMatricesData(); + } catch (...) { + destroy(); + throw; + } } template DTensor::DTensor(const std::vector &data, size_t m, size_t n, size_t k, StorageMode mode) { + size_t size = checkedTensorSize(m, n, k); + validateStorageMode(mode); + if (data.size() != size) throw std::invalid_argument("tensor data has wrong size"); m_numRows = m; m_numCols = n; m_numMats = k; - size_t size = m * n * k; - allocateOnDevice(size); - upload(data, mode); - /* Initialise m_d_ptrMatrices */ - initialisePointersToMatricesData(); + try { + allocateOnDevice(size); + upload(data, mode); + initialisePointersToMatricesData(); + } catch (...) { + destroy(); + throw; + } } template @@ -889,16 +1073,28 @@ DTensor::DTensor(const DTensor &other) { m_numRows = other.m_numRows; m_numCols = other.m_numCols; m_idxStream = other.m_idxStream; - allocateOnDevice(m_numRows * m_numCols * m_numMats); - gpuErrChk(cudaMemcpy(m_d_data, other.raw(), m_numRows * m_numCols * m_numMats * sizeof(T), - cudaMemcpyDeviceToDevice)); - /* Initialise m_d_ptrMatrices */ - initialisePointersToMatricesData(); + if (other.numEl() == 0) return; + try { + allocateOnDevice(numEl()); + initialisePointersToMatricesData(); + cudaStream_t stream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(m_idxStream), &stream)); + gpuErrChk(cudaMemcpyAsync(m_d_data, other.raw(), numEl() * sizeof(T), + cudaMemcpyDeviceToDevice, stream)); + // Copies are complete on return, including work previously queued on + // the source's stream, so the source can subsequently be destroyed. + gpuErrChk(cudaStreamSynchronize(stream)); + } catch (...) { + destroy(); + throw; + } } template DTensor::DTensor(const DTensor &other, size_t axis, size_t from, size_t to) { - if (from > to) throw std::invalid_argument("from > to"); + if (axis > 2 || from > to) throw std::invalid_argument("invalid slice axis or range"); + const size_t extent = axis == 0 ? other.m_numRows : axis == 1 ? other.m_numCols : other.m_numMats; + if (to >= extent) throw std::invalid_argument("slice range exceeds axis"); size_t offset = 0, len = to - from + 1; m_d_ptrMatrices = nullptr; if (axis == 2) { @@ -906,7 +1102,7 @@ DTensor::DTensor(const DTensor &other, size_t axis, size_t from, size_t to m_numRows = other.m_numRows; m_numCols = other.m_numCols; m_numMats = len; - m_d_ptrMatrices = other.m_d_ptrMatrices + from; + m_d_ptrMatrices = other.m_d_ptrMatrices ? other.m_d_ptrMatrices + from : nullptr; } else if (axis == 1) { offset = other.m_numRows * from; m_numRows = other.m_numRows; @@ -925,24 +1121,8 @@ DTensor::DTensor(const DTensor &other, size_t axis, size_t from, size_t to } template -DTensor::DTensor(DTensor &&other) { - /* Steal everything from other */ - m_numCols = other.m_numCols; - m_numRows = other.m_numRows; - m_numMats = other.m_numMats; - m_d_data = other.m_d_data; - m_doDestroyData = other.m_doDestroyData; - m_doDestroyPtrMatrices = other.m_doDestroyPtrMatrices; - m_d_ptrMatrices = other.m_d_ptrMatrices; - m_idxStream = other.m_idxStream; - /* Invalidate other */ - other.m_doDestroyPtrMatrices = false; - other.m_doDestroyData = false; - other.m_d_ptrMatrices = nullptr; - other.m_d_data = nullptr; - other.m_numCols = 0; - other.m_numRows = 0; - other.m_numMats = 0; +DTensor::DTensor(DTensor &&other) noexcept { + swap(other); } template @@ -1034,7 +1214,7 @@ inline float DTensor::maxAbs() const { gpuErrChk( cublasIsamax(Session::getInstance().cuBlasHandle(m_idxStream), m_numRows * m_numCols * m_numMats, m_d_data, 1, &idx)); - gpuErrChk(cudaMemcpy(&hostDst, m_d_data + idx - 1, sizeof(float), cudaMemcpyDeviceToHost)); + hostDst = downloadElement(idx - 1); return std::signbit(hostDst) ? -hostDst : hostDst; } @@ -1045,7 +1225,7 @@ inline double DTensor::maxAbs() const { gpuErrChk( cublasIdamax(Session::getInstance().cuBlasHandle(m_idxStream), m_numRows * m_numCols * m_numMats, m_d_data, 1, &idx)); - gpuErrChk(cudaMemcpy(&hostDst, m_d_data + idx - 1, sizeof(double), cudaMemcpyDeviceToHost)); + hostDst = downloadElement(idx - 1); return std::signbit(hostDst) ? -hostDst : hostDst; } @@ -1056,7 +1236,7 @@ inline float DTensor::minAbs() const { gpuErrChk( cublasIsamin(Session::getInstance().cuBlasHandle(m_idxStream), m_numRows * m_numCols * m_numMats, m_d_data, 1, &idx)); - gpuErrChk(cudaMemcpy(&hostDst, m_d_data + idx - 1, sizeof(float), cudaMemcpyDeviceToHost)); + hostDst = downloadElement(idx - 1); return std::signbit(hostDst) ? -hostDst : hostDst; } @@ -1067,7 +1247,7 @@ inline double DTensor::minAbs() const { gpuErrChk( cublasIdamin(Session::getInstance().cuBlasHandle(m_idxStream), m_numRows * m_numCols * m_numMats, m_d_data, 1, &idx)); - gpuErrChk(cudaMemcpy(&hostDst, m_d_data + idx - 1, sizeof(double), cudaMemcpyDeviceToHost)); + hostDst = downloadElement(idx - 1); return std::signbit(hostDst) ? -hostDst : hostDst; } @@ -1111,15 +1291,24 @@ inline void DTensor::allocateOnDevice(size_t size, bool zero) { m_doDestroyData = true; size_t buffer_size = size * sizeof(T); gpuErrChk(Session::getInstance().cudaAllocate((void**)&m_d_data, buffer_size)); - if (zero) gpuErrChk(cudaMemset(m_d_data, 0, buffer_size)); // set to zero all elements + m_dataBytes = buffer_size; + if (zero) { + auto stream = selectedStream(); + gpuErrChk(cudaMemsetAsync(m_d_data, 0, buffer_size, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); + } if (numMats() > 1) { m_doDestroyPtrMatrices = true; cudaStatus = Session::getInstance().cudaAllocate((void **) &m_d_ptrMatrices, numMats() * sizeof(T *)); if (cudaStatus != cudaSuccess) { - gpuErrChk(cudaFree(m_d_data)); // ... free previously allocated memory + gpuErrChk(Session::getInstance().cudaRelease(m_d_data, m_dataBytes)); + m_d_data = nullptr; + m_dataBytes = 0; + m_doDestroyData = false; gpuErrChk(cudaStatus); // ... and memento mori } + m_matrixPointerCapacity = numMats(); } else { m_doDestroyPtrMatrices = false; } @@ -1127,19 +1316,21 @@ inline void DTensor::allocateOnDevice(size_t size, bool zero) { template inline bool DTensor::upload(const std::vector &vec, StorageMode mode) { + validateStorageMode(mode); size_t size = vec.size(); size_t thisSize = m_numRows * m_numCols * m_numMats; // make sure vec is of right size if (size != thisSize) throw std::invalid_argument("[upload] vec has wrong size"); - std::vector vecCm(thisSize); + if (size == 0) return true; + auto stream = selectedStream(); if (mode == StorageMode::rowMajor) { + std::vector vecCm(thisSize); rm2cm(vec, vecCm); + gpuErrChk(cudaMemcpyAsync(m_d_data, vecCm.data(), size * sizeof(T), cudaMemcpyHostToDevice, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); // Keep the conversion buffer alive. } else { - vecCm = vec; - } - if (size <= thisSize) { - size_t buffer_size = size * sizeof(T); - gpuErrChk(cudaMemcpy(m_d_data, vecCm.data(), buffer_size, cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpyAsync(m_d_data, vec.data(), size * sizeof(T), cudaMemcpyHostToDevice, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); } return true; } @@ -1147,10 +1338,13 @@ inline bool DTensor::upload(const std::vector &vec, StorageMode mode) { template inline void DTensor::download(std::vector &vec) const { vec.resize(m_numRows * m_numCols * m_numMats); - gpuErrChk(cudaMemcpy(vec.data(), + if (vec.empty()) return; + auto stream = selectedStream(); + gpuErrChk(cudaMemcpyAsync(vec.data(), m_d_data, m_numRows * m_numCols * m_numMats * sizeof(T), - cudaMemcpyDeviceToHost)); + cudaMemcpyDeviceToHost, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); } template @@ -1167,6 +1361,7 @@ inline T **DTensor::ptrMatrices() const { template<> inline DTensor DTensor::tr() const { DTensor transposes(m_numCols, m_numRows, m_numMats); + transposes.m_idxStream = m_idxStream; float alpha = 1.0f, beta = 0; size_t numElMat = m_numCols * m_numRows; for (size_t i = 0; i < m_numMats; i++) { @@ -1183,6 +1378,7 @@ inline DTensor DTensor::tr() const { template<> inline DTensor DTensor::tr() const { DTensor transposes(m_numCols, m_numRows, m_numMats); + transposes.m_idxStream = m_idxStream; double alpha = 1.0f, beta = 0; size_t numElMat = m_numCols * m_numRows; for (size_t i = 0; i < m_numMats; i++) { @@ -1201,10 +1397,15 @@ inline void DTensor::deviceCopyTo(DTensor &elsewhere) const { if (elsewhere.numEl() < numEl()) { throw std::invalid_argument("[deviceCopyTo] tensor does not fit into destination"); } - gpuErrChk(cudaMemcpy(elsewhere.raw(), - m_d_data, - m_numRows * m_numCols * m_numMats * sizeof(T), - cudaMemcpyDeviceToDevice)); + if (numEl() == 0) return; + auto stream = selectedStream(); + auto destinationStream = elsewhere.selectedStream(); + if (destinationStream != stream) gpuErrChk(cudaStreamSynchronize(destinationStream)); + if (elsewhere.raw() != m_d_data) { + gpuErrChk(cudaMemcpyAsync(elsewhere.raw(), m_d_data, + m_numRows * m_numCols * m_numMats * sizeof(T), cudaMemcpyDeviceToDevice, stream)); + } + gpuErrChk(cudaStreamSynchronize(stream)); } template<> @@ -1218,12 +1419,19 @@ inline DTensor &DTensor::operator*=(double scalar) { template DTensor &DTensor::operator=(const DTensor &other) { - m_numMats = other.m_numMats; - m_numRows = other.m_numRows; - m_numCols = other.m_numCols; - m_doDestroyData = false; - m_d_data = other.m_d_data; - m_idxStream = other.m_idxStream; + if (this != &other) { + DTensor temporary(other); + swap(temporary); + } + return *this; +} + +template +DTensor &DTensor::operator=(DTensor &&other) noexcept { + if (this != &other) { + DTensor temporary(std::move(other)); + swap(temporary); + } return *this; } @@ -1277,10 +1485,10 @@ inline DTensor &DTensor::operator-=(const DTensor &rhs) template inline T DTensor::operator()(size_t i, size_t j, size_t k) const { - T hostDst; + if (i >= m_numRows || j >= m_numCols || k >= m_numMats) + throw std::invalid_argument("tensor index out of range"); size_t offset = i + m_numRows * (j + m_numCols * k); - gpuErrChk(cudaMemcpy(&hostDst, m_d_data + offset, sizeof(T), cudaMemcpyDeviceToHost)); - return hostDst; + return downloadElement(offset); } template<> @@ -1338,7 +1546,8 @@ inline void DTensor::addAB(const DTensor &A, const DTensor } template<> -inline void DTensor::leastSquaresBatched(DTensor &B) { +inline void DTensor::leastSquaresBatchedImpl(DTensor &B, DTensor &infoArray, + double **matrixPointers, double **rhsPointers) { size_t batchSize = numMats(); size_t nColsB = B.numCols(); if (B.numRows() != m_numRows) @@ -1350,15 +1559,16 @@ inline void DTensor::leastSquaresBatched(DTensor &B) { if (m_numCols > m_numRows) throw std::invalid_argument("[Least squares batched] supports square or tall matrices only"); int info = 0; - DTensor infoArray(batchSize); // TODO consider preallocating? + if (infoArray.numEl() < batchSize || !infoArray.raw()) + throw std::invalid_argument("[Least squares batched] status buffer is too small"); gpuErrChk(cublasDgelsBatched(Session::getInstance().cuBlasHandle(m_idxStream), CUBLAS_OP_N, m_numRows, m_numCols, nColsB, - m_d_ptrMatrices, + matrixPointers, m_numRows, - B.m_d_ptrMatrices, + rhsPointers, m_numRows, &info, infoArray.raw(), @@ -1366,7 +1576,8 @@ inline void DTensor::leastSquaresBatched(DTensor &B) { } template<> -inline void DTensor::leastSquaresBatched(DTensor &B) { +inline void DTensor::leastSquaresBatchedImpl(DTensor &B, DTensor &infoArray, + float **matrixPointers, float **rhsPointers) { size_t batchSize = numMats(); size_t nColsB = B.numCols(); if (B.numRows() != m_numRows) @@ -1378,48 +1589,91 @@ inline void DTensor::leastSquaresBatched(DTensor &B) { if (m_numCols > m_numRows) throw std::invalid_argument("[Least squares batched] supports square or tall matrices only"); int info = 0; - DTensor infoArray(batchSize); // TODO consider preallocating? + if (infoArray.numEl() < batchSize || !infoArray.raw()) + throw std::invalid_argument("[Least squares batched] status buffer is too small"); gpuErrChk(cublasSgelsBatched(Session::getInstance().cuBlasHandle(m_idxStream), CUBLAS_OP_N, m_numRows, m_numCols, nColsB, - m_d_ptrMatrices, + matrixPointers, m_numRows, - B.m_d_ptrMatrices, + rhsPointers, m_numRows, &info, infoArray.raw(), batchSize)); } +template +inline void DTensor::leastSquaresBatched(DTensor &B, DTensor &infoArray) { + if (numMats() == 1 && (!m_d_ptrMatrices || !B.m_d_ptrMatrices)) + throw std::invalid_argument("[Least squares batched] batch of one requires LeastSquaresBatchedWorkspace"); + leastSquaresBatchedImpl(B, infoArray, m_d_ptrMatrices, B.m_d_ptrMatrices); +} + +template +inline void DTensor::leastSquaresBatched(DTensor &B, LeastSquaresBatchedWorkspace &workspace) { + if (numRows() != workspace.m_rows || numCols() != workspace.m_cols || numMats() != workspace.m_batch) + throw std::invalid_argument("[Least squares batched] workspace shape does not match A"); + if (B.numRows() != numRows() || B.numCols() != 1 || B.numMats() != numMats()) + throw std::invalid_argument("[Least squares batched] incompatible rhs"); + T **matrixPointers = m_d_ptrMatrices, **rhsPointers = B.m_d_ptrMatrices; + workspace.m_info.setStreamIdx(streamIdx()); + if (workspace.m_singleMatrixPointers) { + cudaStream_t stream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(streamIdx()), &stream)); + matrixPointers = workspace.m_singleMatrixPointers->raw(); + rhsPointers = matrixPointers + 1; + k_setSingleMatrixPointers<<<1, 1, 0, stream>>>(matrixPointers, raw(), B.raw()); + gpuErrChk(cudaGetLastError()); + } + leastSquaresBatchedImpl(B, workspace.m_info, matrixPointers, rhsPointers); +} + +template +inline void DTensor::leastSquaresBatched(DTensor &B) { + LeastSquaresBatchedWorkspace workspace(*this); + leastSquaresBatched(B, workspace); +} + +template +inline void DTensor::getRowsInto(DTensor &destination, size_t rowsFrom, size_t rowsTo, size_t matIdx) const { + if (rowsFrom > rowsTo || rowsTo >= numRows() || matIdx >= numMats()) + throw std::invalid_argument("invalid row or matrix range"); + const size_t count = rowsTo - rowsFrom + 1; + if (destination.numRows() != count || destination.numCols() != numCols() || destination.numMats() != 1) + throw std::invalid_argument("row destination has incompatible shape"); + // CUDA copies do not support overlapping source and destination regions. + const auto sourceBegin = reinterpret_cast(raw()); + const auto destinationBegin = reinterpret_cast(destination.raw()); + if (sourceBegin < destinationBegin + destination.numEl() * sizeof(T) && + destinationBegin < sourceBegin + numEl() * sizeof(T)) + throw std::invalid_argument("row destination overlaps source storage"); + cudaStream_t stream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(m_idxStream), &stream)); + gpuErrChk(cudaMemcpy2DAsync(destination.raw(), count * sizeof(T), + raw() + rowsFrom + matIdx * numCols() * numRows(), numRows() * sizeof(T), + count * sizeof(T), numCols(), cudaMemcpyDeviceToDevice, stream)); +} + template<> inline DTensor DTensor::getRows(size_t rowsFrom, size_t rowsTo, size_t matIdx) const { - size_t rowsRangeLength = rowsTo - rowsFrom + 1; - size_t n = numCols(), m = numRows(); - DTensor rowsOnly(rowsRangeLength, numCols(), 1); - for (size_t i = 0; i < rowsRangeLength; i++) { - gpuErrChk(cublasDcopy(Session::getInstance().cuBlasHandle(m_idxStream), - n, // # values to copy - raw() + rowsFrom + i + matIdx * n * m, m, - rowsOnly.raw() + i, - rowsRangeLength)); - } + if (rowsFrom > rowsTo || rowsTo >= numRows() || matIdx >= numMats()) + throw std::invalid_argument("invalid row or matrix range"); + DTensor rowsOnly(rowsTo - rowsFrom + 1, numCols(), 1); + rowsOnly.m_idxStream = m_idxStream; + getRowsInto(rowsOnly, rowsFrom, rowsTo, matIdx); return rowsOnly; } template<> inline DTensor DTensor::getRows(size_t rowsFrom, size_t rowsTo, size_t matIdx) const { - size_t rowsRangeLength = rowsTo - rowsFrom + 1; - size_t n = numCols(), m = numRows(); - DTensor rowsOnly(rowsRangeLength, numCols(), 1); - for (size_t i = 0; i < rowsRangeLength; i++) { - gpuErrChk(cublasScopy(Session::getInstance().cuBlasHandle(m_idxStream), - n, // # values to copy - raw() + rowsFrom + i + matIdx * n * m, m, - rowsOnly.raw() + i, - rowsRangeLength)); - } + if (rowsFrom > rowsTo || rowsTo >= numRows() || matIdx >= numMats()) + throw std::invalid_argument("invalid row or matrix range"); + DTensor rowsOnly(rowsTo - rowsFrom + 1, numCols(), 1); + rowsOnly.m_idxStream = m_idxStream; + getRowsInto(rowsOnly, rowsFrom, rowsTo, matIdx); return rowsOnly; } @@ -1490,6 +1744,17 @@ __global__ void k_countNonzeroSingularValues(const T *d_array, size_t n, unsigne } } +TEMPLATE_WITH_TYPE_T +TEMPLATE_CONSTRAINT_REQUIRES_FPX +__global__ void k_countRelativeSingularValues(const T *singularValues, size_t n, + unsigned int *count, T rtol, T atol) { + const size_t idx = threadIdx.x + size_t(blockIdx.x) * blockDim.x; + // gesvd returns singular values in descending order, separately per matrix. + const T relativeThreshold = rtol * singularValues[0]; + const T threshold = atol > relativeThreshold ? atol : relativeThreshold; + if (idx < n && singularValues[idx] > threshold) atomicAdd(count, 1); +} + /** * Singular value decomposition (SVD) needs a workspace to be setup for cuSolver before factorisation. @@ -1502,21 +1767,21 @@ TEMPLATE_CONSTRAINT_REQUIRES_FPX class Svd : public IStatus { private: int m_lwork = -1; ///< Size of workspace needed for SVD - DTensor *m_tensor = nullptr; ///< Pointer to original matrix to be factorised + std::unique_ptr > m_ownedTensor; ///< Optional private input copy + DTensor *m_tensor = nullptr; ///< Observer of caller input or m_ownedTensor std::shared_ptr > m_Vtr; ///< Matrix V' or right singular vectors std::shared_ptr > m_S; ///< Diagonal matrix S or singular values std::shared_ptr > m_U; ///< Matrix U or left singular vectors std::unique_ptr > m_workspace; ///< Workspace for SVD std::shared_ptr > m_rank; ///< Rank of original matrix bool m_computeU = false; ///< Whether to compute U - bool m_destroyMatrix = true; ///< Whether to sacrifice original matrix /** * Ensures tensor to factorise contains exactly one matrix, and that matrix is tall. * @param mat matrix to be factorised */ void checkMatrix(DTensor &tensor) const { - if (tensor.numRows() < tensor.numCols()) { + if (!tensor.numRows() || !tensor.numCols() || !tensor.numMats() || tensor.numRows() < tensor.numCols()) { throw std::invalid_argument("[svd] your matrix is fat (no offence)"); } }; @@ -1539,8 +1804,13 @@ public: bool computeU = false, bool destroyMatrix = true) : IStatus(mat.numMats()) { checkMatrix(mat); - m_destroyMatrix = destroyMatrix; - m_tensor = (destroyMatrix) ? &mat : new DTensor(mat); + if (destroyMatrix) { + m_tensor = &mat; + } else { + m_ownedTensor = std::make_unique >(mat); + m_tensor = m_ownedTensor.get(); + } + m_info->setStreamIdx(mat.streamIdx()); m_computeU = computeU; size_t m = mat.numRows(); size_t n = mat.numCols(); @@ -1551,7 +1821,11 @@ public: m_Vtr = std::make_shared >(n, n, nMats); m_S = std::make_shared >(minMN, 1, nMats); m_rank = std::make_unique >(1, 1, nMats, true); + m_rank->setStreamIdx(mat.streamIdx()); + m_S->setStreamIdx(mat.streamIdx()); + m_Vtr->setStreamIdx(mat.streamIdx()); if (computeU) m_U = std::make_shared >(m, m, nMats); + if (computeU) m_U->setStreamIdx(mat.streamIdx()); } /** @@ -1582,28 +1856,58 @@ public: return m_U; } - /** - * Destroyer. - */ - ~Svd() { - m_lwork = -1; - if (!m_destroyMatrix && m_tensor) delete m_tensor; - } + // Member ownership also releases the private copy if construction throws. + ~Svd() = default; /** * Computes the rank of the original matrices and returns * a (1, 1, nMats)-tensor. + * Each call recomputes the counts on the input matrix's stream. The result + * is reusable internal storage; copy/download it before a later call if a + * snapshot is needed, and order cross-stream consumers explicitly. * @param epsilon any numerical value less than epsilon is considered zero * @return rank of original matrices */ DTensor const &rank(T epsilon = 1e-6) const { + cudaStream_t stream; + const size_t streamIdx = m_tensor->streamIdx(); + gpuErrChk(cusolverDnGetStream(Session::getInstance().cuSolverHandle(streamIdx), &stream)); + m_rank->setStreamIdx(streamIdx); + gpuErrChk(cudaMemsetAsync(m_rank->raw(), 0, m_rank->numEl() * sizeof(unsigned int), stream)); size_t numElS = m_S->numCols() * m_S->numRows(); for (size_t i = 0; i < m_rank->numMats(); i++) { DTensor Si(*m_S, 2, i, i); DTensor rankI(*m_rank, 2, i, i); - k_countNonzeroSingularValues<<>>(Si.raw(), numElS, + k_countNonzeroSingularValues<<>>(Si.raw(), numElS, rankI.raw(), epsilon); + gpuErrChk(cudaGetLastError()); + } + return *m_rank; + } + + /** + * Count singular values strictly greater than max(atol, rtol*sigmaMax), + * using each matrix's own largest singular value on the device. + * Requires a successful factorise(). Tolerances must be finite and nonnegative. + * atol=0 gives scale-invariant relative rank within floating-point limits; + * a positive atol intentionally adds an absolute cutoff. rtol is explicit. + * Shares reusable result storage with rank(): later calls overwrite it. + * Work is enqueued on the input matrix's stream without host downloads. + */ + DTensor const &rankRelative(T rtol, T atol = T(0)) const { + if (!std::isfinite(rtol) || !std::isfinite(atol) || rtol < T(0) || atol < T(0)) + throw std::invalid_argument("rank tolerances must be finite and nonnegative"); + cudaStream_t stream; + const size_t streamIdx = m_tensor->streamIdx(); + gpuErrChk(cusolverDnGetStream(Session::getInstance().cuSolverHandle(streamIdx), &stream)); + m_rank->setStreamIdx(streamIdx); + gpuErrChk(cudaMemsetAsync(m_rank->raw(), 0, m_rank->numEl() * sizeof(unsigned int), stream)); + const size_t numElS = m_S->numRows() * m_S->numCols(); + for (size_t i = 0; i < m_rank->numMats(); ++i) { + k_countRelativeSingularValues<<>>( + m_S->raw() + i * numElS, numElS, m_rank->raw() + i, rtol, atol); + gpuErrChk(cudaGetLastError()); } return *m_rank; } @@ -1612,12 +1916,12 @@ public: template<> inline void Svd::computeWorkspaceSize(size_t m, size_t n) { - gpuErrChk(cusolverDnSgesvd_bufferSize(Session::getInstance().cuSolverHandle(), m, n, &m_lwork)); + gpuErrChk(cusolverDnSgesvd_bufferSize(Session::getInstance().cuSolverHandle(m_tensor->streamIdx()), m, n, &m_lwork)); } template<> inline void Svd::computeWorkspaceSize(size_t m, size_t n) { - gpuErrChk(cusolverDnDgesvd_bufferSize(Session::getInstance().cuSolverHandle(), m, n, &m_lwork)); + gpuErrChk(cusolverDnDgesvd_bufferSize(Session::getInstance().cuSolverHandle(m_tensor->streamIdx()), m, n, &m_lwork)); } @@ -1634,7 +1938,7 @@ inline void Svd::factorise() { if (m_computeU) Ui = std::make_unique >(*m_U, 2, i, i); gpuErrChk( - cusolverDnDgesvd(Session::getInstance().cuSolverHandle(), + cusolverDnDgesvd(Session::getInstance().cuSolverHandle(m_tensor->streamIdx()), (m_computeU) ? 'A' : 'N', 'A', m, n, Ai.raw(), m, @@ -1661,7 +1965,7 @@ inline void Svd::factorise() { if (m_computeU) Ui = std::make_unique >(*m_U, 2, i, i); gpuErrChk( - cusolverDnSgesvd(Session::getInstance().cuSolverHandle(), + cusolverDnSgesvd(Session::getInstance().cuSolverHandle(m_tensor->streamIdx()), (m_computeU) ? 'A' : 'N', 'A', m, n, Ai.raw(), m, @@ -1701,9 +2005,11 @@ private: public: CholeskyFactoriser(DTensor &A) : IStatus() { + if (!A.numRows() || !A.numCols() || !A.numMats()) throw std::invalid_argument("[Cholesky] empty tensor"); if (A.numMats() > 1) throw std::invalid_argument("[Cholesky] 3D tensors require `CholeskyBatchFactoriser`"); if (A.numRows() != A.numCols()) throw std::invalid_argument("[Cholesky] Matrix A must be square"); m_matrix = &A; + m_info->setStreamIdx(A.streamIdx()); computeWorkspaceSize(); m_workspace = std::make_unique >(m_workspaceSize); } @@ -1726,7 +2032,7 @@ public: template<> inline void CholeskyFactoriser::computeWorkspaceSize() { size_t n = m_matrix->numRows(); - gpuErrChk(cusolverDnDpotrf_bufferSize(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnDpotrf_bufferSize(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, n, nullptr, n, &m_workspaceSize)); } @@ -1734,7 +2040,7 @@ inline void CholeskyFactoriser::computeWorkspaceSize() { template<> inline void CholeskyFactoriser::computeWorkspaceSize() { size_t n = m_matrix->numRows(); - gpuErrChk(cusolverDnSpotrf_bufferSize(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnSpotrf_bufferSize(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, n, nullptr, n, &m_workspaceSize)); } @@ -1742,7 +2048,7 @@ inline void CholeskyFactoriser::computeWorkspaceSize() { template<> inline void CholeskyFactoriser::factorise() { size_t n = m_matrix->numRows(); - gpuErrChk(cusolverDnDpotrf(Session::getInstance().cuSolverHandle(), CUBLAS_FILL_MODE_LOWER, n, + gpuErrChk(cusolverDnDpotrf(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, n, m_matrix->raw(), n, m_workspace->raw(), m_workspaceSize, @@ -1753,7 +2059,7 @@ inline void CholeskyFactoriser::factorise() { template<> inline void CholeskyFactoriser::factorise() { size_t n = m_matrix->numRows(); - gpuErrChk(cusolverDnSpotrf(Session::getInstance().cuSolverHandle(), CUBLAS_FILL_MODE_LOWER, n, + gpuErrChk(cusolverDnSpotrf(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, n, m_matrix->raw(), n, m_workspace->raw(), m_workspaceSize, @@ -1763,7 +2069,7 @@ inline void CholeskyFactoriser::factorise() { template<> inline void CholeskyFactoriser::solve(DTensor &rhs) { size_t n = m_matrix->numRows(); - gpuErrChk(cusolverDnDpotrs(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnDpotrs(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, n, 1, m_matrix->raw(), n, @@ -1774,7 +2080,7 @@ inline void CholeskyFactoriser::solve(DTensor &rhs) { template<> inline void CholeskyFactoriser::solve(DTensor &rhs) { size_t n = m_matrix->numRows(); - gpuErrChk(cusolverDnSpotrs(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnSpotrs(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, n, 1, m_matrix->raw(), n, @@ -1797,7 +2103,7 @@ TEMPLATE_WITH_TYPE_T TEMPLATE_CONSTRAINT_REQUIRES_FPX class QRFactoriser : public IStatus { private: - int m_workspaceSize = 0; ///< Size of workspace needed for LS + int m_workspaceSize = 0; ///< Maximum workspace element count for geqrf and both ormqr calls std::unique_ptr > m_householder; ///< For storing householder reflectors std::unique_ptr > m_workspace; ///< Workspace for LS DTensor *m_matrix; ///< Lhs matrix template. Do not destroy! @@ -1809,16 +2115,19 @@ private: public: QRFactoriser(DTensor &A) : IStatus() { + if (!A.numRows() || !A.numCols() || !A.numMats()) throw std::invalid_argument("[QR] empty tensor"); if (A.numMats() > 1) throw std::invalid_argument("[QR] 3D tensors require `leastSquaresBatched`"); if (A.numRows() < A.numCols()) throw std::invalid_argument("[QR] Matrix A must be tall or square"); m_matrix = &A; + m_info->setStreamIdx(A.streamIdx()); + m_householder = std::make_unique >(m_matrix->numCols()); computeWorkspaceSize(); m_workspace = std::make_unique >(m_workspaceSize); - m_householder = std::make_unique >(m_matrix->numCols()); } /** * Factorise matrix. + * Work is enqueued on the bound matrix's selected stream. */ void factorise(); @@ -1827,6 +2136,8 @@ public: * A is the matrix that is factorised and b is the provided matrix. * A and b must have compatible dimensions (same number of rows and matrices=1). * A must be tall or square (m>=n). + * Both solver and BLAS work use A's stream, independent of b.streamIdx(). + * Keep inputs and this workspace alive until that stream completes. * @param b provided matrix */ void leastSquares(DTensor &b); @@ -1834,6 +2145,7 @@ public: /** * Populate the given tensors with Q and R. * Caution! This is an inefficient method: only to be used for debugging. + * It synchronizes A's stream before reading the triangular factor on host. * * @param Q matrix Q (preallocated) * @param R matrix R (preallocated) @@ -1847,27 +2159,37 @@ template<> inline void QRFactoriser::computeWorkspaceSize() { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); - gpuErrChk(cusolverDnDgeqrf_bufferSize(Session::getInstance().cuSolverHandle(), - m, n, - nullptr, m, - &m_workspaceSize)); + auto handle = Session::getInstance().cuSolverHandle(m_matrix->streamIdx()); + int factoriseSize = 0, solveSize = 0, extractSize = 0; + gpuErrChk(cusolverDnDgeqrf_bufferSize(handle, m, n, m_matrix->raw(), m, &factoriseSize)); + // Buffer-size queries do not modify C; use A as a sufficiently sized placeholder. + gpuErrChk(cusolverDnDormqr_bufferSize(handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_T, + m, 1, n, m_matrix->raw(), m, m_householder->raw(), m_matrix->raw(), m, &solveSize)); + gpuErrChk(cusolverDnDormqr_bufferSize(handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_N, + m, n, n, m_matrix->raw(), m, m_householder->raw(), m_matrix->raw(), m, &extractSize)); + m_workspaceSize = std::max({factoriseSize, solveSize, extractSize}); } template<> inline void QRFactoriser::computeWorkspaceSize() { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); - gpuErrChk(cusolverDnSgeqrf_bufferSize(Session::getInstance().cuSolverHandle(), - m, n, - nullptr, m, - &m_workspaceSize)); + auto handle = Session::getInstance().cuSolverHandle(m_matrix->streamIdx()); + int factoriseSize = 0, solveSize = 0, extractSize = 0; + gpuErrChk(cusolverDnSgeqrf_bufferSize(handle, m, n, m_matrix->raw(), m, &factoriseSize)); + // Buffer-size queries do not modify C; use A as a sufficiently sized placeholder. + gpuErrChk(cusolverDnSormqr_bufferSize(handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_T, + m, 1, n, m_matrix->raw(), m, m_householder->raw(), m_matrix->raw(), m, &solveSize)); + gpuErrChk(cusolverDnSormqr_bufferSize(handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_N, + m, n, n, m_matrix->raw(), m, m_householder->raw(), m_matrix->raw(), m, &extractSize)); + m_workspaceSize = std::max({factoriseSize, solveSize, extractSize}); } template<> inline void QRFactoriser::factorise() { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); - gpuErrChk(cusolverDnDgeqrf(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnDgeqrf(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), m, n, m_matrix->raw(), m, m_householder->raw(), @@ -1880,7 +2202,7 @@ template<> inline void QRFactoriser::factorise() { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); - gpuErrChk(cusolverDnSgeqrf(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnSgeqrf(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), m, n, m_matrix->raw(), m, m_householder->raw(), @@ -1893,7 +2215,7 @@ inline void QRFactoriser::leastSquares(DTensor &rhs) { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); double alpha = 1.; - gpuErrChk(cusolverDnDormqr(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnDormqr(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_SIDE_LEFT, CUBLAS_OP_T, m, 1, n, m_matrix->raw(), m, m_householder->raw(), @@ -1912,7 +2234,7 @@ inline void QRFactoriser::leastSquares(DTensor &rhs) { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); float alpha = 1.; - gpuErrChk(cusolverDnSormqr(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnSormqr(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_SIDE_LEFT, CUBLAS_OP_T, m, 1, n, m_matrix->raw(), m, m_householder->raw(), @@ -1930,68 +2252,78 @@ template<> inline void QRFactoriser::getQR(DTensor &Q, DTensor &R) { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); - if (Q.numRows() != m || Q.numCols() != n) + if (Q.numRows() != m || Q.numCols() != n || Q.numMats() != 1) throw std::invalid_argument("[QR] invalid shape of Q."); - if (R.numRows() != n || R.numCols() != n) + if (R.numRows() != n || R.numCols() != n || R.numMats() != 1) throw std::invalid_argument("[QR] invalid shape of R."); // Initialize Q to 1's on diagonal std::vector vecQ(m * n, 0.); for (size_t r = 0; r < m; r++) { for (size_t c = 0; c < n; c++) { - if (r == c) { vecQ[r * n + c] = 1.; } + if (r == c) { vecQ[r + c * m] = 1.; } } } - Q.upload(vecQ, rowMajor); + cudaStream_t stream; + gpuErrChk(cusolverDnGetStream(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), &stream)); + gpuErrChk(cudaMemcpyAsync(Q.raw(), vecQ.data(), vecQ.size() * sizeof(double), cudaMemcpyHostToDevice, stream)); // Apply Householder reflectors to compute Q - gpuErrChk(cusolverDnDormqr(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnDormqr(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_SIDE_LEFT, CUBLAS_OP_N, m, n, n, m_matrix->raw(), m, m_householder->raw(), Q.raw(), m, m_workspace->raw(), m_workspaceSize, m_info->raw())); + // This debugging helper reads R on the host; finish the selected stream first. + Session::getInstance().synchronizeStream(m_matrix->streamIdx()); // Extract upper triangular R std::vector vecR(n * n, 0.); for (size_t r = 0; r < n; r++) { for (size_t c = r; c < n; c++) { - vecR[r * n + c] = (*m_matrix)(r, c); + vecR[r + c * n] = (*m_matrix)(r, c); } } - R.upload(vecR, rowMajor); + gpuErrChk(cudaMemcpyAsync(R.raw(), vecR.data(), vecR.size() * sizeof(double), cudaMemcpyHostToDevice, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); // Complete host-buffer use before return. } template<> inline void QRFactoriser::getQR(DTensor &Q, DTensor &R) { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); - if (Q.numRows() != m || Q.numCols() != n) + if (Q.numRows() != m || Q.numCols() != n || Q.numMats() != 1) throw std::invalid_argument("[QR] invalid shape of Q."); - if (R.numRows() != n || R.numCols() != n) + if (R.numRows() != n || R.numCols() != n || R.numMats() != 1) throw std::invalid_argument("[QR] invalid shape of R."); // Initialize Q to 1's on diagonal std::vector vecQ(m * n, 0.); for (size_t r = 0; r < m; r++) { for (size_t c = 0; c < n; c++) { - if (r == c) { vecQ[r * n + c] = 1.; } + if (r == c) { vecQ[r + c * m] = 1.; } } } - Q.upload(vecQ, rowMajor); + cudaStream_t stream; + gpuErrChk(cusolverDnGetStream(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), &stream)); + gpuErrChk(cudaMemcpyAsync(Q.raw(), vecQ.data(), vecQ.size() * sizeof(float), cudaMemcpyHostToDevice, stream)); // Apply Householder reflectors to compute Q - gpuErrChk(cusolverDnSormqr(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnSormqr(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_SIDE_LEFT, CUBLAS_OP_N, m, n, n, m_matrix->raw(), m, m_householder->raw(), Q.raw(), m, m_workspace->raw(), m_workspaceSize, m_info->raw())); + // This debugging helper reads R on the host; finish the selected stream first. + Session::getInstance().synchronizeStream(m_matrix->streamIdx()); // Extract upper triangular R std::vector vecR(n * n, 0.); for (size_t r = 0; r < n; r++) { for (size_t c = r; c < n; c++) { - vecR[r * n + c] = (*m_matrix)(r, c); + vecR[r + c * n] = (*m_matrix)(r, c); } } - R.upload(vecR, rowMajor); + gpuErrChk(cudaMemcpyAsync(R.raw(), vecR.data(), vecR.size() * sizeof(float), cudaMemcpyHostToDevice, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); // Complete host-buffer use before return. } @@ -2011,6 +2343,7 @@ class Nullspace { private: std::unique_ptr > m_nullspace; ///< Stores all nullspace matrices (N) std::unique_ptr > m_projOp; ///< Stores all projection operators (N*N') + std::unique_ptr > m_projectionScratch; ///< Constructor-allocated vector output, separate from GEMM input public: /** @@ -2036,6 +2369,11 @@ public: * b = (b1, ..., bk) onto the nullspace of the original tensor. * That is, projection zi = Ni * Ni' * bi. * The projection zi is stored in bi. + * b must have shape (original numCols, 1, original numMats). + * GEMM and copy-back are queued on b's selected stream without allocations. + * Keep this object and b alive until that stream completes. Reuse on one + * stream is ordered; finish prior work before reusing this object on another. + * @throws std::invalid_argument if b has an incompatible shape */ void project(DTensor &b); }; @@ -2045,9 +2383,10 @@ template TEMPLATE_CONSTRAINT_REQUIRES_FPX inline Nullspace::Nullspace(DTensor &a) { size_t m = a.numRows(), n = a.numCols(), nMats = a.numMats(); - if (m > n) throw std::invalid_argument("[nullspace] I was expecting a square or fat matrix"); + if (!m || !n || !nMats || m > n) throw std::invalid_argument("[nullspace] expected a nonempty square or fat matrix"); m_nullspace = std::make_unique >(n, n, nMats, true); m_projOp = std::make_unique >(n, n, nMats, true); + m_projectionScratch = std::make_unique >(n, 1, nMats); auto aTranspose = a.tr(); Svd svd(aTranspose, true); svd.factorise(); @@ -2076,12 +2415,24 @@ inline Nullspace::Nullspace(DTensor &a) { DTensor Ntr = currNsDst.tr(); currProjOpDst.addAB(currNsDst, Ntr); } + // Setup is complete before project() can consume the projector on any stream. + cudaStream_t setupStream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(m_projOp->streamIdx()), &setupStream)); + gpuErrChk(cudaStreamSynchronize(setupStream)); } template TEMPLATE_CONSTRAINT_REQUIRES_FPX inline void Nullspace::project(DTensor &b) { - b.addAB(*m_projOp, b, 1, 0); + if (b.numRows() != m_projectionScratch->numRows() || b.numCols() != 1 + || b.numMats() != m_projectionScratch->numMats()) + throw std::invalid_argument("[nullspace] expected one vector per original matrix"); + m_projectionScratch->setStreamIdx(b.streamIdx()); + m_projectionScratch->addAB(*m_projOp, b, 1, 0); + cudaStream_t stream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(b.streamIdx()), &stream)); + gpuErrChk(cudaMemcpyAsync(b.raw(), m_projectionScratch->raw(), b.numEl() * sizeof(T), + cudaMemcpyDeviceToDevice, stream)); } @@ -2100,10 +2451,28 @@ TEMPLATE_CONSTRAINT_REQUIRES_FPX class CholeskyBatchFactoriser : public IStatus { private: DTensor *m_matrix; ///< Matrix to factorise or lower-triangular decomposed matrix + std::unique_ptr> m_singleMatrixPointers; ///< Two one-entry arrays, only for batches of one size_t m_numRows = 0; ///< Number of rows in `m_matrix` size_t m_numMats = 0; ///< Number of matrices in `m_matrix` bool m_factorisationDone = false; ///< Whether `m_matrix` holds original or lower-triangular matrix + void prepareSingleMatrixPointers(DTensor *rhs = nullptr) { + if (!m_singleMatrixPointers) return; + cudaStream_t stream; + gpuErrChk(cusolverDnGetStream(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), &stream)); + k_setSingleMatrixPointers<<<1, 1, 0, stream>>>(m_singleMatrixPointers->raw(), + m_matrix->raw(), rhs ? rhs->raw() : nullptr); + gpuErrChk(cudaGetLastError()); + } + + T **matrixPointers() const { + return m_singleMatrixPointers ? m_singleMatrixPointers->raw() : m_matrix->ptrMatrices(); + } + + T **rhsPointers(DTensor &rhs) const { + return m_singleMatrixPointers ? m_singleMatrixPointers->raw() + 1 : rhs.ptrMatrices(); + } + public: CholeskyBatchFactoriser() = delete; @@ -2114,10 +2483,12 @@ public: */ CholeskyBatchFactoriser(DTensor &A, bool factorised = false) : IStatus(A.numMats()), m_factorisationDone(factorised) { - if (A.numRows() != A.numCols()) throw std::invalid_argument("[CholeskyBatch] A must be square"); + if (!A.numRows() || !A.numMats() || A.numRows() != A.numCols()) throw std::invalid_argument("[CholeskyBatch] A must be nonempty and square"); m_matrix = &A; + m_info->setStreamIdx(A.streamIdx()); m_numRows = A.numRows(); m_numMats = A.numMats(); + if (m_numMats == 1) m_singleMatrixPointers = std::make_unique>(2); } /** @@ -2135,10 +2506,11 @@ public: template<> inline void CholeskyBatchFactoriser::factorise() { if (m_factorisationDone) return; - gpuErrChk(cusolverDnDpotrfBatched(Session::getInstance().cuSolverHandle(), + prepareSingleMatrixPointers(); + gpuErrChk(cusolverDnDpotrfBatched(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, m_numRows, - m_matrix->ptrMatrices(), + matrixPointers(), m_numRows, m_info->raw(), m_numMats)); @@ -2148,10 +2520,11 @@ inline void CholeskyBatchFactoriser::factorise() { template<> inline void CholeskyBatchFactoriser::factorise() { if (m_factorisationDone) return; - gpuErrChk(cusolverDnSpotrfBatched(Session::getInstance().cuSolverHandle(), + prepareSingleMatrixPointers(); + gpuErrChk(cusolverDnSpotrfBatched(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, m_numRows, - m_matrix->ptrMatrices(), + matrixPointers(), m_numRows, m_info->raw(), m_numMats)); @@ -2165,13 +2538,14 @@ inline void CholeskyBatchFactoriser::solve(DTensor &b) { throw std::invalid_argument("[CholeskyBatchSolve] A and b incompatible"); } if (b.numCols() != 1) throw std::invalid_argument("[CholeskyBatchSolve] only supports `b` with one column"); - gpuErrChk(cusolverDnDpotrsBatched(Session::getInstance().cuSolverHandle(), + prepareSingleMatrixPointers(&b); + gpuErrChk(cusolverDnDpotrsBatched(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, m_numRows, 1, ///< only supports rhs = 1 - m_matrix->ptrMatrices(), + matrixPointers(), m_numRows, - b.ptrMatrices(), + rhsPointers(b), m_numRows, m_info->raw(), m_numMats)); @@ -2184,13 +2558,14 @@ inline void CholeskyBatchFactoriser::solve(DTensor &b) { throw std::invalid_argument("[CholeskyBatchSolve] A and b incompatible"); } if (b.numCols() != 1) throw std::invalid_argument("[CholeskyBatchSolve] only supports `b` with one column"); - gpuErrChk(cusolverDnSpotrsBatched(Session::getInstance().cuSolverHandle(), + prepareSingleMatrixPointers(&b); + gpuErrChk(cusolverDnSpotrsBatched(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, m_numRows, 1, ///< only supports rhs = 1 - m_matrix->ptrMatrices(), + matrixPointers(), m_numRows, - b.ptrMatrices(), + rhsPointers(b), m_numRows, m_info->raw(), m_numMats)); @@ -2234,7 +2609,7 @@ public: * @throws std::invalid_argument if a.numMats() > 1 */ GivensAnnihilator(DTensor &a) { - if (a.numMats() > 1) { + if (!a.numRows() || !a.numCols() || a.numMats() != 1) { throw std::invalid_argument("[GivensAnnihilator] tensors (numMats > 1) not supported"); } m_matrix = &a; @@ -2248,7 +2623,7 @@ public: * @param a */ void setMatrix(DTensor &a) { - if (a.numMats() > 1) { + if (!a.numRows() || !a.numCols() || a.numMats() != 1) { throw std::invalid_argument("[GivensAnnihilator] tensors (numMats > 1) not supported"); } m_matrix = &a; @@ -2275,6 +2650,13 @@ __global__ void k_givensAnnihilateRHypot(const T *data, size_t nRows) { T xij = data[i + j * nRows]; T xkj = data[k + j * nRows]; + if (xij == T(0) && xkj == T(0)) { + // The target is already zero: leave both rows unchanged. Avoid 0 * infinity. + res[0] = T(0); + res[1] = T(1); // cos + res[2] = T(0); // -sin + return; + } res[0] = rhypot(xij, xkj); res[1] = xij * (*res); // cos res[2] = xkj * (*res); // -sin @@ -2293,22 +2675,25 @@ inline void GivensAnnihilator::annihilate(size_t i, size_t k, size_t j) { * Pass cosine and sine as device pointers * (Avoid having to download first) */ - gpuErrChk( - cublasSetPointerMode(Session::getInstance().cuBlasHandle(m_matrix->streamIdx()), CUBLAS_POINTER_MODE_DEVICE)); + auto handle = Session::getInstance().cuBlasHandle(m_matrix->streamIdx()); + cudaStream_t stream; + cublasPointerMode_t previousMode; + gpuErrChk(cublasGetStream(handle, &stream)); + gpuErrChk(cublasGetPointerMode(handle, &previousMode)); + gpuErrChk(cublasSetPointerMode(handle, CUBLAS_POINTER_MODE_DEVICE)); /* Useful definitions */ T *aux = m_d_rhyp_cos_sin->raw(); - T *matData = m_matrix->raw(); /* Call kernel to determine 1/sqrt(Ai^2 + Ak^2) */ - k_givensAnnihilateRHypot<<<1, 1>>>(m_matrix->raw(), aux, i, k, j, nR); + k_givensAnnihilateRHypot<<<1, 1, 0, stream>>>(m_matrix->raw(), aux, i, k, j, nR); + gpuErrChk(cudaGetLastError()); /* Apply Givens rotation */ m_matrix->applyLeftGivensRotation(i, k, aux + 1, aux + 2); - /* Change back to default behaviour */ - gpuErrChk( - cublasSetPointerMode(Session::getInstance().cuBlasHandle(m_matrix->streamIdx()), CUBLAS_POINTER_MODE_HOST)); + /* Restore the caller's pointer mode. */ + gpuErrChk(cublasSetPointerMode(handle, previousMode)); } diff --git a/python/README.md b/python/README.md index ebde65a..e41320c 100644 --- a/python/README.md +++ b/python/README.md @@ -25,4 +25,54 @@ import gputils_api as g x = g.read_array_from_gputils_binary_file('my_data.bt') ``` +### Tests + +From the repository root, run all Python tests with: + +```bash +./scripts/run_python_tests.sh +``` + +The script uses `build/test-venv/bin/python` if available, otherwise `python3`. +Set `PYTHON` to choose an interpreter with NumPy and pytest installed. It honours +`GPUTILS_BUILD_DIR` and `GPUTILS_FILE_EXCHANGE` and accepts extra pytest arguments, +for example `./scripts/run_python_tests.sh -k round_trip`. + +To invoke pytest directly: + +```bash +PYTHONPATH=python python -m pytest python/test -v +``` + +Tests use temporary directories, deterministic values, and full-array comparisons. +They cover shapes and precisions, non-contiguous inputs, independent binary-layout +fixtures, malformed files, and filesystem errors. Empty and zero-dimensional arrays +are expected to raise `ValueError`; one-element vectors remain valid inputs. + +The reader validates the complete 24-byte header, positive dimensions, NumPy size +limits, and exact payload byte count before loading array data. It rejects both +truncated payloads and trailing bytes with `ValueError`. The writer validates inputs +before opening the output file, so invalid inputs do not create or overwrite files. +Object-containing and zero-size dtypes are rejected by both functions. + +Actual C++/Python file-exchange tests need the test-only `gputils_file_exchange` +executable. Build with `GPUTILS_BUILD_TEST=ON`; the tests find it automatically at +`build/test/gputils_file_exchange`. For another build directory, set: + +```bash +GPUTILS_FILE_EXCHANGE=/absolute/path/to/build/test/gputils_file_exchange \ + PYTHONPATH=python python -m pytest python/test -v +``` + +Without a built helper, only the C++ exchange tests are skipped. An explicitly +configured missing helper is an error. `scripts/build_and_test.sh` builds the helper +and sets its path, so these tests run as part of the combined suite. + +The format stores three little-endian `uint64` dimensions, followed by consecutive +column-major matrices. It does not store a dtype: pass the matching `dt` to the +reader. Same-width dtype mistakes cannot reliably be detected. Python can preserve +explicit payload byte orders when supplied the matching dtype; C++ interoperability +tests use little-endian float32/float64 data on a little-endian host. Read results +always have three axes, padding vector and matrix shapes with singleton dimensions. + diff --git a/python/gputils_api/gputils_api.py b/python/gputils_api/gputils_api.py index 3fef87f..e291fe3 100644 --- a/python/gputils_api/gputils_api.py +++ b/python/gputils_api/gputils_api.py @@ -1,24 +1,48 @@ +import os +import struct + import numpy as np + +def _fixed_size_dtype(dt): + dtype = np.dtype(dt) + if dtype.hasobject or dtype.itemsize == 0: + raise ValueError("The binary format requires a fixed-size dtype without objects") + return dtype + + def read_array_from_gputils_binary_file(path, dt=np.dtype('d')): """ Reads an array from a bt file :param path: path to file :param dt: numpy-compatible data type :raises ValueError: if the file name specified `path` does not have the .bt extension + :raises ValueError: if the header, dimensions, dtype, or payload size are invalid """ if not path.endswith(".bt"): raise ValueError("The file must have the .bt extension") + dtype = _fixed_size_dtype(dt) with open(path, 'rb') as f: - nr = int.from_bytes(f.read(8), byteorder='little', signed=False) # read number of rows - nc = int.from_bytes(f.read(8), byteorder='little', signed=False) # read number of columns - nm = int.from_bytes(f.read(8), byteorder='little', signed=False) # read number of matrices - dat = np.fromfile(f, dtype=np.dtype(dt)) # read data - - if nm >= 2: # if we actually have a 3D tensor (not a matrix or a vector) - dat = dat.reshape((nm, nc, nr)).swapaxes(0, 2) # I'll explain this to you when you grow up - else: - dat = dat.reshape((nr, nc, nm)) # reshape + header = f.read(24) + if len(header) != 24: + raise ValueError("The binary tensor header must contain exactly 24 bytes") + # Now decode the 24-byte binary header into three integers... + nr, nc, nm = struct.unpack(' np.iinfo(np.intp).max: + raise ValueError("Tensor dimensions exceed NumPy size limits") + if os.fstat(f.fileno()).st_size - 24 != payload_size: + raise ValueError("The binary payload size does not match its dimensions and dtype") + + dat = np.fromfile(f, dtype=dtype, count=count) + if dat.size != count or f.read(1): + raise ValueError("The binary payload changed size while being read") + + # Each matrix is column-major, with matrices stored consecutively. + dat = dat.reshape((nm, nc, nr)).swapaxes(0, 2) return dat @@ -28,7 +52,8 @@ def write_array_to_gputils_binary_file(x, path): :param x: numpy array to save to file :param path: path to file - :raises ValueError: if `x` has more than 3 dimensions + :raises ValueError: if `x` is empty or does not have between 1 and 3 dimensions + :raises ValueError: if `x` has a zero-size or object-containing dtype :raises ValueError: if the file name specified `path` does not have the .bt extension """ @@ -36,8 +61,11 @@ def write_array_to_gputils_binary_file(x, path): raise ValueError("The file must have the .bt extension") x_shape = x.shape x_dims = len(x_shape) - if x_dims >= 4: - raise ValueError("given array cannot have more than 3 dimensions") + if not 1 <= x_dims <= 3: + raise ValueError("The array must have between 1 and 3 dimensions") + if x.size == 0: + raise ValueError("Tensor dimensions must be positive") + _fixed_size_dtype(x.dtype) nr = x_shape[0] nc = x_shape[1] if x_dims >= 2 else 1 nm = x_shape[2] if x_dims == 3 else 1 @@ -49,4 +77,4 @@ def write_array_to_gputils_binary_file(x, path): f.write(nr.to_bytes(8, 'little')) # write number of rows f.write(nc.to_bytes(8, 'little')) # write number of columns f.write(nm.to_bytes(8, 'little')) # write number of matrices - x.tofile(f) # write data \ No newline at end of file + x.tofile(f) # write data diff --git a/python/pytest.ini b/python/pytest.ini new file mode 100644 index 0000000..aa1640f --- /dev/null +++ b/python/pytest.ini @@ -0,0 +1,4 @@ +[pytest] +python_files = test.py test_*.py *_test.py +testpaths = test +pythonpath = . diff --git a/python/test/test.py b/python/test/test.py index 5364191..eef12ee 100644 --- a/python/test/test.py +++ b/python/test/test.py @@ -1,83 +1,303 @@ import os +from pathlib import Path +import struct +import subprocess +import sys +import tempfile import unittest +from unittest import mock + import numpy as np import gputils_api as gpuapi -class GputilApiTestCase(unittest.TestCase): - - _B = np.array([ - [[1, 2], [3, 4], [5, 6]], - [[7, 8], [9, 10], [11, 12]], - [[13, 14], [15, 16], [17, 18]] - ], dtype=np.dtype('d')) - - @staticmethod - def local_abs_path(): - cwd = os.getcwd() - return cwd.split('open-codegen')[0] - - @classmethod - def setUpClass(cls): - n = 5 - base_dir = GputilApiTestCase.local_abs_path() - eye_d = np.eye(n, dtype=np.dtype('d')) - gpuapi.write_array_to_gputils_binary_file(eye_d, os.path.join(base_dir, 'eye_d.bt')) - - eye_f = np.eye(n, dtype=np.dtype('f')) - gpuapi.write_array_to_gputils_binary_file(eye_f, os.path.join(base_dir, 'eye_f.bt')) - - xd = np.random.randn(2, 4, 6).astype('d') - xd[1, 2, 3] = -12.3 - gpuapi.write_array_to_gputils_binary_file(xd, os.path.join(base_dir, 'rand_246_d.bt')) +DTYPES = (np.float32, np.float64) - xf = np.random.randn(2, 4, 6).astype('f') - xf[1, 2, 3] = float(-12.3) - gpuapi.write_array_to_gputils_binary_file(xf, os.path.join(base_dir, 'rand_246_f.bt')) - a = np.linspace(-100, 100, 4 * 5).reshape((4, 5)).astype('d') - gpuapi.write_array_to_gputils_binary_file(a, os.path.join(base_dir, 'a_d.bt')) +class FileTestCase(unittest.TestCase): + def setUp(self): + directory = tempfile.TemporaryDirectory(prefix='gputils-io-') + self.addCleanup(directory.cleanup) + self.directory = Path(directory.name) + self.path = str(self.directory / 'tensor.bt') - gpuapi.write_array_to_gputils_binary_file(cls._B, os.path.join(base_dir, 'b_d.bt')) + def assert_round_trip(self, original): + gpuapi.write_array_to_gputils_binary_file(original, self.path) + actual = gpuapi.read_array_from_gputils_binary_file(self.path, original.dtype) + expected_shape = original.shape + (1,) * (3 - original.ndim) + self.assertEqual(actual.shape, expected_shape) + self.assertEqual(actual.dtype, original.dtype) + np.testing.assert_array_equal(actual, original.reshape(expected_shape)) + return actual - def __test_read_eye(self, dt): - base_dir = GputilApiTestCase.local_abs_path() - path = os.path.join(base_dir, f'eye_{dt}.bt') - r = gpuapi.read_array_from_gputils_binary_file(path, dt=np.dtype(dt)) - err = r[:, :, 0] - np.eye(5) - err_norm = np.linalg.norm(err, np.inf) - self.assertTrue(err_norm < 1e-12) +class GputilApiTestCase(FileTestCase): def test_read_eye_d(self): - self.__test_read_eye('d') + self.assert_round_trip(np.eye(5, dtype=np.float64)) def test_read_eye_f(self): - self.__test_read_eye('f') - - def __test_read_rand(self, dt): - base_dir = GputilApiTestCase.local_abs_path() - path = os.path.join(base_dir, f'rand_246_{dt}.bt') - r = gpuapi.read_array_from_gputils_binary_file(path, dt=np.dtype(dt)) - r_shape = r.shape - self.assertEqual(2, r_shape[0]) - self.assertEqual(4, r_shape[1]) - self.assertEqual(6, r_shape[2]) - e = np.abs(r[1, 2, 3] + 12.3) - self.assertTrue(e < 1e-6) + self.assert_round_trip(np.eye(5, dtype=np.float32)) + + def test_single_matrix_round_trip(self): + self.assert_round_trip(np.array([[1., 2., 3.], [4., 5., 6.]], dtype=np.float64)) def test_read_rand_d(self): - self.__test_read_rand('d') + # Deterministic values and a full comparison replace random spot checks. + self.assert_round_trip(np.arange(48, dtype=np.float64).reshape(2, 4, 6) / 4 - 5) def test_read_rand_f(self): - self.__test_read_rand('f') + self.assert_round_trip(np.arange(48, dtype=np.float32).reshape(2, 4, 6) / 4 - 5) def test_read_a_tensor_3d(self): - self.__test_read_rand('f') - base_dir = GputilApiTestCase.local_abs_path() - path = os.path.join(base_dir, 'b_d.bt') - b = gpuapi.read_array_from_gputils_binary_file(path) - self.assertEqual(GputilApiTestCase._B.shape, b.shape) - self.assertLess(np.linalg.norm((b - GputilApiTestCase._B).reshape(-1, 1)), 1e-6) + self.assert_round_trip(np.arange(18, dtype=np.float64).reshape(3, 3, 2) - 9) + + def test_round_trip_shapes_and_precisions(self): + for dtype in DTYPES: + for shape in ((1,), (5,), (1, 1), (1, 5), (5, 1), (2, 3), (2, 3, 1), (2, 3, 2)): + with self.subTest(dtype=dtype, shape=shape): + data = np.arange(np.prod(shape), dtype=dtype).reshape(shape) / 4 - 3 + self.assert_round_trip(data) + + def test_noncontiguous_and_fortran_inputs(self): + for dtype in DTYPES: + matrix = np.arange(48, dtype=dtype).reshape(6, 8) + tensor = np.arange(120, dtype=dtype).reshape(4, 6, 5) + inputs = (matrix.T, matrix[::2, 1::2], matrix[::-1, ::-1], + tensor[::2, 1::2, ::-1], np.asfortranarray(matrix)) + for index, data in enumerate(inputs): + with self.subTest(dtype=dtype, layout=index): + self.assert_round_trip(data) + + def test_writer_matches_independent_binary_layout(self): + for dtype, code in ((np.float32, 'f'), (np.float64, 'd')): + for shape in ((2, 3, 1), (2, 3, 2)): + with self.subTest(dtype=dtype, shape=shape): + data = np.arange(np.prod(shape), dtype=dtype).reshape(shape) / 4 - 2 + expected = struct.pack('f4', 'f8', '/dev/null 2>&1; then + echo "Python with NumPy is required. Install python/ in a virtual environment and set PYTHON to its interpreter." >&2 + exit 1 +fi + +# Extra arguments are forwarded to CMake, e.g. -DSM_ARCH=120 -DCPPVERSION=20. +cmake -S "$repo_root" -B "$build_dir" -DCMAKE_BUILD_TYPE=RelWithDebInfo "$@" -DGPUTILS_BUILD_TEST=ON +cmake --build "$build_dir" --parallel "$jobs" +ctest --test-dir "$build_dir/test" --output-on-failure --no-tests=error + +# The Python tests write fixtures to their working directory. Keep those in +# the build tree rather than the source checkout. +mkdir -p "$build_dir/python-tests" +( + cd -- "$build_dir/python-tests" + PYTHONPATH="$repo_root/python${PYTHONPATH:+:$PYTHONPATH}" \ + GPUTILS_FILE_EXCHANGE="$build_dir/test/gputils_file_exchange" \ + "$python" "$repo_root/python/test/test.py" -v +) diff --git a/scripts/check_memory.sh b/scripts/check_memory.sh new file mode 100755 index 0000000..bb6b0d9 --- /dev/null +++ b/scripts/check_memory.sh @@ -0,0 +1,45 @@ +#!/usr/bin/env bash +# Build/test, then check the full CUDA suite and main executable for leaks. +set -euo pipefail + +repo_root="$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")/.." && pwd)" +build_dir="${GPUTILS_BUILD_DIR:-$repo_root/build}" +case "$build_dir" in /*) ;; *) build_dir="$repo_root/$build_dir" ;; esac +sanitizer="${COMPUTE_SANITIZER:-/usr/local/cuda/bin/compute-sanitizer}" +repeats="${SANITIZER_REPEATS:-3}" +if ! command -v "$sanitizer" >/dev/null 2>&1; then + echo "Compute Sanitizer not found; set COMPUTE_SANITIZER to its path." >&2 + exit 1 +fi +if [[ ! "$repeats" =~ ^[1-9][0-9]*$ ]]; then + echo "SANITIZER_REPEATS must be a positive integer." >&2 + exit 1 +fi + +bash "$repo_root/scripts/build_and_test.sh" \ + '-DCMAKE_CUDA_FLAGS_RELWITHDEBINFO=-O2 -g -lineinfo -DNDEBUG' "$@" +mkdir -p "$build_dir/sanitizer" +cd -- "$build_dir/sanitizer" +options=(--tool memcheck --leak-check full --check-device-heap yes --error-exitcode 1) + +# One process, shuffled repetitions: catches leaks and stale ownership retained +# across tests/repeated use as well as invalid device accesses. pipefail preserves +# sanitizer or test failures while tee keeps complete logs. +"$sanitizer" "${options[@]}" "$build_dir/test/gputils_test" \ + --gtest_shuffle --gtest_random_seed=1 --gtest_repeat="$repeats" \ + 2>&1 | tee tests-memcheck.log +"$sanitizer" "${options[@]}" "$build_dir/gputils_main" \ + 2>&1 | tee main-memcheck.log + +# Fail closed if instrumentation did not produce a successful completion report. +for report in tests-memcheck.log main-memcheck.log; do + if ! grep -Eq 'ERROR SUMMARY: 0 errors' "$report"; then + echo "Missing clean sanitizer summary: $report" >&2 + exit 1 + fi + if ! grep -Eq 'LEAK SUMMARY: 0 bytes leaked in 0 allocations' "$report"; then + echo "Leak check failed: $report" >&2 + exit 1 + fi +done +echo "CUDA memory checks passed; logs: $build_dir/sanitizer" diff --git a/scripts/run_python_tests.sh b/scripts/run_python_tests.sh new file mode 100755 index 0000000..201a352 --- /dev/null +++ b/scripts/run_python_tests.sh @@ -0,0 +1,24 @@ +#!/usr/bin/env bash +# Run all Python tests with pytest; forward any extra pytest arguments. +set -euo pipefail + +repo_root="$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")/.." && pwd)" +build_dir="${GPUTILS_BUILD_DIR:-$repo_root/build}" +case "$build_dir" in /*) ;; *) build_dir="$repo_root/$build_dir" ;; esac +python="${PYTHON:-python3}" +if [[ -z "${PYTHON:-}" && -x "$build_dir/test-venv/bin/python" ]]; then + python="$build_dir/test-venv/bin/python" +fi + +if ! "$python" -c 'import numpy, pytest' >/dev/null 2>&1; then + echo "Python with NumPy and pytest is required. Install them in a virtual environment and set PYTHON to its interpreter." >&2 + exit 1 +fi + +# Reuse a built helper for the C++ exchange tests; honour an explicit override. +# The Python suite reports those tests as skipped if no helper is available. +if [[ -z "${GPUTILS_FILE_EXCHANGE:-}" && -x "$build_dir/test/gputils_file_exchange" ]]; then + export GPUTILS_FILE_EXCHANGE="$build_dir/test/gputils_file_exchange" +fi +export PYTHONPATH="$repo_root/python${PYTHONPATH:+:$PYTHONPATH}" +exec "$python" -m pytest "$repo_root/python/test" -v "$@" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4abfa1c..8cf248c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -13,6 +13,9 @@ add_executable(gputils_test) target_sources(gputils_test # add files PRIVATE testTensor.cu + testFileIO.cu + hostAllocationFailure.cpp + sessionAcrossTranslationUnits.cu ) target_link_libraries(gputils_test PRIVATE @@ -28,4 +31,8 @@ target_include_directories(gputils_test ) include(GoogleTest) gtest_discover_tests(gputils_test) + +add_executable(gputils_file_exchange exchangeTensorFiles.cu) +target_link_libraries(gputils_file_exchange PRIVATE gputils_compiler_flags cublas cusolver cudadevrt) +target_include_directories(gputils_file_exchange PRIVATE "${PROJECT_SOURCE_DIR}/include") # ---- diff --git a/test/exchangeTensorFiles.cu b/test/exchangeTensorFiles.cu new file mode 100644 index 0000000..0d51cca --- /dev/null +++ b/test/exchangeTensorFiles.cu @@ -0,0 +1,45 @@ +// Test-only bridge: exercise the real C++ reader/writer from Python tests. +#include "tensor.cuh" +#include +#include +#include +#include + +template +void exchangeFile(const std::string &mode, const std::string &path, size_t rows, size_t cols, size_t mats) { + std::vector expected(rows * cols * mats); + for (size_t k = 0; k < mats; ++k) + for (size_t c = 0; c < cols; ++c) + for (size_t r = 0; r < rows; ++r) + expected[r + rows * (c + cols * k)] = T(r + 10 * c + 100 * k) + T(.25); + if (mode == "write") { + DTensor tensor(expected, rows, cols, mats); + tensor.saveToFile(path); + } else if (mode == "verify") { + auto tensor = DTensor::parseFromFile(path); + if (tensor.numRows() != rows || tensor.numCols() != cols || tensor.numMats() != mats) + throw std::runtime_error("C++ reader returned an incorrect shape"); + std::vector actual; + tensor.download(actual); + if (actual != expected) throw std::runtime_error("C++ reader returned incorrect element values"); + } else { + throw std::invalid_argument("mode must be write or verify"); + } +} + +int main(int argc, char **argv) { + try { + if (argc != 7) throw std::invalid_argument("usage: gputils_file_exchange write|verify float|double path rows cols mats"); + const size_t rows = std::stoull(argv[4]), cols = std::stoull(argv[5]), mats = std::stoull(argv[6]); + // This executable serves tiny deterministic test fixtures only. + if (!rows || !cols || !mats || rows > 16 || cols > 16 || mats > 16) + throw std::invalid_argument("test fixture dimensions must be in [1,16]"); + if (std::string(argv[2]) == "float") exchangeFile(argv[1], argv[3], rows, cols, mats); + else if (std::string(argv[2]) == "double") exchangeFile(argv[1], argv[3], rows, cols, mats); + else throw std::invalid_argument("dtype must be float or double"); + return 0; + } catch (const std::exception &error) { + std::cerr << error.what() << '\n'; + return 1; + } +} diff --git a/test/hostAllocationFailure.cpp b/test/hostAllocationFailure.cpp new file mode 100644 index 0000000..89ec1e5 --- /dev/null +++ b/test/hostAllocationFailure.cpp @@ -0,0 +1,36 @@ +#include "hostAllocationFailure.h" +#include +#include + +namespace { +thread_local gputils_test::AllocationFailurePredicate failurePredicate = nullptr; +} + +gputils_test::FailHostAllocation::FailHostAllocation(AllocationFailurePredicate predicate) noexcept + : previous(failurePredicate) { + failurePredicate = predicate; +} + +gputils_test::FailHostAllocation::~FailHostAllocation() { + failurePredicate = previous; +} + +void *operator new(std::size_t bytes) { + if (failurePredicate && failurePredicate(bytes)) { + // Disarm before throwing so cleanup and exception reporting can allocate. + failurePredicate = nullptr; + throw std::bad_alloc(); + } + for (;;) { + if (void *memory = std::malloc(bytes == 0 ? 1 : bytes)) return memory; + const auto handler = std::get_new_handler(); + if (!handler) throw std::bad_alloc(); + handler(); + } +} + +void *operator new[](std::size_t bytes) { return ::operator new(bytes); } +void operator delete(void *memory) noexcept { std::free(memory); } +void operator delete[](void *memory) noexcept { std::free(memory); } +void operator delete(void *memory, std::size_t) noexcept { std::free(memory); } +void operator delete[](void *memory, std::size_t) noexcept { std::free(memory); } diff --git a/test/hostAllocationFailure.h b/test/hostAllocationFailure.h new file mode 100644 index 0000000..780750c --- /dev/null +++ b/test/hostAllocationFailure.h @@ -0,0 +1,22 @@ +#ifndef GPUTILS_TEST_HOST_ALLOCATION_FAILURE_H +#define GPUTILS_TEST_HOST_ALLOCATION_FAILURE_H + +#include + +namespace gputils_test { +using AllocationFailurePredicate = bool (*)(std::size_t) noexcept; + +// Test-executable-only, thread-local, one-shot interception of host operator new. +// The predicate must not allocate. Production targets do not link this helper. +class FailHostAllocation { +public: + explicit FailHostAllocation(AllocationFailurePredicate predicate) noexcept; + ~FailHostAllocation(); + FailHostAllocation(const FailHostAllocation &) = delete; + FailHostAllocation &operator=(const FailHostAllocation &) = delete; +private: + AllocationFailurePredicate previous; +}; +} + +#endif diff --git a/test/sessionAcrossTranslationUnits.cu b/test/sessionAcrossTranslationUnits.cu new file mode 100644 index 0000000..64e4d74 --- /dev/null +++ b/test/sessionAcrossTranslationUnits.cu @@ -0,0 +1,19 @@ +#include "sessionAcrossTranslationUnits.h" +#include "tensor.cuh" + +void gputils_test::configureStreamsFromOtherTranslationUnit(std::size_t count) { + Session::setStreams(count); +} + +Session &gputils_test::sessionFromOtherTranslationUnit() { + return Session::getInstance(); +} + +bool gputils_test::executeOnSecondStreamFromOtherTranslationUnit() { + DTensor value(std::vector{3,4},2); + value.setStreamIdx(1); + value *= 2; + std::vector actual; + value.download(actual); + return actual == std::vector{6,8}; +} diff --git a/test/sessionAcrossTranslationUnits.h b/test/sessionAcrossTranslationUnits.h new file mode 100644 index 0000000..9ea9b0c --- /dev/null +++ b/test/sessionAcrossTranslationUnits.h @@ -0,0 +1,13 @@ +#ifndef GPUTILS_TEST_SESSION_ACROSS_TRANSLATION_UNITS_H +#define GPUTILS_TEST_SESSION_ACROSS_TRANSLATION_UNITS_H + +#include +class Session; + +namespace gputils_test { +void configureStreamsFromOtherTranslationUnit(std::size_t count); +Session &sessionFromOtherTranslationUnit(); +bool executeOnSecondStreamFromOtherTranslationUnit(); +} + +#endif diff --git a/test/testFileIO.cu b/test/testFileIO.cu new file mode 100644 index 0000000..f7b7af0 --- /dev/null +++ b/test/testFileIO.cu @@ -0,0 +1,192 @@ +#include +#include "tensor.cuh" +#include +#include +#include +#include +#include +#include +#include + +namespace { +class TemporaryIODirectory { + std::filesystem::path m_path; +public: + TemporaryIODirectory() { + static std::atomic counter{0}; + const auto stamp = std::chrono::steady_clock::now().time_since_epoch().count(); + for (int attempt = 0; attempt < 100; ++attempt) { + auto path = std::filesystem::temp_directory_path() + / ("gputils-io-" + std::to_string(stamp) + "-" + std::to_string(counter++)); + if (std::filesystem::create_directory(path)) { m_path = path; return; } + } + throw std::runtime_error("cannot create unique test directory"); + } + ~TemporaryIODirectory() { + std::error_code error; + std::filesystem::remove_all(m_path, error); + } + std::filesystem::path file(const std::string &name) const { return m_path / name; } +}; + +template +void appendNative(std::vector &bytes, T value) { + const auto *begin = reinterpret_cast(&value); + bytes.insert(bytes.end(), begin, begin + sizeof(T)); +} + +void writeBytes(const std::filesystem::path &path, const std::vector &bytes) { + std::ofstream output(path, std::ios::binary); + ASSERT_TRUE(output.is_open()); + output.write(bytes.data(), static_cast(bytes.size())); + output.close(); + ASSERT_TRUE(output.good()); +} + +template +void roundTripShapes() { + TemporaryIODirectory directory; + const std::array, 5> shapes{{{1,1,1}, {5,1,1}, {1,5,1}, {2,3,1}, {2,3,2}}}; + for (const auto &shape : shapes) { + std::vector expected(shape[0] * shape[1] * shape[2]); + for (size_t i = 0; i < expected.size(); ++i) expected[i] = T(int(i) - 5); + DTensor original(expected, shape[0], shape[1], shape[2]); + for (const auto *name : {"tensor.bt", "tensor.dtensor"}) { + SCOPED_TRACE(name); + const auto path = directory.file(name).string(); + original.saveToFile(path); + auto actual = DTensor::parseFromFile(path); + EXPECT_EQ(actual.numRows(), shape[0]); + EXPECT_EQ(actual.numCols(), shape[1]); + EXPECT_EQ(actual.numMats(), shape[2]); + std::vector values; + actual.download(values); + EXPECT_EQ(values, expected); + } + } +} + +template +void textPreservesPrecision() { + TemporaryIODirectory directory; + const auto path = directory.file("precision.dtensor").string(); + const std::vector expected{T(0), T(-1), T(1) / T(3), std::nextafter(T(1), T(2)), + std::numeric_limits::min(), -std::numeric_limits::min(), + std::numeric_limits::max(), -std::numeric_limits::max()}; + DTensor tensor(expected, 2, 4); + tensor.saveToFile(path); + auto actual = DTensor::parseFromFile(path); + std::vector values; + actual.download(values); + EXPECT_EQ(values, expected); +} + +template +void binaryMatchesIndependentLayout() { + TemporaryIODirectory directory; + const auto path = directory.file("layout.bt"); + std::vector expected; + appendNative(expected, 2); + appendNative(expected, 3); + appendNative(expected, 2); + // A pair of matrices [[1,2,3],[4,5,6]] and its negative, stored column-major. + const std::vector values{1,4,2,5,3,6,-1,-4,-2,-5,-3,-6}; + for (T value : values) appendNative(expected, value); + DTensor tensor(values, 2, 3, 2); + tensor.saveToFile(path.string()); + std::ifstream input(path, std::ios::binary); + const std::vector actual{std::istreambuf_iterator(input), std::istreambuf_iterator()}; + EXPECT_EQ(actual, expected); + // The reader is also tested directly against bytes independent of the writer. + writeBytes(path, expected); + auto parsed = DTensor::parseFromFile(path.string()); + EXPECT_EQ(parsed(0, 1, 0), T(2)); + EXPECT_EQ(parsed(1, 2, 1), T(-6)); + std::vector downloaded; + parsed.download(downloaded); + EXPECT_EQ(downloaded, values); +} +} // namespace + +TEST(FileIOTest, shapesRoundTripFloat) { roundTripShapes(); } +TEST(FileIOTest, shapesRoundTripDouble) { roundTripShapes(); } +TEST(FileIOTest, shapesRoundTripInt) { roundTripShapes(); } +TEST(FileIOTest, textPreservesFloatPrecision) { textPreservesPrecision(); } +TEST(FileIOTest, textPreservesDoublePrecision) { textPreservesPrecision(); } +TEST(FileIOTest, binaryLayoutFloat) { binaryMatchesIndependentLayout(); } +TEST(FileIOTest, binaryLayoutDouble) { binaryMatchesIndependentLayout(); } + +TEST(FileIOTest, textPreservesIntegerLimits) { + TemporaryIODirectory directory; + const auto path = directory.file("integers.dtensor").string(); + const std::vector expected{std::numeric_limits::min(), -12345, 0, std::numeric_limits::max()}; + DTensor tensor(expected, 2, 2); + tensor.saveToFile(path); + auto parsed = DTensor::parseFromFile(path); + std::vector actual; + parsed.download(actual); + EXPECT_EQ(actual, expected); +} + +TEST(FileIOTest, rejectsEveryTruncatedBinaryHeaderLength) { + TemporaryIODirectory directory; + const auto path = directory.file("header.bt"); + std::vector header; + for (uint64_t dimension : {2, 3, 1}) appendNative(header, dimension); + for (size_t length = 0; length < header.size(); ++length) { + SCOPED_TRACE(length); + writeBytes(path, std::vector(header.begin(), header.begin() + length)); + EXPECT_THROW(DTensor::parseFromFile(path.string()), std::invalid_argument); + } +} + +TEST(FileIOTest, rejectsInvalidBinaryDimensionsAndExtraBytes) { + TemporaryIODirectory directory; + const auto path = directory.file("invalid.bt"); + const std::array, 5> shapes{{{0,1,1}, {2,0,1}, {2,3,0}, + {uint64_t(INT_MAX) + 1,1,1}, {UINT64_MAX,UINT64_MAX,2}}}; + const size_t before = Session::getInstance().totalAllocatedBytes(); + for (const auto &shape : shapes) { + std::vector bytes; + for (auto dimension : shape) appendNative(bytes, dimension); + writeBytes(path, bytes); + EXPECT_THROW(DTensor::parseFromFile(path.string()), std::invalid_argument); + } + std::vector bytes; + for (uint64_t dimension : {2, 3, 1}) appendNative(bytes, dimension); + for (double value : {1.,4.,2.,5.,3.,6.}) appendNative(bytes, value); + bytes.push_back('\xff'); + writeBytes(path, bytes); + EXPECT_THROW(DTensor::parseFromFile(path.string()), std::invalid_argument); + EXPECT_EQ(Session::getInstance().totalAllocatedBytes(), before); +} + +TEST(FileIOTest, reportsOutputPathFailures) { + TemporaryIODirectory directory; + DTensor tensor(std::vector{1,2,3,4}, 2, 2); + for (const auto *name : {"tensor.bt", "tensor.dtensor"}) { + const auto missingParent = directory.file("missing") / name; + EXPECT_THROW(tensor.saveToFile(missingParent.string()), std::runtime_error); + const auto folder = directory.file(name); + ASSERT_TRUE(std::filesystem::create_directory(folder)); + EXPECT_THROW(tensor.saveToFile(folder.string()), std::runtime_error); + } +} + +TEST(FileIOTest, overwriteTruncatesOldPayload) { + TemporaryIODirectory directory; + for (const auto *name : {"tensor.bt", "tensor.dtensor"}) { + const auto path = directory.file(name).string(); + DTensor large(std::vector(100, 7), 10, 10); + large.saveToFile(path); + DTensor small(std::vector{1,2}, 2); + small.saveToFile(path); + auto parsed = DTensor::parseFromFile(path); + EXPECT_EQ(parsed.numRows(), 2u); + EXPECT_EQ(parsed.numCols(), 1u); + EXPECT_EQ(parsed.numMats(), 1u); + std::vector actual; + parsed.download(actual); + EXPECT_EQ(actual, (std::vector{1,2})); + } +} diff --git a/test/testTensor.cu b/test/testTensor.cu index bff6ac8..c399701 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -1,6 +1,12 @@ #include #include "../include/tensor.cuh" +#include #include +#include +#include +#include +#include "hostAllocationFailure.h" +#include "sessionAcrossTranslationUnits.h" #define PRECISION_LOW 1e-4 #define PRECISION_HIGH 1e-10 @@ -23,6 +29,344 @@ protected: #define TENSOR_DATA_234APB {8, -4, 12, 6, 6, 17, 41, 7, 5, 5, 19, 17, 14, 13, 5, 11, -8, -4, 6, 8, 8, -2, 8, 13} #define TENSOR_DATA_234AMB {-6, 8, -6, 2, 4, -5, -27, 9, 13, 11, -5, 3, -4, -5, 1, -7, 10, 2, 2, -2, 0, 8, 0, 3}; +/* --------------------------------------- + * Session streams + * + * Configure before main, without initializing CUDA, so filtered/shuffled tests + * have the same stream setup and do not depend on test execution order. + * --------------------------------------- */ + +static const bool testStreamsConfigured = [] { + gputils_test::configureStreamsFromOtherTranslationUnit(2); + return true; +}(); + +TEST(SessionTest, configurationIsSharedAcrossTranslationUnits) { + // The other CUDA translation unit configured the count before main. This + // translation unit must see that count when it first initializes Session. + auto &local = Session::getInstance(); + ASSERT_EQ(2,local.numStreams()); + auto &other = gputils_test::sessionFromOtherTranslationUnit(); + EXPECT_EQ(&local,&other); + EXPECT_EQ(2,other.numStreams()); + EXPECT_EQ(local.cuBlasHandle(1),other.cuBlasHandle(1)); + const size_t before = local.totalAllocatedBytes(); + EXPECT_TRUE(gputils_test::executeOnSecondStreamFromOtherTranslationUnit()); + EXPECT_EQ(before,local.totalAllocatedBytes()); + EXPECT_THROW(gputils_test::configureStreamsFromOtherTranslationUnit(3),std::logic_error); + EXPECT_THROW(gputils_test::configureStreamsFromOtherTranslationUnit(0),std::invalid_argument); + EXPECT_EQ(2,local.numStreams()); +} + +TEST(SessionTest, streamSelectionAndExecution) { + EXPECT_THROW(Session::setStreams(0), std::invalid_argument); + + DTensor A({1, 2, 3, 4}, 2, 2); + DTensor B({5, 6, 7, 8}, 2, 2); + DTensor C(2, 2, 1, true); + A.setStreamIdx(1); + B.setStreamIdx(1); + C.setStreamIdx(1); + + EXPECT_EQ(1, A.streamIdx()); + EXPECT_EQ(1, B.streamIdx()); + EXPECT_EQ(1, C.streamIdx()); + EXPECT_THROW(A.setStreamIdx(2), std::invalid_argument); + EXPECT_THROW(Session::setStreams(3), std::logic_error); + + C.addAB(A, B); + auto selectedRow = C.getRows(1, 1, 0); + EXPECT_EQ(1, selectedRow.streamIdx()); + Session::getInstance().synchronizeStream(1); + std::vector rowResult; + selectedRow.download(rowResult); + EXPECT_EQ((std::vector{34, 46}), rowResult); + std::vector result; + C.download(result); + EXPECT_EQ((std::vector{23, 34, 31, 46}), result); + + DTensor lhs({1, 0, 0, 1, 1, 0, 0, 1}, 2, 2, 2); + DTensor rhs(std::vector{2, 3, 4, 5}, 2, 1, 2); + DTensor status(2); // The solver uses lhs's stream, not the buffer's. + lhs.setStreamIdx(1); + rhs.setStreamIdx(1); + lhs.leastSquaresBatched(rhs, status); + Session::getInstance().synchronizeStream(1); + rhs.download(result); + EXPECT_EQ((std::vector{2, 3, 4, 5}), result); + std::vector codes; + status.download(codes); + EXPECT_EQ((std::vector{0, 0}), codes); +} + +// Temporarily replace handle 1's streams with a nonblocking stream to remove +// implicit legacy-default-stream ordering from the regression tests. +class NonblockingTestStream { +public: + cudaStream_t stream, oldBlasStream, oldSolverStream; + cublasHandle_t blas; + cusolverDnHandle_t solver; + cublasPointerMode_t oldMode; + explicit NonblockingTestStream(size_t index = 1) + : blas(Session::getInstance().cuBlasHandle(index)), solver(Session::getInstance().cuSolverHandle(index)) { + gpuErrChk(cublasGetStream(blas, &oldBlasStream)); + gpuErrChk(cusolverDnGetStream(solver, &oldSolverStream)); + gpuErrChk(cublasGetPointerMode(blas, &oldMode)); + gpuErrChk(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking)); + gpuErrChk(cublasSetStream(blas, stream)); + gpuErrChk(cusolverDnSetStream(solver, stream)); + } + ~NonblockingTestStream() { + gpuErrChk(cudaStreamSynchronize(stream)); + gpuErrChk(cublasSetPointerMode(blas, oldMode)); + gpuErrChk(cublasSetStream(blas, oldBlasStream)); + gpuErrChk(cusolverDnSetStream(solver, oldSolverStream)); + gpuErrChk(cudaStreamDestroy(stream)); + } +}; + +static void CUDART_CB delayTestInput(void *) { + // No CUDA API calls are permitted inside a CUDA host callback. + std::this_thread::sleep_for(std::chrono::milliseconds(20)); +} + +template +void choleskyUsesSelectedStream(T epsilon) { + NonblockingTestStream execution; + for (int mode=0;mode<3;++mode) { + const size_t batch = mode == 0 ? 1 : 2; + std::vector matrix = mode == 2 + ? std::vector{2,1,0,T(std::sqrt(2.0)), 3,1,0,2} + : std::vector{4,2,2,3, 9,3,3,5}; + matrix.resize(4*batch); + std::vector rhs{8,8,0,12}; + rhs.resize(2*batch); + DTensor producerA(matrix,2,2,batch), producerB(rhs,2,1,batch); + DTensor A(2,2,batch,true), B(2,1,batch,true), captured(2,1,batch); + A.setStreamIdx(1); + auto exercise = [&](auto &factoriser) { + EXPECT_EQ(1,factoriser.info().streamIdx()); + gpuErrChk(cudaLaunchHostFunc(execution.stream,delayTestInput,nullptr)); + gpuErrChk(cudaMemcpyAsync(A.raw(),producerA.raw(),matrix.size()*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + gpuErrChk(cudaMemcpyAsync(B.raw(),producerB.raw(),rhs.size()*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + if (mode != 2) factoriser.factorise(); + factoriser.solve(B); + gpuErrChk(cudaMemcpyAsync(captured.raw(),B.raw(),rhs.size()*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + std::vector codes; + factoriser.info().download(codes); // Must wait for the selected solver stream. + EXPECT_EQ(0,codes[0]); // potrsBatched reports one call-level status. + std::vector actual; + captured.download(actual); + EXPECT_NEAR(T(1),actual[0],epsilon); + EXPECT_NEAR(T(2),actual[1],epsilon); + if (batch == 2) { + EXPECT_NEAR(T(-1),actual[2],epsilon); + EXPECT_NEAR(T(3),actual[3],epsilon); + } + EXPECT_EQ(0,B.streamIdx()); // The matrix, not the RHS, selects execution. + }; + if (mode == 0) { + CholeskyFactoriser factoriser(A); + exercise(factoriser); + } else { + CholeskyBatchFactoriser factoriser(A,mode == 2); + exercise(factoriser); + } + } +} + +TEST_F(TensorTest, choleskyOrdersWorkOnSelectedStream) { + choleskyUsesSelectedStream(PRECISION_LOW); + choleskyUsesSelectedStream(PRECISION_HIGH); +} + +template +void transfersUseSelectedStream() { + NonblockingTestStream execution; + const std::vector data{1,2,3,4}; + DTensor producer(data,2,2), tensor(2,2,1,true); + tensor.setStreamIdx(1); + auto prepare = [&] { + gpuErrChk(cudaLaunchHostFunc(execution.stream,delayTestInput,nullptr)); + gpuErrChk(cudaMemcpyAsync(tensor.raw(),producer.raw(),4*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + }; + prepare(); + std::vector actual; + tensor.download(actual); // No explicit synchronization by the caller. + EXPECT_EQ(data,actual); + prepare(); + EXPECT_EQ(T(4),tensor(1,1)); + for (auto mode : {columnMajor,rowMajor}) { + gpuErrChk(cudaLaunchHostFunc(execution.stream,delayTestInput,nullptr)); + gpuErrChk(cudaMemsetAsync(tensor.raw(),0,4*sizeof(T),execution.stream)); + tensor.upload(data,mode); + gpuErrChk(cudaStreamSynchronize(execution.stream)); + tensor.download(actual); + EXPECT_EQ((mode == columnMajor ? data : std::vector{1,3,2,4}),actual); + } + prepare(); + auto transposed = tensor.tr(); + EXPECT_EQ(1,transposed.streamIdx()); + transposed.download(actual); + EXPECT_EQ((std::vector{1,3,2,4}),actual); + EXPECT_EQ(T(4),tensor.maxAbs()); + EXPECT_EQ(T(1),tensor.minAbs()); +} + +TEST_F(TensorTest, hostTransfersAndTransposeUseSelectedStream) { + transfersUseSelectedStream(); + transfersUseSelectedStream(); +} + +template +void deviceCopiesUseSelectedStreams() { + NonblockingTestStream execution; + DTensor source(4,1,1,true), destination(6), producer(std::vector{1,2,3,4},4); + source.setStreamIdx(1); + cudaStream_t destinationStream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(),&destinationStream)); + gpuErrChk(cudaLaunchHostFunc(destinationStream,delayTestInput,nullptr)); + gpuErrChk(cudaMemsetAsync(destination.raw(),0,6*sizeof(T),destinationStream)); + gpuErrChk(cudaLaunchHostFunc(execution.stream,delayTestInput,nullptr)); + gpuErrChk(cudaMemcpyAsync(source.raw(),producer.raw(),4*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + source.deviceCopyTo(destination); + // A consumer on the destination's stream can run immediately after return. + destination *= T(2); + std::vector actual; + destination.download(actual); + EXPECT_EQ((std::vector{2,4,6,8,0,0}),actual); + EXPECT_EQ(1,source.streamIdx()); + EXPECT_EQ(0,destination.streamIdx()); + destination.setStreamIdx(1); + source.deviceCopyTo(destination); + destination.download(actual); + EXPECT_EQ((std::vector{1,2,3,4,0,0}),actual); +} + +TEST_F(TensorTest, deviceCopiesOrderBothStreams) { + deviceCopiesUseSelectedStreams(); + deviceCopiesUseSelectedStreams(); +} + +template +void qrUsesSelectedStream(T epsilon) { + NonblockingTestStream execution; + DTensor source(std::vector{1,0,1,0,1,1}, 3, 2); + DTensor sourceB(std::vector{2,3,5}, 3); + DTensor A(3,2,1,true), B(3,1,1,true), captured(3), Q(3,2), R(2,2); + A.setStreamIdx(1); + QRFactoriser qr(A); + gpuErrChk(cudaLaunchHostFunc(execution.stream, delayTestInput, nullptr)); + gpuErrChk(cudaMemcpyAsync(A.raw(), source.raw(), 6*sizeof(T), cudaMemcpyDeviceToDevice, execution.stream)); + gpuErrChk(cudaMemcpyAsync(B.raw(), sourceB.raw(), 3*sizeof(T), cudaMemcpyDeviceToDevice, execution.stream)); + qr.factorise(); + qr.leastSquares(B); + // A downstream consumer must see the completed solution on this stream. + gpuErrChk(cudaMemcpyAsync(captured.raw(), B.raw(), 3*sizeof(T), cudaMemcpyDeviceToDevice, execution.stream)); + gpuErrChk(cudaStreamSynchronize(execution.stream)); + std::vector solution; + captured.download(solution); + EXPECT_NEAR(2, solution[0], epsilon); + EXPECT_NEAR(3, solution[1], epsilon); + EXPECT_EQ(0, qr.info()(0)); + + qr.getQR(Q,R); + std::vector q, r, expected; + Q.download(q); R.download(r); source.download(expected); + for (size_t c=0;c<2;++c) for (size_t row=0;row<3;++row) { + T actual=0; + for (size_t k=0;k<2;++k) actual += q[row+3*k]*r[k+2*c]; + EXPECT_NEAR(expected[row+3*c], actual, epsilon); + } +} + +TEST_F(TensorTest, qrOrdersWorkOnSelectedStream) { + qrUsesSelectedStream(PRECISION_LOW); + qrUsesSelectedStream(PRECISION_HIGH); +} + +template +void givensUsesSelectedStream(T epsilon) { + NonblockingTestStream execution; + DTensor source(std::vector{3,4,0,5},2,2); + DTensor A(std::vector{1,0,0,1},2,2), captured(2,2); + A.setStreamIdx(1); + GivensAnnihilator rotation(A); + for (auto mode : {CUBLAS_POINTER_MODE_HOST, CUBLAS_POINTER_MODE_DEVICE}) { + A.upload(std::vector{1,0,0,1}); + gpuErrChk(cublasSetPointerMode(execution.blas,mode)); + gpuErrChk(cudaLaunchHostFunc(execution.stream,delayTestInput,nullptr)); + gpuErrChk(cudaMemcpyAsync(A.raw(),source.raw(),4*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + rotation.annihilate(0,1,0); + gpuErrChk(cudaMemcpyAsync(captured.raw(),A.raw(),4*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + gpuErrChk(cudaStreamSynchronize(execution.stream)); + cublasPointerMode_t actualMode; + gpuErrChk(cublasGetPointerMode(execution.blas,&actualMode)); + EXPECT_EQ(mode,actualMode); + std::vector actual; + captured.download(actual); + const std::vector expected{5,0,4,3}; + for(size_t i=0;i(PRECISION_LOW); + givensUsesSelectedStream(PRECISION_HIGH); +} + +template +void streamSelectionPreservesStorage() { + static_assert(std::is_same_v &>().setStreamIdx(0)), DTensor &>); + static_assert(std::is_same_v &&>().setStreamIdx(0)), DTensor>); + auto &session = Session::getInstance(); + const size_t before = session.totalAllocatedBytes(); + const std::vector data{1, 2, 3, 4, 5, 6, 7, 8}; + { + DTensor source(data, 2, 2, 2); + auto *raw = source.raw(); + auto *matrices = source.ptrMatrices(); + const size_t allocated = session.totalAllocatedBytes(); + auto &reference = source.setStreamIdx(0); + EXPECT_EQ(&source, &reference); + EXPECT_EQ(raw, source.raw()); + EXPECT_EQ(matrices, source.ptrMatrices()); + EXPECT_EQ(allocated, session.totalAllocatedBytes()); + + EXPECT_THROW(std::move(source).setStreamIdx(session.numStreams()), std::invalid_argument); + EXPECT_EQ(raw, source.raw()); + EXPECT_EQ(0, source.streamIdx()); + + auto moved = std::move(source).setStreamIdx(0); + EXPECT_EQ(raw, moved.raw()); + EXPECT_EQ(matrices, moved.ptrMatrices()); + EXPECT_EQ(nullptr, source.raw()); + EXPECT_EQ(0, source.numEl()); + EXPECT_EQ(allocated, session.totalAllocatedBytes()); + std::vector downloaded; + moved.download(downloaded); + EXPECT_EQ(data, downloaded); + } + EXPECT_EQ(before, session.totalAllocatedBytes()); + { + auto chained = DTensor(data, 2, 2, 2).setStreamIdx(0).setStreamIdx(0); + EXPECT_EQ(2, chained.numRows()); + EXPECT_EQ(2, chained.numCols()); + EXPECT_EQ(2, chained.numMats()); + EXPECT_EQ(0, chained.streamIdx()); + std::vector downloaded; + chained.download(downloaded); + EXPECT_EQ(data, downloaded); + } + EXPECT_EQ(before, session.totalAllocatedBytes()); +} + +TEST_F(TensorTest, streamSelectionPreservesStorage) { + streamSelectionPreservesStorage(); + streamSelectionPreservesStorage(); + streamSelectionPreservesStorage(); +} + /* --------------------------------------- * Zero Tensor (Constructor) * --------------------------------------- */ @@ -45,6 +389,71 @@ TEST_F(TensorTest, tensorConstructionZero) { tensorConstructionZero(); } +namespace { +thread_local size_t failureDeviceBytes = 0; +thread_local size_t failureHostBytes = 0; +thread_local size_t observedDeviceBytes = 0; +thread_local size_t observedHostBytes = 0; + +bool failAfterDeviceAllocation(size_t bytes) noexcept { + const size_t allocated = Session::getInstance().totalAllocatedBytes(); + if (bytes != failureHostBytes || allocated != failureDeviceBytes) return false; + observedDeviceBytes = allocated; + observedHostBytes = bytes; + return true; +} +} + +template +void allocatingConstructorsReleaseStorageOnHostFailure() { + auto &session = Session::getInstance(); // Initialize CUDA before arming failure. + const size_t initialBytes = session.totalAllocatedBytes(); + const std::vector data(30,T(7)); + // Warm up CUDA/library initialization and all transfer paths before testing. + { DTensor warmup(data,2,3,5,rowMajor); } + for (int scenario=0;scenario<4;++scenario) { + SCOPED_TRACE("allocation-failure scenario " + std::to_string(scenario)); + const bool sized = scenario < 2; + const size_t mats = scenario == 3 ? 1 : 5; + const size_t elements = 6*mats; + const std::vector input(elements,T(7)); + const size_t baseline = session.totalAllocatedBytes(); + failureDeviceBytes = baseline + elements*sizeof(T) + (mats > 1 ? mats*sizeof(T *) : 0); + // Sized constructors fail in the pointer-table host vector. Vector + // constructors fail in the row-major conversion vector. + failureHostBytes = sized ? mats*sizeof(T *) : elements*sizeof(T); + observedDeviceBytes = observedHostBytes = 0; + bool caught = false; + { + gputils_test::FailHostAllocation failure(failAfterDeviceAllocation); + try { + if (sized) { DTensor tensor(2,3,mats,scenario == 1); } + else { DTensor tensor(input,2,3,mats,rowMajor); } + } catch (const std::bad_alloc &) { caught = true; } + } + EXPECT_TRUE(caught); + EXPECT_EQ(failureHostBytes,observedHostBytes); + EXPECT_EQ(failureDeviceBytes,observedDeviceBytes); // GPU allocations existed at failure. + EXPECT_GT(observedDeviceBytes,baseline); + EXPECT_EQ(baseline,session.totalAllocatedBytes()); + // The original exception is recoverable and CUDA remains usable. + { + DTensor recovered(input,2,3,mats,rowMajor); + std::vector actual; + recovered.download(actual); + EXPECT_EQ(input,actual); // Constant input has identical row/column layout. + } + EXPECT_EQ(baseline,session.totalAllocatedBytes()); + } + EXPECT_EQ(initialBytes,session.totalAllocatedBytes()); +} + +TEST_F(TensorTest, constructorsReleaseStorageOnHostAllocationFailure) { + allocatingConstructorsReleaseStorageOnHostFailure(); + allocatingConstructorsReleaseStorageOnHostFailure(); + allocatingConstructorsReleaseStorageOnHostFailure(); +} + /* --------------------------------------- * Row- and column-major data * --------------------------------------- */ @@ -191,14 +600,34 @@ TEST_F(TensorTest, parseTensorFromFileBinary) { * --------------------------------------- */ TEST_F(TensorTest, parseTensorFromBinaryPython) { - std::string fName = "../../python/b_d.bt"; - DTensor b = DTensor::parseFromFile(fName); + // Write the binary layout emitted by write_array_to_gputils_binary_file() + // directly, so this test does not depend on the Python test suite having + // run first or on its working directory. + const auto fName = std::filesystem::path(testing::TempDir()) / "gputils_python_fixture.bt"; + const uint64_t numRows = 3, numCols = 3, numMats = 2; + const std::vector pythonBinaryData = { + 1, 7, 13, 3, 9, 15, 5, 11, 17, + 2, 8, 14, 4, 10, 16, 6, 12, 18 + }; + { + std::ofstream file(fName, std::ios::binary); + ASSERT_TRUE(file.is_open()); + file.write(reinterpret_cast(&numRows), sizeof(numRows)); + file.write(reinterpret_cast(&numCols), sizeof(numCols)); + file.write(reinterpret_cast(&numMats), sizeof(numMats)); + file.write(reinterpret_cast(pythonBinaryData.data()), + static_cast(pythonBinaryData.size() * sizeof(double))); + ASSERT_TRUE(file.good()); + } + + DTensor b = DTensor::parseFromFile(fName.string()); for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < 3; j++) { EXPECT_NEAR(1 + 2*j + 6*i, b(i, j, 0), PRECISION_HIGH); EXPECT_NEAR(2 + 2*j + 6*i, b(i, j, 1), PRECISION_HIGH); } } + std::filesystem::remove(fName); } @@ -337,7 +766,8 @@ TEMPLATE_WITH_TYPE_T void tensorSlicingConstructorAxis0() { std::vector data = TENSOR_DATA_234A; DTensor tenz(data, 2, 3, 4); - DTensor tenzSlice(tenz, 0, 2, 3); // elements 2..3 + DTensor column(tenz, 1, 1, 1); + DTensor tenzSlice(column, 0, 0, 1); // both rows of the second column EXPECT_EQ(2, tenzSlice.numRows()); EXPECT_EQ(1, tenzSlice.numCols()); EXPECT_EQ(1, tenzSlice.numMats()); @@ -675,10 +1105,17 @@ void tensorAssignmentOperator() { DTensor tenz(data, 2, 3, 4); DTensor other; other = tenz; - EXPECT_EQ(tenz.raw(), other.raw()); + EXPECT_NE(tenz.raw(), other.raw()); + EXPECT_NE(tenz.ptrMatrices(), other.ptrMatrices()); EXPECT_EQ(2, other.numRows()); EXPECT_EQ(3, other.numCols()); EXPECT_EQ(4, other.numMats()); + std::vector actual; + other.download(actual); + EXPECT_EQ(data, actual); + tenz.upload(std::vector(data.size(), T(0))); + other.download(actual); + EXPECT_EQ(data, actual); } TEST_F(TensorTest, tensorAssignmentOperator) { @@ -687,6 +1124,124 @@ TEST_F(TensorTest, tensorAssignmentOperator) { tensorAssignmentOperator(); } +template +void tensorAssignmentOwnership() { + static_assert(std::is_nothrow_move_constructible_v>); + static_assert(std::is_nothrow_move_assignable_v>); + auto &session = Session::getInstance(); + const size_t baseline = session.totalAllocatedBytes(); + const std::vector data{1,2,3,4,5,6,7,8}; + const size_t bytes = data.size()*sizeof(T) + 2*sizeof(T *); + { + DTensor source(data,2,2,2), target(13); + target = source; + EXPECT_EQ(baseline+2*bytes,session.totalAllocatedBytes()); + const DTensor &self = target; + auto *raw = target.raw(); + target = self; + EXPECT_EQ(raw,target.raw()); + EXPECT_EQ(baseline+2*bytes,session.totalAllocatedBytes()); + std::vector matrices(2); + gpuErrChk(cudaMemcpy(matrices.data(),target.ptrMatrices(),2*sizeof(T *),cudaMemcpyDeviceToHost)); + EXPECT_EQ(target.raw(),matrices[0]); + EXPECT_EQ(target.raw()+4,matrices[1]); + + auto *sourceRaw = source.raw(); + auto *sourceMatrices = source.ptrMatrices(); + target = std::move(source); + EXPECT_EQ(sourceRaw,target.raw()); + EXPECT_EQ(sourceMatrices,target.ptrMatrices()); + EXPECT_EQ(nullptr,source.raw()); + EXPECT_EQ(nullptr,source.ptrMatrices()); + EXPECT_EQ(0,source.numEl()); + EXPECT_EQ(baseline+bytes,session.totalAllocatedBytes()); + auto &moveSelf = target; + target = std::move(moveSelf); + EXPECT_EQ(sourceRaw,target.raw()); + std::vector actual; + target.download(actual); + EXPECT_EQ(data,actual); + source = target; // A moved-from tensor can be assigned again. + EXPECT_NE(source.raw(),target.raw()); + } + EXPECT_EQ(baseline,session.totalAllocatedBytes()); + { + DTensor target(13); + { DTensor source(data,2,2,2); target = source; } + EXPECT_EQ(baseline+bytes,session.totalAllocatedBytes()); + std::vector actual; + target.download(actual); + EXPECT_EQ(data,actual); // Source destruction cannot invalidate the copy. + DTensor empty; + DTensor emptyCopy(empty); + EXPECT_EQ(nullptr,emptyCopy.raw()); + target = empty; + EXPECT_EQ(nullptr,target.raw()); + EXPECT_EQ(0,target.numEl()); + EXPECT_EQ(baseline,session.totalAllocatedBytes()); + target = DTensor(data,2,2,2); + EXPECT_EQ(baseline+bytes,session.totalAllocatedBytes()); + target = std::move(emptyCopy); + EXPECT_EQ(nullptr,target.raw()); + EXPECT_EQ(baseline,session.totalAllocatedBytes()); + } + { + DTensor owner(data,2,2,2); + DTensor view(owner,2,0,0); + view = owner; // Detach the view without freeing the owner's storage. + EXPECT_NE(owner.raw(),view.raw()); + EXPECT_EQ(baseline+2*bytes,session.totalAllocatedBytes()); + DTensor source(data,2,2,2); + DTensor moveDestination(owner,2,0,0); + auto *raw = source.raw(); + moveDestination = std::move(source); + EXPECT_EQ(raw,moveDestination.raw()); + std::vector actual; + owner.download(actual); + EXPECT_EQ(data,actual); + } + EXPECT_EQ(baseline,session.totalAllocatedBytes()); + { + DTensor owner(data,2,2,2); + DTensor view(owner,2,1,1); + owner = view; // Copy before releasing the storage being viewed. + std::vector actual; + owner.download(actual); + EXPECT_EQ((std::vector{5,6,7,8}),actual); + EXPECT_EQ(1,owner.numMats()); + } + EXPECT_EQ(baseline,session.totalAllocatedBytes()); +} + +TEST_F(TensorTest, assignmentReleasesAndTransfersOwnership) { + tensorAssignmentOwnership(); + tensorAssignmentOwnership(); + tensorAssignmentOwnership(); +} + +TEST_F(TensorTest, copiesRespectSelectedStream) { + NonblockingTestStream execution; + const std::vector data{1,2,3,4,5,6,7,8}; + DTensor producer(data,2,2,2), source(2,2,2,true), assigned(3); + source.setStreamIdx(1); + gpuErrChk(cudaLaunchHostFunc(execution.stream,delayTestInput,nullptr)); + gpuErrChk(cudaMemcpyAsync(source.raw(),producer.raw(),data.size()*sizeof(double),cudaMemcpyDeviceToDevice,execution.stream)); + DTensor constructed(source); + assigned = source; + EXPECT_EQ(1,constructed.streamIdx()); + EXPECT_EQ(1,assigned.streamIdx()); + std::vector actual; + constructed.download(actual); + EXPECT_EQ(data,actual); + assigned.download(actual); + EXPECT_EQ(data,actual); + DTensor moved; + auto *pointer = assigned.raw(); + moved = std::move(assigned); + EXPECT_EQ(1,moved.streamIdx()); + EXPECT_EQ(pointer,moved.raw()); +} + /* --------------------------------------- * Tensor times-equals scalar * --------------------------------------- */ @@ -918,6 +1473,11 @@ void tensorGetRows() { std::vector actual1(6); Ar1.download(actual1); EXPECT_EQ(expected1, actual1); + + auto fullMatrix = A.getRows(0, 2, 1); + std::vector actualFull; + fullMatrix.download(actualFull); + EXPECT_EQ((std::vector{5, 6, 7, 8, 9, 10, 11, 12, 13}), actualFull); } TEST_F(TensorTest, tensorGetRows) { @@ -925,6 +1485,69 @@ TEST_F(TensorTest, tensorGetRows) { tensorGetRows(); } +template +void tensorGetRowsInto() { + std::vector data(24); + DTensor source(4,3,2); + for (size_t count : {size_t(1),size_t(2),size_t(4)}) { + DTensor destination(count,3); + auto *pointer = destination.raw(); + const size_t allocated = Session::getInstance().totalAllocatedBytes(); + const size_t first = count == 4 ? 0 : 1; + for (size_t repeat=0;repeat<3;++repeat) { + std::iota(data.begin(),data.end(),T(10*repeat)); + source.upload(data); + source.getRowsInto(destination,first,first+count-1,repeat%2); + Session::getInstance().synchronizeStream(source.streamIdx()); + EXPECT_EQ(pointer,destination.raw()); + EXPECT_EQ(allocated,Session::getInstance().totalAllocatedBytes()); + std::vector actual; + destination.download(actual); + for (size_t col=0;col<3;++col) for(size_t row=0;row destination(2,3), wrongRows(1,3), wrongColumns(2,2), wrongBatch(2,3,2), empty; + EXPECT_THROW(source.getRowsInto(destination,2,1,0),std::invalid_argument); + EXPECT_THROW(source.getRowsInto(destination,0,4,0),std::invalid_argument); + EXPECT_THROW(source.getRowsInto(destination,0,1,2),std::invalid_argument); + EXPECT_THROW(source.getRowsInto(wrongRows,0,1,0),std::invalid_argument); + EXPECT_THROW(source.getRowsInto(wrongColumns,0,1,0),std::invalid_argument); + EXPECT_THROW(source.getRowsInto(wrongBatch,0,1,0),std::invalid_argument); + EXPECT_THROW(source.getRowsInto(empty,0,1,0),std::invalid_argument); + DTensor overlapping(source,2,0,0); + EXPECT_THROW(source.getRowsInto(overlapping,0,3,0),std::invalid_argument); + EXPECT_THROW(overlapping.getRowsInto(overlapping,0,3,0),std::invalid_argument); +} + +TEST_F(TensorTest, getRowsIntoReusesDestination) { + tensorGetRowsInto(); + tensorGetRowsInto(); + tensorGetRowsInto(); +} + +template +void getRowsIntoUsesSelectedStream() { + NonblockingTestStream execution; + DTensor producer(std::vector{1,2,3,4,5,6},3,2); + DTensor source(3,2,1,true), destination(2,2), captured(2,2); + source.setStreamIdx(1); + gpuErrChk(cudaLaunchHostFunc(execution.stream,delayTestInput,nullptr)); + gpuErrChk(cudaMemcpyAsync(source.raw(),producer.raw(),6*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + source.getRowsInto(destination,1,2,0); + gpuErrChk(cudaMemcpyAsync(captured.raw(),destination.raw(),4*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + gpuErrChk(cudaStreamSynchronize(execution.stream)); + EXPECT_EQ(0,destination.streamIdx()); // Caller configuration is preserved. + std::vector actual; + captured.download(actual); + EXPECT_EQ((std::vector{2,3,5,6}),actual); +} + +TEST_F(TensorTest, getRowsIntoOrdersWorkOnSelectedStream) { + getRowsIntoUsesSelectedStream(); + getRowsIntoUsesSelectedStream(); +} + /* --------------------------------------- * Tensor: transpose @@ -991,6 +1614,233 @@ TEST_F(TensorTest, tensorBytesAllocatedDeallocated) { tensorBytesAllocatedDeallocated(); } +TEST(AllocationCounterTest, handlesLargeCountsAndRejectsOverflowAndUnderflow) { + gputils_detail::AllocationCounter counter; + const size_t large = static_cast(INT_MAX) + 4096; + counter.recordAllocation(large); + counter.recordAllocation(17); + EXPECT_EQ(large+17,counter.bytes()); + counter.recordDeallocation(large); + EXPECT_EQ(17,counter.bytes()); + EXPECT_THROW(counter.recordDeallocation(18),std::logic_error); + EXPECT_EQ(17,counter.bytes()); + counter.recordDeallocation(17); + counter.recordAllocation(std::numeric_limits::max()); + EXPECT_THROW(counter.recordAllocation(1),std::overflow_error); + EXPECT_EQ(std::numeric_limits::max(),counter.bytes()); + counter.recordDeallocation(std::numeric_limits::max()); + EXPECT_EQ(0,counter.bytes()); +} + +template +void reshapeAccountsForOwnedCapacity() { + auto &session = Session::getInstance(); + const size_t baseline = session.totalAllocatedBytes(); + const size_t dataBytes = 24*sizeof(T); + std::vector data(24); + std::iota(data.begin(),data.end(),T(1)); + { + DTensor tensor(data,2,3,4); + auto *raw = tensor.raw(); + auto *pointers = tensor.ptrMatrices(); + EXPECT_EQ(baseline+dataBytes+4*sizeof(T *),session.totalAllocatedBytes()); + tensor.reshape(4,3,2); + tensor.reshape(8,3,1); + tensor.reshape(2,4,3); + EXPECT_EQ(raw,tensor.raw()); + EXPECT_EQ(pointers,tensor.ptrMatrices()); + EXPECT_EQ(baseline+dataBytes+4*sizeof(T *),session.totalAllocatedBytes()); + std::vector matrices(3); + gpuErrChk(cudaMemcpy(matrices.data(),tensor.ptrMatrices(),3*sizeof(T *),cudaMemcpyDeviceToHost)); + EXPECT_EQ(raw,matrices[0]); + EXPECT_EQ(raw+8,matrices[1]); + EXPECT_EQ(raw+16,matrices[2]); + + tensor.reshape(2,2,6); + pointers = tensor.ptrMatrices(); + EXPECT_EQ(baseline+dataBytes+6*sizeof(T *),session.totalAllocatedBytes()); + tensor.reshape(4,3,2); + tensor.reshape(2,3,4); + EXPECT_EQ(pointers,tensor.ptrMatrices()); + EXPECT_EQ(baseline+dataBytes+6*sizeof(T *),session.totalAllocatedBytes()); + EXPECT_THROW(tensor.reshape(2,3,5),std::invalid_argument); + std::vector actual; + tensor.download(actual); + EXPECT_EQ(data,actual); + + DTensor moved(std::move(tensor)); + DTensor assigned(7); + assigned = std::move(moved); + EXPECT_EQ(baseline+dataBytes+6*sizeof(T *),session.totalAllocatedBytes()); + DTensor other(5); + assigned.swap(other); + other.reshape(24,1,1); + EXPECT_EQ(baseline+dataBytes+6*sizeof(T *)+5*sizeof(T),session.totalAllocatedBytes()); + DTensor copy(other); // Copy owns only the storage its current shape needs. + EXPECT_EQ(baseline+2*dataBytes+6*sizeof(T *)+5*sizeof(T),session.totalAllocatedBytes()); + } + EXPECT_EQ(baseline,session.totalAllocatedBytes()); + { + DTensor owner(data,2,3,4); + auto *ownerPointers = owner.ptrMatrices(); + { + DTensor view(owner,2,1,2); + view.reshape(3,2,2); // Own a new pointer array, but keep borrowing data. + EXPECT_NE(ownerPointers+1,view.ptrMatrices()); + EXPECT_EQ(baseline+dataBytes+6*sizeof(T *),session.totalAllocatedBytes()); + std::vector matrices(2); + gpuErrChk(cudaMemcpy(matrices.data(),view.ptrMatrices(),2*sizeof(T *),cudaMemcpyDeviceToHost)); + EXPECT_EQ(owner.raw()+6,matrices[0]); + EXPECT_EQ(owner.raw()+12,matrices[1]); + } + EXPECT_EQ(ownerPointers,owner.ptrMatrices()); + EXPECT_EQ(baseline+dataBytes+4*sizeof(T *),session.totalAllocatedBytes()); + } + EXPECT_EQ(baseline,session.totalAllocatedBytes()); +} + +TEST_F(TensorTest, reshapeTracksAllocationCapacity) { + reshapeAccountsForOwnedCapacity(); + reshapeAccountsForOwnedCapacity(); + reshapeAccountsForOwnedCapacity(); +} + +/* --------------------------------------- + * Input validation + * --------------------------------------- */ + +TEST_F(TensorTest, sliceConstructorRejectsInvalidAxisAndRanges) { + DTensor tensor(2, 3, 2); + + EXPECT_THROW(DTensor(tensor, 3, 0, 0), std::invalid_argument); + EXPECT_THROW(DTensor(tensor, 0, 1, 0), std::invalid_argument); + EXPECT_THROW(DTensor(tensor, 0, 0, 2), std::invalid_argument); + EXPECT_THROW(DTensor(tensor, 1, 0, 3), std::invalid_argument); + EXPECT_THROW(DTensor(tensor, 2, 0, 2), std::invalid_argument); +} + +TEST_F(TensorTest, givensRotationsRejectBatches) { + const double c = 1.0; + const double minusS = 0.0; + DTensor batchedMatrix(2, 3, 2); + + EXPECT_THROW(batchedMatrix.applyRightGivensRotation(0, 1, &c, &minusS), std::invalid_argument); + EXPECT_THROW(batchedMatrix.applyLeftGivensRotation(0, 1, &c, &minusS), std::invalid_argument); +} + +TEST_F(TensorTest, tensorOperationsRejectIncompatibleDimensions) { + DTensor twoByTwo(2, 2); + DTensor twoByThree(2, 3); + DTensor tooSmall(3); + + EXPECT_THROW(twoByTwo.dotF(twoByThree), std::invalid_argument); + EXPECT_THROW(twoByTwo.deviceCopyTo(tooSmall), std::invalid_argument); +} + +TEST_F(TensorTest, tensorSetupRejectsInvalidDimensionsAndRanges) { + EXPECT_THROW(DTensor(0, 1, 1), std::invalid_argument); + EXPECT_THROW(DTensor(static_cast(INT_MAX) + 1, 1, 1), std::invalid_argument); + EXPECT_THROW(DTensor(static_cast(INT_MAX), static_cast(INT_MAX), 2), std::invalid_argument); + EXPECT_THROW(DTensor((std::vector{1}), 2), std::invalid_argument); + EXPECT_THROW(DTensor::createRandomTensor(2, 1, 1, 2.0, 1.0), std::invalid_argument); + DTensor tensor(2, 3, 2); + EXPECT_THROW(tensor.reshape(static_cast(INT_MAX), static_cast(INT_MAX), 2), std::invalid_argument); + EXPECT_THROW(tensor.getRows(1, 0, 0), std::invalid_argument); + EXPECT_THROW(tensor.getRows(0, 2, 0), std::invalid_argument); + EXPECT_THROW(tensor.getRows(0, 1, 2), std::invalid_argument); + EXPECT_THROW(tensor(2, 0, 0), std::invalid_argument); + EXPECT_THROW(tensor(0, 3, 0), std::invalid_argument); + EXPECT_THROW(tensor(0, 0, 2), std::invalid_argument); + EXPECT_THROW(tensor.upload(std::vector(12), static_cast(99)), std::invalid_argument); +} + +TEST_F(TensorTest, factoriserSetupRejectsEmptyAndIncompatibleMatrices) { + DTensor empty; + DTensor wide(2, 3); + DTensor nonsquare(3, 2, 2); + EXPECT_THROW(Svd svd(empty), std::invalid_argument); + EXPECT_THROW(CholeskyFactoriser cholesky(empty), std::invalid_argument); + EXPECT_THROW(QRFactoriser qr(wide), std::invalid_argument); + EXPECT_THROW(CholeskyBatchFactoriser cholesky(nonsquare), std::invalid_argument); + EXPECT_THROW(GivensAnnihilator givens(empty), std::invalid_argument); +} + +TEST_F(TensorTest, fileParsingRejectsMalformedDataAndHeaders) { + const auto path = std::filesystem::path(testing::TempDir()) / "gputils_invalid_input.dtensor"; + for (const auto &content : {"2\n1\n1\n1\n", "2\n1\n1\n1\ninvalid\n", "1\n1\n1\n1\n2\n", "2\n"}) { + { std::ofstream file(path); ASSERT_TRUE(file.is_open()); file << content; } + EXPECT_THROW(DTensor::parseFromFile(path.string()), std::invalid_argument); + } + std::filesystem::remove(path); + const auto binaryPath = std::filesystem::path(testing::TempDir()) / "gputils_invalid_header.bt"; + { std::ofstream file(binaryPath, std::ios::binary); ASSERT_TRUE(file.is_open()); file.put('\1'); } + EXPECT_THROW(DTensor::parseFromFile(binaryPath.string()), std::invalid_argument); + std::filesystem::remove(binaryPath); +} + +TEST_F(TensorTest, leastSquaresRejectsIncompatibleDimensions) { + DTensor tallMatrix(3, 2, 2); + DTensor wrongRows(2, 1, 2); + DTensor wrongColumns(3, 2, 2); + DTensor wrongBatchSize(3, 1, 1); + DTensor wideMatrix(2, 3, 2); + DTensor compatibleRhs(2, 1, 2); + + EXPECT_THROW(tallMatrix.leastSquaresBatched(wrongRows), std::invalid_argument); + EXPECT_THROW(tallMatrix.leastSquaresBatched(wrongColumns), std::invalid_argument); + EXPECT_THROW(tallMatrix.leastSquaresBatched(wrongBatchSize), std::invalid_argument); + EXPECT_THROW(wideMatrix.leastSquaresBatched(compatibleRhs), std::invalid_argument); +} + +/* --------------------------------------- + * Matrix multiplication coefficients + * --------------------------------------- */ + +TEST_F(TensorTest, tensorAddABUsesAlphaBetaForASingleMatrix) { + DTensor A({1, 2, 3, 4}, 2, 2); + DTensor B({5, 6, 7, 8}, 2, 2); + DTensor C({1, 2, 3, 4}, 2, 2); + + C.addAB(A, B, 2.0, -1.0); + + std::vector result; + C.download(result); + EXPECT_EQ((std::vector{45, 66, 59, 88}), result); +} + +/* --------------------------------------- + * Malformed input files + * --------------------------------------- */ + +TEST_F(TensorTest, parseTensorRejectsTruncatedBinaryFile) { + const auto path = std::filesystem::path(testing::TempDir()) / "gputils_truncated_fixture.bt"; + const uint64_t rows = 2, cols = 2, mats = 1; + const double onlyOneElement = 1.0; + { + std::ofstream file(path, std::ios::binary); + ASSERT_TRUE(file.is_open()); + file.write(reinterpret_cast(&rows), sizeof(rows)); + file.write(reinterpret_cast(&cols), sizeof(cols)); + file.write(reinterpret_cast(&mats), sizeof(mats)); + file.write(reinterpret_cast(&onlyOneElement), sizeof(onlyOneElement)); + } + + EXPECT_THROW(DTensor::parseFromFile(path.string()), std::invalid_argument); + std::filesystem::remove(path); +} + +TEST_F(TensorTest, parseTensorRejectsMalformedTextFile) { + const auto path = std::filesystem::path(testing::TempDir()) / "gputils_malformed_fixture.dtensor"; + { + std::ofstream file(path); + ASSERT_TRUE(file.is_open()); + file << "not-a-dimension\n2\n1\n1\n2\n3\n4\n"; + } + + EXPECT_THROW(DTensor::parseFromFile(path.string()), std::invalid_argument); + std::filesystem::remove(path); +} + /* ================================================================================================ @@ -1038,6 +1888,134 @@ TEST_F(LeastSquaresTest, tensorLS1) { tensorLeastSquares1(PRECISION_HIGH); } +// A batch of one is valid even though ordinary single-matrix tensors do not +// allocate a device array of matrix pointers. Exercise the public solver API. +template +void leastSquaresSingleMatrixBatch(T epsilon) { + NonblockingTestStream execution; + auto &session = Session::getInstance(); + for (size_t batch : {size_t(1), size_t(2)}) { + for (size_t rows : {size_t(2), size_t(3)}) { + std::vector data(rows * 2 * batch, T(0)), rhs(rows * batch, T(0)); + for (size_t k = 0; k < batch; ++k) { + data[k * rows * 2] = T(2); + data[k * rows * 2 + rows + 1] = T(4); + rhs[k * rows] = T(2); rhs[k * rows + 1] = T(8); + if (rows == 3) { + data[k * rows * 2 + 2] = data[k * rows * 2 + rows + 2] = T(1); + rhs[k * rows + 2] = T(3); + } + } + DTensor A(data, rows, 2, batch), otherA(data, rows, 2, batch); + DTensor b(rhs, rows, 1, batch), otherB(rhs, rows, 1, batch); + A.setStreamIdx(1); otherA.setStreamIdx(1); b.setStreamIdx(1); otherB.setStreamIdx(1); + const size_t before = session.totalAllocatedBytes(); + { + LeastSquaresBatchedWorkspace workspace(A); + const size_t allocated = session.totalAllocatedBytes(); + EXPECT_EQ(allocated - before, batch * sizeof(int) + (batch == 1 ? 2 * sizeof(T *) : 0)); + for (int repeat = 0; repeat < 3; ++repeat) { + DTensor &lhs = repeat % 2 ? otherA : A; + DTensor &result = repeat % 2 ? otherB : b; + lhs.upload(data); result.upload(rhs); + lhs.leastSquaresBatched(result, workspace); + std::vector solution; + result.download(solution); + for (size_t k = 0; k < batch; ++k) { + ASSERT_EQ(workspace.info()(k), 0); + EXPECT_NEAR(solution[k * rows], T(1), epsilon); + EXPECT_NEAR(solution[k * rows + 1], T(2), epsilon); + } + EXPECT_EQ(session.totalAllocatedBytes(), allocated); + } + } + EXPECT_EQ(session.totalAllocatedBytes(), before); + if (batch == 1) { + EXPECT_EQ(A.ptrMatrices(), nullptr); + EXPECT_EQ(b.ptrMatrices(), nullptr); + DTensor status(1); + EXPECT_THROW(A.leastSquaresBatched(b, status), std::invalid_argument); + } + A.upload(data); b.upload(rhs); + A.leastSquaresBatched(b); // Convenience overload constructs its own temporary workspace. + std::vector solution; + b.download(solution); + EXPECT_NEAR(solution[0], T(1), epsilon); + EXPECT_NEAR(solution[1], T(2), epsilon); + } + } +} + +TEST_F(LeastSquaresTest, singleMatrixBatchFloat) { + leastSquaresSingleMatrixBatch(PRECISION_LOW); +} + +TEST_F(LeastSquaresTest, singleMatrixBatchDouble) { + leastSquaresSingleMatrixBatch(PRECISION_HIGH); +} + +TEST_F(LeastSquaresTest, workspaceRejectsIncompatibleShapes) { + DTensor A(3, 2), b(3), differentRows(4, 2), differentCols(3, 1), differentBatch(3, 2, 2); + DTensor shortB(2), matrixB(3, 2), batchedB(3, 1, 2), wide(2, 3), empty; + const size_t before = Session::getInstance().totalAllocatedBytes(); + EXPECT_THROW(LeastSquaresBatchedWorkspace workspace(wide), std::invalid_argument); + EXPECT_THROW(LeastSquaresBatchedWorkspace workspace(empty), std::invalid_argument); + EXPECT_EQ(Session::getInstance().totalAllocatedBytes(), before); + LeastSquaresBatchedWorkspace workspace(A); + EXPECT_THROW(differentRows.leastSquaresBatched(b, workspace), std::invalid_argument); + EXPECT_THROW(differentCols.leastSquaresBatched(b, workspace), std::invalid_argument); + EXPECT_THROW(differentBatch.leastSquaresBatched(b, workspace), std::invalid_argument); + EXPECT_THROW(A.leastSquaresBatched(shortB, workspace), std::invalid_argument); + EXPECT_THROW(A.leastSquaresBatched(matrixB, workspace), std::invalid_argument); + EXPECT_THROW(A.leastSquaresBatched(batchedB, workspace), std::invalid_argument); +} + +template +void leastSquaresReusesStatusBuffer(T epsilon) { + auto &session = Session::getInstance(); + for (size_t rows : {size_t(2), size_t(3)}) { + const std::vector matrix = rows == 2 + ? std::vector{1,0,0,1, 2,0,0,2} + : std::vector{1,0,1,0,1,1, 2,0,2,0,2,2}; + DTensor A(matrix, rows, 2, 2); + DTensor B(rows, 1, 2); + DTensor status(3); // Excess capacity is allowed and left untouched. + DTensor tooSmall(1); + DTensor empty; + EXPECT_THROW(A.leastSquaresBatched(B, tooSmall), std::invalid_argument); + EXPECT_THROW(A.leastSquaresBatched(B, empty), std::invalid_argument); + const auto *statusPointer = status.raw(); + const size_t allocated = session.totalAllocatedBytes(); + for (int repeat = 0; repeat < 3; ++repeat) { + T x = T(2+repeat), y = T(3-repeat); + std::vector rhs = rows == 2 + ? std::vector{x,y,2*x,2*y} + : std::vector{x,y,x+y,2*x,2*y,2*(x+y)}; + A.upload(matrix); // gels overwrites its inputs, so restore each run. + B.upload(rhs); + status.upload(std::vector{-1,-1,123}); + A.leastSquaresBatched(B, status); + session.synchronizeStream(A.streamIdx()); + EXPECT_EQ(allocated, session.totalAllocatedBytes()); + EXPECT_EQ(statusPointer, status.raw()); + std::vector solution; + B.download(solution); + for (size_t batch = 0; batch < 2; ++batch) { + EXPECT_NEAR(x, solution[batch*rows], epsilon); + EXPECT_NEAR(y, solution[batch*rows+1], epsilon); + } + std::vector codes; + status.download(codes); + EXPECT_EQ((std::vector{0,0,123}), codes); + } + } +} + +TEST_F(LeastSquaresTest, reusableStatusBuffer) { + leastSquaresReusesStatusBuffer(PRECISION_LOW); + leastSquaresReusesStatusBuffer(PRECISION_HIGH); +} + /* ================================================================================================ * SVD TESTS @@ -1051,6 +2029,74 @@ protected: } }; +template +void svdReleasesOwnedCopyOnConstructionFailure() { + auto &session = Session::getInstance(); + const std::vector data{1,0,0,0,2,0, 3,0,0,0,4,0}; + DTensor source(data,3,2,2); + auto *sourcePointer = source.raw(); + int workspaceElements = 0; + if constexpr (std::is_same_v) { + gpuErrChk(cusolverDnSgesvd_bufferSize(session.cuSolverHandle(),3,2,&workspaceElements)); + } else { + gpuErrChk(cusolverDnDgesvd_bufferSize(session.cuSolverHandle(),3,2,&workspaceElements)); + } + ASSERT_GT(workspaceElements,0); + { Svd warmup(source,true,false); warmup.factorise(); } + const size_t copyBytes = data.size()*sizeof(T)+2*sizeof(T *); + const size_t infoBytes = 2*sizeof(int)+2*sizeof(int *); + const size_t outputBytes = 8*sizeof(T)+2*sizeof(T *); + for (bool destroyMatrix : {false,true}) { + SCOPED_TRACE(destroyMatrix ? "borrowed input" : "owned copy"); + const size_t baseline = session.totalAllocatedBytes(); + failureDeviceBytes = baseline+infoBytes+(destroyMatrix ? 0 : copyBytes) + +size_t(workspaceElements)*sizeof(T)+outputBytes; + failureHostBytes = 2*sizeof(T *); + observedDeviceBytes = observedHostBytes = 0; + bool caught = false; + { + gputils_test::FailHostAllocation failure(failAfterDeviceAllocation); + try { Svd broken(source,true,destroyMatrix); } + catch (const std::bad_alloc &) { caught = true; } + } + EXPECT_TRUE(caught); + // Failure occurs in V' pointer-table setup, after the optional input + // copy and workspace have completed construction and allocated on GPU. + EXPECT_EQ(failureDeviceBytes,observedDeviceBytes); + EXPECT_EQ(failureHostBytes,observedHostBytes); + EXPECT_EQ(baseline,session.totalAllocatedBytes()); + EXPECT_EQ(sourcePointer,source.raw()); + std::vector actual; + source.download(actual); + EXPECT_EQ(data,actual); + { + Svd recovered(source,true,false); + recovered.factorise(); + EXPECT_EQ(0,recovered.info()(0,0,0)); + EXPECT_EQ(0,recovered.info()(0,0,1)); + } + EXPECT_EQ(baseline,session.totalAllocatedBytes()); + source.download(actual); + EXPECT_EQ(data,actual); + } + const size_t beforeBorrow = session.totalAllocatedBytes(); + { + Svd borrowed(source,false,true); + borrowed.factorise(); + } + EXPECT_EQ(beforeBorrow,session.totalAllocatedBytes()); + EXPECT_EQ(sourcePointer,source.raw()); + source.upload(data); // Successful borrowed-input destruction must not free source. + std::vector actual; + source.download(actual); + EXPECT_EQ(data,actual); +} + +TEST_F(SvdTest, constructionFailureReleasesOwnedCopy) { + svdReleasesOwnedCopyOnConstructionFailure(); + svdReleasesOwnedCopyOnConstructionFailure(); +} + /* --------------------------------------- * Computation of singular values * and matrix rank @@ -1198,6 +2244,168 @@ TEST_F(SvdTest, singularValuesRankMultipleMatrices) { singularValuesRankMultipleMatrices(PRECISION_HIGH); } +template +void rankRepeatedCallsRecomputeCounts() { + SCOPED_TRACE((std::is_same_v ? "float" : "double")); + // Diagonal matrices have known singular values, well away from thresholds. + const std::vector data{ + 5,0,0, 0,2,0, 0,0,T(0.01), + 4,0,0, 0,T(0.05),0, 0,0,0 + }; + DTensor matrices(data,3,3,2); + Svd svd(matrices); + svd.factorise(); + Session::getInstance().synchronizeAllStreams(); + EXPECT_EQ(0,svd.info()(0)); + EXPECT_EQ(0,svd.info()(0,0,1)); + + const std::vector thresholds{T(0.1),T(0.1),T(0.001),T(0.1),T(10)}; + const std::vector> expected{{2,1},{2,1},{3,2},{2,1},{0,0}}; + for (size_t call=0;call actual; + svd.rank(thresholds[call]).download(actual); + EXPECT_EQ(expected[call],actual); + } +} + +TEST_F(SvdTest, rankRepeatedCallsRecomputeCounts) { + rankRepeatedCallsRecomputeCounts(); + rankRepeatedCallsRecomputeCounts(); +} + +template +void rankCallsAreOrderedOnSelectedStream() { + NonblockingTestStream execution; + DTensor producer(std::vector{5,0,0,2, 4,0,0,T(0.05)},2,2,2); + DTensor matrices(2,2,2,true); + matrices.setStreamIdx(1); + Svd svd(matrices); + DTensor first(1,1,2), second(1,1,2), third(1,1,2); + gpuErrChk(cudaLaunchHostFunc(execution.stream,delayTestInput,nullptr)); + gpuErrChk(cudaMemcpyAsync(matrices.raw(),producer.raw(),8*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + svd.factorise(); + const auto &rank1 = svd.rank(T(0.1)); + EXPECT_EQ(1,rank1.streamIdx()); + gpuErrChk(cudaMemcpyAsync(first.raw(),rank1.raw(),2*sizeof(unsigned int),cudaMemcpyDeviceToDevice,execution.stream)); + const auto &rank2 = svd.rank(T(0.001)); + gpuErrChk(cudaMemcpyAsync(second.raw(),rank2.raw(),2*sizeof(unsigned int),cudaMemcpyDeviceToDevice,execution.stream)); + const auto &rank3 = svd.rank(T(10)); + gpuErrChk(cudaMemcpyAsync(third.raw(),rank3.raw(),2*sizeof(unsigned int),cudaMemcpyDeviceToDevice,execution.stream)); + // No synchronization between calls: reset, count and downstream copies + // must stay ordered on the same stream. + gpuErrChk(cudaStreamSynchronize(execution.stream)); + std::vector actual; + first.download(actual); + EXPECT_EQ((std::vector{2,1}),actual); + second.download(actual); + EXPECT_EQ((std::vector{2,2}),actual); + third.download(actual); + EXPECT_EQ((std::vector{0,0}),actual); +} + +TEST_F(SvdTest, rankCallsUseSelectedStream) { + rankCallsAreOrderedOnSelectedStream(); + rankCallsAreOrderedOnSelectedStream(); +} + +template +void relativeRankRespectsPerMatrixScale() { + const std::vector scales{std::ldexp(T(1),-20),std::ldexp(T(1),20),std::ldexp(T(1),-10),T(1)}; + std::vector data(36,T(0)); + for (size_t batch=0;batch<3;++batch) { + data[9*batch] = T(8)*scales[batch]; + data[9*batch+4] = (batch == 2 ? T(0.25) : T(2))*scales[batch]; + data[9*batch+8] = (batch == 2 ? T(0) : T(0.125))*scales[batch]; + } + DTensor matrices(data,3,3,4); + Svd svd(matrices); + svd.factorise(); + const size_t allocated = Session::getInstance().totalAllocatedBytes(); + std::vector actual; + for (int repeat=0;repeat<2;++repeat) { + svd.rankRelative(T(0.1)).download(actual); + EXPECT_EQ((std::vector{2,2,1,0}),actual); + svd.rankRelative(T(0.1),T(1)).download(actual); + EXPECT_EQ((std::vector{0,2,0,0}),actual); + svd.rankRelative(T(0),T(0)).download(actual); + EXPECT_EQ((std::vector{3,3,2,0}),actual); + // The old absolute API intentionally remains sensitive to scaling. + svd.rank(T(0.1)).download(actual); + EXPECT_EQ((std::vector{0,3,0,0}),actual); + } + svd.rankRelative(T(0.1)).download(actual); + EXPECT_EQ((std::vector{2,2,1,0}),actual); + EXPECT_EQ(allocated,Session::getInstance().totalAllocatedBytes()); +} + +TEST_F(SvdTest, relativeRankUsesPerMatrixScaleAndAbsoluteCutoff) { + relativeRankRespectsPerMatrixScale(); + relativeRankRespectsPerMatrixScale(); +} + +template +void relativeRankChecksBoundariesAndTolerances() { + DTensor matrix(std::vector{8,0,0,0,2,0,0,0,T(0.125)},3,3); + Svd svd(matrix); + svd.factorise(); + std::vector actual; + svd.rankRelative(T(0.25)).download(actual); // threshold=2: equality excluded + EXPECT_EQ((std::vector{1}),actual); + svd.rankRelative(T(0),T(0.125)).download(actual); + EXPECT_EQ((std::vector{2}),actual); + svd.rankRelative(T(1)).download(actual); // threshold=sigmaMax + EXPECT_EQ((std::vector{0}),actual); + svd.rankRelative(T(0),std::numeric_limits::max()).download(actual); + EXPECT_EQ((std::vector{0}),actual); + const auto &previous = svd.rankRelative(T(0.25)); + for (T invalid : {T(-1),std::numeric_limits::infinity(),std::numeric_limits::quiet_NaN()}) { + EXPECT_THROW(svd.rankRelative(invalid),std::invalid_argument); + EXPECT_THROW(svd.rankRelative(T(0),invalid),std::invalid_argument); + } + previous.download(actual); + EXPECT_EQ((std::vector{1}),actual); // Invalid calls do not clear results. +} + +TEST_F(SvdTest, relativeRankValidatesTolerancesAndUsesStrictCutoff) { + relativeRankChecksBoundariesAndTolerances(); + relativeRankChecksBoundariesAndTolerances(); +} + +template +void relativeRankUsesSelectedStream() { + NonblockingTestStream execution; + DTensor producer(std::vector{8,0,0,2,4,0,0,T(0.125)},2,2,2); + DTensor matrix(2,2,2,true); + matrix.setStreamIdx(1); + Svd svd(matrix); + DTensor first(2),second(2),third(2); + gpuErrChk(cudaLaunchHostFunc(execution.stream,delayTestInput,nullptr)); + gpuErrChk(cudaMemcpyAsync(matrix.raw(),producer.raw(),8*sizeof(T),cudaMemcpyDeviceToDevice,execution.stream)); + svd.factorise(); + const auto &rank1 = svd.rankRelative(T(0.25)); + EXPECT_EQ(1,rank1.streamIdx()); + gpuErrChk(cudaMemcpyAsync(first.raw(),rank1.raw(),2*sizeof(unsigned int),cudaMemcpyDeviceToDevice,execution.stream)); + const auto &rank2 = svd.rankRelative(T(0.01)); + gpuErrChk(cudaMemcpyAsync(second.raw(),rank2.raw(),2*sizeof(unsigned int),cudaMemcpyDeviceToDevice,execution.stream)); + const auto &rank3 = svd.rank(T(0.2)); + gpuErrChk(cudaMemcpyAsync(third.raw(),rank3.raw(),2*sizeof(unsigned int),cudaMemcpyDeviceToDevice,execution.stream)); + gpuErrChk(cudaStreamSynchronize(execution.stream)); + std::vector actual; + first.download(actual); + EXPECT_EQ((std::vector{1,1}),actual); + second.download(actual); + EXPECT_EQ((std::vector{2,2}),actual); + third.download(actual); + EXPECT_EQ((std::vector{2,1}),actual); +} + +TEST_F(SvdTest, relativeRankOrdersRepeatedCallsOnSelectedStream) { + relativeRankUsesSelectedStream(); + relativeRankUsesSelectedStream(); +} + /* ================================================================================================ * CHOLESKY TESTS * ================================================================================================ */ @@ -1210,6 +2418,52 @@ protected: } }; +template +void choleskySingleMatrixBatch(T epsilon) { + NonblockingTestStream execution; + auto &session = Session::getInstance(); + for (bool factorised : {false, true}) { + DTensor A(factorised ? std::vector{2, 0, 0, 3} : std::vector{4, 0, 0, 9}, 2, 2, 1); + DTensor b(2), otherB(2); + A.setStreamIdx(1); b.setStreamIdx(1); otherB.setStreamIdx(1); + ASSERT_EQ(A.ptrMatrices(), nullptr); + ASSERT_EQ(b.ptrMatrices(), nullptr); + const size_t before = session.totalAllocatedBytes(); + { + CholeskyBatchFactoriser chol(A, factorised); + const size_t allocated = session.totalAllocatedBytes(); + EXPECT_EQ(allocated - before, sizeof(int) + 2 * sizeof(T *)); + if (!factorised) { + chol.factorise(); + ASSERT_EQ(chol.info()(0), 0); + } + for (int repeat = 0; repeat < 3; ++repeat) { + DTensor &rhs = repeat % 2 ? otherB : b; + rhs.upload(std::vector{4, 18}); + chol.solve(rhs); + std::vector solution; + rhs.download(solution); + ASSERT_EQ(chol.info()(0), 0); + ASSERT_EQ(solution.size(), 2u); + EXPECT_NEAR(solution[0], T(1), epsilon); + EXPECT_NEAR(solution[1], T(2), epsilon); + EXPECT_EQ(session.totalAllocatedBytes(), allocated); + } + } + EXPECT_EQ(session.totalAllocatedBytes(), before); + EXPECT_EQ(A.ptrMatrices(), nullptr); + EXPECT_EQ(b.ptrMatrices(), nullptr); + } +} + +TEST_F(CholeskyTest, singleMatrixBatchFloat) { + choleskySingleMatrixBatch(PRECISION_LOW); +} + +TEST_F(CholeskyTest, singleMatrixBatchDouble) { + choleskySingleMatrixBatch(PRECISION_HIGH); +} + /* --------------------------------------- * Cholesky factorisation @@ -1457,6 +2711,64 @@ TEST_F(QRTest, qrFactorisation) { qrFactorisation(PRECISION_HIGH); } +// These shapes exercise ormqr workspace requirements exceeding geqrf's on CUDA 13. +// Check both callers of ormqr, reusing the same factorisation and workspace. +template +void qrLargeWorkspace(size_t m, size_t n, T tolerance) { + std::vector data(m * n, T(0)); + std::vector expected(n), rhs(m, T(0)); + for (size_t c = 0; c < n; ++c) { + expected[c] = c % 2 ? T(-1) : T(1); + for (size_t r = c; r < m; r += n) { + data[r + c * m] = T(2); + rhs[r] += T(2) * expected[c]; + if (r + 1 < m) { + data[r + 1 + c * m] = T(0.25); + rhs[r + 1] += T(0.25) * expected[c]; + } + } + } + DTensor original(data, m, n, 1, columnMajor); + DTensor factors(original); + factors.setStreamIdx(1); + QRFactoriser qr(factors); + qr.factorise(); + ASSERT_EQ(qr.info()(0), 0); + + DTensor b(rhs, m); + b.setStreamIdx(1); + qr.leastSquares(b); + ASSERT_EQ(qr.info()(0), 0); + std::vector solution; + b.download(solution); + for (size_t c = 0; c < n; ++c) + ASSERT_NEAR(solution[c], expected[c], tolerance) << "column " << c; + + DTensor Q(m, n), R(n, n), reconstructed(m, n); + Q.setStreamIdx(1); R.setStreamIdx(1); reconstructed.setStreamIdx(1); + qr.getQR(Q, R); + ASSERT_EQ(qr.info()(0), 0); + reconstructed.addAB(Q, R); + reconstructed -= original; + EXPECT_LT(reconstructed.normF() / original.normF(), tolerance); +} + +TEST_F(QRTest, largeSquareWorkspaceFloat) { + qrLargeWorkspace(512, 512, 1e-4f); +} + +TEST_F(QRTest, largeSquareWorkspaceDouble) { + qrLargeWorkspace(512, 512, 1e-10); +} + +TEST_F(QRTest, largeTallWorkspaceFloat) { + qrLargeWorkspace(4096, 256, 1e-4f); +} + +TEST_F(QRTest, largeTallWorkspaceDouble) { + qrLargeWorkspace(4096, 256, 1e-10); +} + /* --------------------------------------- * QR factorisation * - tall and skinny matrix @@ -1662,6 +2974,157 @@ TEST_F(NullspaceTest, projectOnNullspaceTensor) { projectOnNullspaceTensor(PRECISION_HIGH); } +template +void denseNullspaceProjectionMatchesCpu(size_t n, size_t batches, T tolerance) { + // For A = v^T with nonzero column vector v, ker(A) has the dense projector + // P = I - v*v^T/(v^T*v). The CPU reference needs no GPU factorisation. + std::vector rows(n * batches), input(n * batches); + for (size_t k = 0; k < batches; ++k) { + for (size_t i = 0; i < n; ++i) { + const size_t index = i + k * n; + rows[index] = T(1) + T((i + 2 * k) % 7) / T(8); + if (k % 2 && i % 2) rows[index] = -rows[index]; + input[index] = T(16) + T(int(i % 17) - 8) / T(4) + T(k); + } + } + DTensor A(rows, 1, n, batches); + Nullspace nullspace(A); + + // Independent out-of-place control uses the same N*N^T construction as + // Nullspace, but keeps GEMM input/output storage separate. + auto transpose = nullspace.nullspace().tr(); + DTensor projector(n, n, batches); + projector.addAB(nullspace.nullspace(), transpose); + DTensor source(n, 1, batches), control(n, 1, batches), projected(n, 1, batches); + const auto maxError = [](const std::vector &actual, const std::vector &expected) { + T error = 0; + for (size_t i = 0; i < actual.size(); ++i) { + if (!std::isfinite(actual[i])) return std::numeric_limits::infinity(); + error = std::max(error, std::abs(actual[i] - expected[i])); + } + return error; + }; + + for (int repeat = 0; repeat < 3; ++repeat) { + auto current = input; + for (size_t i = 0; i < current.size(); ++i) + current[i] += T(repeat) * T(int(i % 5) - 2) / T(2); + std::vector expected(current.size()); + for (size_t k = 0; k < batches; ++k) { + long double numerator = 0, denominator = 0; + for (size_t i = 0; i < n; ++i) { + const size_t index = i + k * n; + numerator += static_cast(rows[index]) * current[index]; + denominator += static_cast(rows[index]) * rows[index]; + } + for (size_t i = 0; i < n; ++i) { + const size_t index = i + k * n; + expected[index] = T(static_cast(current[index]) + - static_cast(rows[index]) * numerator / denominator); + } + } + source.upload(current); + control.addAB(projector, source); + std::vector actual; + control.download(actual); + ASSERT_EQ(actual.size(), expected.size()); + ASSERT_LE(maxError(actual, expected), tolerance) + << "Out-of-place control must match the CPU reference before testing project()"; + + projected.upload(current); + for (int application = 0; application < 2; ++application) { + // The second application checks P(Pb) = Pb, without restoring b. + nullspace.project(projected); + projected.download(actual); + ASSERT_EQ(actual.size(), expected.size()); + EXPECT_LE(maxError(actual, expected), tolerance) + << "n=" << n << ", batches=" << batches << ", input=" << repeat + << ", projection=" << application; + } + } +} + +// Larger dense projectors exercise GEMM tiling that small vector tests can miss. +// Float tolerances allow basis-construction roundoff at these sizes; the control +// checks that this error stays well below the observed order-one aliasing errors. +TEST_F(NullspaceTest, denseProjectionRepeatedFloat) { + denseNullspaceProjectionMatchesCpu(4096, 1, 1e-2f); +} + +TEST_F(NullspaceTest, denseProjectionRepeatedDouble) { + denseNullspaceProjectionMatchesCpu(4096, 1, 1e-10); +} + +TEST_F(NullspaceTest, batchedDenseProjectionRepeatedFloat) { + denseNullspaceProjectionMatchesCpu(2048, 3, 1e-2f); +} + +TEST_F(NullspaceTest, batchedDenseProjectionRepeatedDouble) { + denseNullspaceProjectionMatchesCpu(2048, 3, 1e-10); +} + +template +void nullspaceProjectionReusesScratchOnSelectedStream(T tolerance) { + NonblockingTestStream setup(0); + NonblockingTestStream execution; + auto &session = Session::getInstance(); + for (size_t batches : {size_t(1), size_t(3)}) { + const size_t n = 4; + std::vector input(n * batches), expected(n * batches); + for (size_t k = 0; k < batches; ++k) + for (size_t i = 0; i < n; ++i) { + input[i + k * n] = T(k + 1) * T(i + 1); + expected[i + k * n] = T(k + 1) * (T(i) - T(1.5)); + } + DTensor A(std::vector(n * batches, T(1)), 1, n, batches); + DTensor producer(input, n, 1, batches), vectors(n, 1, batches, true), captured(n, 1, batches); + vectors.setStreamIdx(1); captured.setStreamIdx(1); + T *originalPointer = vectors.raw(); + const size_t before = session.totalAllocatedBytes(); + { + Nullspace nullspace(A); // Construct on stream 0, consume on nonblocking stream 1. + const size_t allocated = session.totalAllocatedBytes(); + for (int repeat = 0; repeat < 3; ++repeat) { + gpuErrChk(cudaLaunchHostFunc(execution.stream, delayTestInput, nullptr)); + gpuErrChk(cudaMemcpyAsync(vectors.raw(), producer.raw(), input.size() * sizeof(T), + cudaMemcpyDeviceToDevice, execution.stream)); + nullspace.project(vectors); + nullspace.project(vectors); // Reuse scratch without a host synchronization. + gpuErrChk(cudaMemcpyAsync(captured.raw(), vectors.raw(), input.size() * sizeof(T), + cudaMemcpyDeviceToDevice, execution.stream)); + std::vector actual; + captured.download(actual); + ASSERT_EQ(actual.size(), expected.size()); + for (size_t i = 0; i < actual.size(); ++i) EXPECT_NEAR(actual[i], expected[i], tolerance); + EXPECT_EQ(vectors.raw(), originalPointer); + EXPECT_EQ(vectors.streamIdx(), 1u); + EXPECT_EQ(session.totalAllocatedBytes(), allocated); + } + } + EXPECT_EQ(session.totalAllocatedBytes(), before); + } +} + +TEST_F(NullspaceTest, projectionReusesScratchOnSelectedStreamFloat) { + nullspaceProjectionReusesScratchOnSelectedStream(PRECISION_LOW); +} + +TEST_F(NullspaceTest, projectionReusesScratchOnSelectedStreamDouble) { + nullspaceProjectionReusesScratchOnSelectedStream(PRECISION_HIGH); +} + +TEST_F(NullspaceTest, projectionRejectsIncompatibleShapes) { + DTensor A(std::vector{1,1,1}, 1, 3); + Nullspace nullspace(A); + DTensor shortVector(2), matrix(3, 2), wrongBatch(3, 1, 2), empty; + const size_t before = Session::getInstance().totalAllocatedBytes(); + EXPECT_THROW(nullspace.project(shortVector), std::invalid_argument); + EXPECT_THROW(nullspace.project(matrix), std::invalid_argument); + EXPECT_THROW(nullspace.project(wrongBatch), std::invalid_argument); + EXPECT_THROW(nullspace.project(empty), std::invalid_argument); + EXPECT_EQ(Session::getInstance().totalAllocatedBytes(), before); +} + /* ================================================================================================ * GIVENSANNIHILATOR TESTS @@ -1705,6 +3168,30 @@ TEST_F(GivensAnnihilatorTest, givensAnnihilateElement) { givensAnnihilateElement(PRECISION_HIGH); } +template +void givensZeroPairLeavesMatrixUnchanged() { + // Column-major [[0,2], [0,-3], [7,11]]: the selected pair is already zero. + // An identity rotation must also preserve nonzero entries in the other column. + const std::vector expected{0, 0, 7, 2, -3, 11}; + DTensor matrix(expected, 3, 2); + GivensAnnihilator annihilator(matrix); + annihilator.annihilate(0, 1, 0); + std::vector actual; + matrix.download(actual); + ASSERT_EQ(actual.size(), expected.size()); + for (size_t index = 0; index < actual.size(); ++index) + ASSERT_TRUE(std::isfinite(actual[index])) << "Non-finite element at flat index " << index; + EXPECT_EQ(actual, expected); +} + +TEST_F(GivensAnnihilatorTest, zeroPairLeavesMatrixUnchangedFloat) { + givensZeroPairLeavesMatrixUnchanged(); +} + +TEST_F(GivensAnnihilatorTest, zeroPairLeavesMatrixUnchangedDouble) { + givensZeroPairLeavesMatrixUnchanged(); +} + /* --------------------------------------- * GivensAnnihilator: correctness