GPUtils

1. DTensor

The DTensor class is for manipulating data on a GPU. It manages their memory and facilitates various algebraic operations.

A tensor has three axes: [rows (m) x columns (n) x matrices (k)]. An (m,n,1)-tensor stores a matrix, and an (m,1,1)-tensor stores a vector.

Explicit tensor dimensions must be positive and fit in the integer dimensions used by CUDA's BLAS and solver interfaces. A default-constructed tensor is an uninitialized placeholder. Tensor sizes and allocation sizes are checked for overflow. Slices are non-owning views: their source storage must remain alive and must not be reallocated while the view is used.

Allocating constructors and the deep-copy constructor clean up acquired device storage if a C++ exception occurs during construction, then rethrow the original exception. For example, a host conversion-buffer allocation failure propagates std::bad_alloc without retaining the constructor's GPU allocations. This does not change gpuErrChk's existing process-termination policy for CUDA errors.

Numerical operations intended for repeated use may leave operand and index validation to the caller. In particular, addAB requires compatible matrix and batch dimensions and correctly sized output, += and -= require matching tensor shapes, and Givens rotation indices must be in range. Validate these preconditions during application setup before entering time-critical code.

We first need to decide on a data type between float or double. We will use float in the following examples.

1.1. Vectors

The simplest way to create an empty DTensor object is by constructing a vector:

size_t n = 100;
DTensor myTensor(n);

Important

This creates an n-dimensional vector as an (n,1,1)-tensor on the device.

A DTensor can be instantiated from host memory:

std::vector<float> h_a{4., -5., 6., 9., 8., 5., 9., -10.2, 9., 11.};
DTensor<float> myTensor(h_a, h_a.size());
std::cout << myTensor << "\n";

Caution

Printing a DTensor to std::cout will slow down your program (it requires the data to be downloaded from the device). Printing was designed for quick debugging.

We will often need to create slices (or shallow copies) of a DTensor given a range of values. We can then do:

size_t axis = 0;  // rows=0, cols=1, mats=2
size_t from = 3;
size_t to = 5;
DTensor<float> mySlice(myTensor, axis, from, to);
std::cout << mySlice << "\n";

Sometimes we need to reuse an already allocated DTensor by uploading new data from the host by using the method upload. Here is a short example:

std::vector<float> h_a{1., 2., 3.};  // host data a
DTensor<float> myVec(h_a, 3);  // create vector in tensor on device
std::vector<float> h_b{4., -5., 6.};  // host data b
myVec.upload(h_b);
std::cout << myVec << "\n";

We can upload some host data to a particular position of a DTensor as follows:

std::vector<float> hostData{1., 2., 3.};
// here, `true` tells the constructor to set all allocated elements to zero
DTensor<float> x(7, 1, 1, true);  // x = [0, 0, 0, 0, 0, 0, 0]'
DTensor<float> mySlice(x, 0, 3, 5); 
mySlice.upload(hostData);
std::cout << x << "\n";  // x = [0, 0, 0, 1, 2, 3, 0]'

If necessary, the data can be downloaded from the device to the host using download.

Very often we will also need to copy data from an existing DTensor to another DTensor (without passing through the host). To do this we can use deviceCopyTo. Here is an example:

DTensor<float> x(10);
DTensor<float> y(10);
x.deviceCopyTo(y);  // x ---> y (device memory to device memory)

For copies between preallocated buffers on the same stream, use deviceCopyToAsync to enqueue the copy without an explicit synchronization:

DTensor<float> source(std::vector<float>(10, 1.f), 10);
DTensor<float> destination(10);
source.deviceCopyToAsync(destination);
destination *= 2.f; // Ordered after the copy on the same selected stream.
std::vector<float> result;
destination.download(result); // Completes the stream before the CPU uses result.

An explicit-stream overload is also available:

cudaStream_t stream;
gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(source.streamIdx()), &stream));
source.deviceCopyToAsync(destination, stream);

For a nonempty source, both tensors must select the same actual CUDA stream; the explicit overload requires it to match stream too. Mismatches throw std::invalid_argument. The operation checks destination capacity, copies only the source-sized prefix, and leaves any destination tail unchanged. Identical data pointers are a no-op; partially overlapping copy regions are rejected. No storage is allocated or resized, and neither stream index changes.

The backing allocations must remain valid until the copy completes. Queue source producers and destination consumers on the same stream, or establish CUDA-event dependencies for accesses from other streams. Changing a stream index does not itself order existing work. Synchronize before the CPU needs completion or before releasing/replacing storage that is still in use. Deferred CUDA errors are reported at the relevant completion boundary. deviceCopyTo, copy construction, and copy assignment retain their existing completion guarantees.

For asynchronous host transfers, RegisteredHostVector<T> owns a fixed-size std::vector<T> and registers its existing allocation with CUDA. Construct it by moving an existing vector (no element copy), or by specifying an element count for a value-initialized buffer. Registration happens during construction; reuse the buffers for repeated transfers.

std::vector<float> values{1, 2, 3, 4}; // Column-major matrix data.
RegisteredHostVector<float> input(std::move(values));
RegisteredHostVector<float> output(4);
DTensor<float> matrix(2, 2);

cudaStream_t stream;
gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(matrix.streamIdx()), &stream));

matrix.uploadAsync(input);
matrix *= 2.f; // Ordered after the upload on the same stream.
matrix.downloadAsync(output);

// Independent CPU work can run here.
gpuErrChk(cudaStreamSynchronize(stream));
std::cout << output[0] << "\n"; // 2; output is now ready for CPU access.

Both methods also accept an explicit stream, for example matrix.downloadAsync(output, stream). It must match the tensor's actual selected CUDA stream. Host buffers must contain exactly matrix.numEl() elements; size and stream mismatches throw std::invalid_argument. Empty tensor/buffer pairs are a no-op. Transfers preserve shape, storage, and stream indices and perform no allocation or explicit synchronization.

These methods use column-major storage, including consecutive column-major matrices for batches. They do not convert row-major input; prepare the layout before uploading, or use the existing blocking upload(values, rowMajor).

Keep both host and device allocations alive until their transfers complete. Do not modify upload inputs or read/modify download outputs while a transfer is pending. Use a CUDA event or synchronize the actual selected stream before accessing/reusing buffers; accesses on other streams require explicit dependencies. The registered vector exposes data(), size(), empty(), and element access, but no resizing operations. It is movable and noncopyable. Move construction preserves its allocation and registration; wait before move-assigning over or destroying a buffer involved in a pending transfer. Destruction unregisters the memory before releasing the vector and introduces no explicit synchronization. Registration failure throws std::runtime_error and releases the adopted vector; registration requires CUDA platform support.

Copy construction and copy assignment create independent owning storage:

DTensor<float> 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 one element of the DTensor to the host. This cannot be used to set a value, so don't do anything like x(0, 0, 0) = 4.5!

Caution

For the love of god, do not put this () operator in a loop.

1.2. Computation of scalar quantities

The following scalar quantities can be computed (internally, we use cublas functions):

• .normF(): the Frobenius norm of a tensor $x$, using nrm2 (i.e., the 2-norm, or Euclidean norm, if $x$ is a vector)
• .sumAbs(): the sum of the absolute of all the elements, using asum (i.e., the 1-norm if $x$ is a vector)

1.3. Some cool operators

We can element-wise add DTensors on the device as follows:

std::vector<float> host_x{1., 2., 3., 4., 5., 6.,  7.};
std::vector<float> host_y{1., 3., 5., 7., 9., 11., 13.};
DTensor<float> x(host_x, host_x.size());
DTensor<float> y(host_y, host_y.size());
x += y;  // x = [2, 5, 8, 11, 14, 17, 20]'
std::cout << x << "\n";

To element-wise subtract y from x we can use x -= y.

We can also scale a DTensor by a scalar with *= (e.g, x *= 5.0f). To negate the values of a DTensor we can do x *= -1.0f.

We can also compute the inner product (as a (1,1,1)-tensor) of two vectors as follows:

std::vector<float> host_x{1., 2., 3., 4., 5., 6.,  7.};
std::vector<float> host_y{1., 3., 5., 7., 9., 11., 13.};
DTensor<float> xtr(host_x, 1, host_x.size());  // column vector
DTensor<float> y(host_y, host_y.size());  // row vector
DTensor<float> innerProduct = x * y;

If necessary, we can also use the following element-wise operations

DTensor<float> x(host_x, host_x.size());  // row vector
auto sum = x + y;
auto diff = x - y;
auto scaledX = 3.0f * x;

1.4. Matrices

To store a matrix in a DTensor we need to provide the data in an array; we can use either column-major (default) or row-major format. TODO implement row-major Suppose we need to store the matrix

$$A = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ 10 & 11 & 12 \\ 13 & 14 & 15 \end{bmatrix},$$

where this data is stored in row-major format. Then, we do

size_t rows = 5;
size_t cols = 3;
std::vector<float> h_data{1.0f, 2.0f, 3.0f,
                          4.0f, 5.0f, 6.0f,
                          7.0f, 8.0f, 9.0f,
                          10.0f, 11.0f, 12.0f,
                          13.0f, 14.0f, 15.0f};
DTensor<float> myTensor(h_data, rows, cols, 1, rowMajor);

Choose rowMajor or columnMajor as appropriate.

We can also preallocate memory for a DTensor as follows:

DTensor<float> a(rows, cols, 1);

Then, we can upload the data as follows:

a.upload(h_data, rowMajor);

The copy constructor has also been implemented; to hard-copy a vector just do DTensor<float> myCopy(existingTensor).

The number of rows and columns of a DTensor can be retrieved using the methods .numRows() and .numCols() respectively.

1.5. More operations

The operators += are -= supported for device matrices.

Matrix-matrix multiplication is as simple as:

size_t m = 2, k = 3, n=5;
std::vector<float> aData{1.0f,  2.0f,  3.0f,
                         4.0f,  5.0f,  6.0f};
std::vector<float> bData{1.0f,  2.0f,  3.0f,  4.0f,  5.0f,
                         6.0f,  7.0f,  8.0f,  9.0f, 10.0f,
                         11.0f, 12.0f, 13.0f, 14.0f, 15.0f};
DTensor<float> A(aData, m, k, 1, rowMajor);
DTensor<float> B(bData, k, n, 1, rowMajor);
auto X = A * B;
std::cout << A << B << X << "\n";

1.6. Tensors

As you would expect, all operations mentioned so far are supported by actual tensors as batched operations (that is, (m,n)-matrix-wise).

Also, we can create the transposes of a DTensor using .tr(). This transposes each (m,n)-matrix and stores it in a new DTensor at the same k-index. Transposition in-place is not possible.

1.7. Least squares

The solution of least squares has been implmented as a tensor method. Say we want to solve A\b using least squares. We first create $A$ and $b$

size_t m = 4;
size_t n = 3;
std::vector<float> aData{1.0f, 2.0f, 4.0f,
                         2.0f, 13.0f, 23.0f,
                         4.0f, 23.0f, 77.0f,
                         6.0f, 7.0f, 8.0f};
std::vector<float> bData{1.0f, 2.0f, 3.0f, 4.0f};
DTensor<float> A(aData, m, n, 1, rowMajor);
DTensor<float> B(bData, m);

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:

LeastSquaresBatchedWorkspace<float> workspace(A);
A.leastSquaresBatched(B, workspace);
Session::getInstance().synchronizeStream(A.streamIdx());
std::vector<int> 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:

DTensor<int> 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<int> 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<T> 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

This particular example demonstrates how the solution may overwrite only part of the given B, as B is a (4,1,1)-tensor and the solution is a (3,1,1)-tensor.

Reusable QR operations

QRFactoriser<T> supports nonempty tall or square single matrices in float and double. Select the matrix's stream before constructing the factoriser. The constructor allocates the reflectors and the maximum workspace required by factorisation, applying the reflectors to a vector, and forming thin Q.

DTensor<double> A(aData, m, n, 1, columnMajor); // m >= n
DTensor<double> rhs(bData, m);
DTensor<double> Q(m, n), R(n, n);
A.setStreamIdx(0);
QRFactoriser<double> qr(A);

qr.factorise();     // overwrites A with Householder reflectors and R
qr.applyQt(rhs);    // overwrites the m-entry vector with full Q^T b
qr.extractR(R);     // writes n-by-n R, including zeros below the diagonal
qr.formThinQ(Q);    // writes the m-by-n orthonormal columns of thin Q
Session::getInstance().synchronizeStream(A.streamIdx());

The first n entries returned by applyQt are thinQ^T b; the remaining m-n entries are coordinates in the orthogonal complement. It performs no triangular solve. leastSquares(rhs) instead applies the reflectors and solves the triangular system in place, requiring full column rank.

The three new operations preserve the stored reflectors and can be repeated in any order after factorise(). They perform no host transfers, GPUtils device allocations, or explicit synchronization. They always enqueue on A's selected cuSOLVER stream, irrespective of the outputs' streamIdx(), and leave output storage, shapes and stream indices unchanged. Order both prior producers and later consumers on other streams with explicit dependencies. Synchronize A's stream before reading outputs through a different stream, and retain all storage and the factoriser until the queued work completes. A workspace cannot be used concurrently on different streams. Select matching streams for a session's cuBLAS and cuSOLVER handles if configuring those handles directly.

Output shapes are checked and outputs must not overlap A. Refill A's existing storage to factorise another matrix; changing its shape, allocation or selected stream requires a new factoriser. getQR(Q, R) is a blocking convenience wrapper for the GPU operations and also rejects overlapping Q and R destinations. The new operations support CUDA Graph capture after library warm-up; construct the factoriser and output tensors before capture. info() reports the last queued cuSOLVER operation's status; extractR leaves it unchanged.

1.8. Saving and loading tensors

Tensor data can be stored in simple text files or binary files. The text-based format has the following structure

number_of_rows
number_of_columns
number_of_matrices
data (one entry per line)

To save a tensor in a file, simply call DTensor::saveToFile(filename).

If the file extension is .bt (binary tensor), the data will be stored in binary format. 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<T>::parseFromFile(filename) can be used. For example:

auto z = DTensor<double>::parseFromFile("path/to/my.dtensor")

If necessary, you can provide a second argument to parseFromFile to specify the order in which the data are stored (the StorageMode).

Soon we will release a Python API for reading and serialising (numpy) arrays to .bt files.

2. Cholesky factorisation and system solution

Warning

This factorisation only works with positive-definite matrices.

Here is an example:

$$A = \begin{bmatrix} 1 & 2 & 4 \\ 2 & 13 & 23 \\ 4 & 23 & 77 \end{bmatrix}.$$

This is how to perform a Cholesky factorisation:

size_t n = 3;
std::vector<float> aData{1.0f, 2.0f, 4.0f,
                         2.0f, 13.0f, 23.0f,
                         4.0f, 23.0f, 77.0f};
DTensor<float> A(aData, n, n, 1, rowMajor);
CholeskyFactoriser<float> cfEngine(A);
status = cfEngine.factorise();

Then, you can solve the system A\b

std::vector<float> bData{1.0f, 2.0f, 3.0f};
DTensor<float> B(bData, n);
cfEngine.solve(B);

The DTensor B will be overwritten with the solution.

3. Singular Value Decomposition

Warning

This implementation only works with square or tall matrices.

Here is an example with the 4-by-3 matrix

$$B = \begin{bmatrix} 1 & 2 & 3 \\ 6 & 7 & 8 \\ 6 & 7 & 8 \\ 6 & 7 & 8 \end{bmatrix}.$$

Evidently, the rank of $B$ is 2, so there will be two nonzero singular values.

This is how to perform an SVD decomposition:

size_t m = 4;
size_t n = 3;
std::vector<float> bData{1.0f, 2.0f, 3.0f,
                         6.0f, 7.0f, 8.0f,
                         6.0f, 7.0f, 8.0f,
                         6.0f, 7.0f, 8.0f};
DTensor<float> B(bData, m, n, 1, rowMajor);
SvdFactoriser<float> svdEngine(B);
status = svdEngine.factorise();

By default, SvdFactoriser will not compute matrix $U$. If you need it, create an instance of SvdFactoriser as follows

SvdFactoriser<float> svdEngine(B, true); // computes U

Note that the default behaviour of .factorise() is to destroy 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:

std::cout << "S = " << svdEngine.singularValues() << "\n";
std::cout << "V' = " << svdEngine.rightSingularVectors() << "\n";

Note that $U$ can be obtained, if it is computed in the first place, by the method .leftSingularVectors() which returns an object of type std::optional<DeviceMatrix<TElement>>. Here is an example:

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:

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<T>::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. The user provides a DTensor made of (padded) matrices. Then, Nullspace computes, possibly pads, and returns the nullspace matrices N = (N1, ..., Nk) in another DTensor.

DTensor<float> paddedMatrices(m, n, k);
Nullspace N(paddedMatrices);  // computes N and NN'
DTensor<float> ns = N.nullspace();  // returns N

Each padded nullspace matrix Ni is orthogonal, and Nullspace further computes and stores the nullspace projection operators NN' = (N1N1', ..., NkNk'). This allows the user to project-in-place onto the nullspace.

DTensor<float> 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

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

/* 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

auto a = DTensor<double>::createRandomTensor(3, 6, 4, -1, 1).setStreamIdx(0);
auto b = DTensor<double>::createRandomTensor(3, 6, 4, -1, 1).setStreamIdx(1);
// do stuff...
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. The asynchronous host-transfer methods below leave completion to the caller.
• 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.
• source.deviceCopyToAsync(destination[, stream]) enqueues on the common selected stream without explicit synchronization. It rejects nonempty copies with incompatible streams and partially overlapping regions. The caller retains the backing allocations and orders accesses from other streams.
• tensor.uploadAsync(host[, stream]) and tensor.downloadAsync(host[, stream]) transfer column-major data using an owning RegisteredHostVector of the exact element count. They enqueue on the tensor's selected stream; an explicit stream must match it. Keep allocations alive and wait before CPU access or buffer reuse.
• 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):

python3 -m venv build/test-venv
build/test-venv/bin/python -m pip install numpy

Run all C++ and Python unit tests:

./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:

./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!