From db3bfb592d3f8a895666670ec7ca552719f6111f Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 21:42:20 +0000 Subject: [PATCH 01/23] tests should all pass now --- test/testTensor.cu | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/test/testTensor.cu b/test/testTensor.cu index bff6ac8..691a107 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -1,6 +1,8 @@ #include #include "../include/tensor.cuh" +#include #include +#include #define PRECISION_LOW 1e-4 #define PRECISION_HIGH 1e-10 @@ -191,14 +193,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); } From d6e8f1ab4450e8fadd43def79d4e3b89cc9b15fc Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 21:48:48 +0000 Subject: [PATCH 02/23] updated .gitignore (cmake files etc) --- .gitignore | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) 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 From dd6200d047f293ae11814ed31e2152ee88fbcd2d Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 21:49:09 +0000 Subject: [PATCH 03/23] CMakeLists: enable_testing() --- CMakeLists.txt | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9b20dbc..e1efa37 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -67,6 +67,9 @@ if(NOT GPUTILS_BUILD_TEST) set(GPUTILS_BUILD_TEST OFF) # Set to ON for local testing (or add `-DGPUTILS_BUILD_TEST=ON` to your CMake profile) endif() if (GPUTILS_BUILD_TEST) + # Register tests in the top-level build tree so `ctest --test-dir build` + # discovers tests added by the test subdirectory. + enable_testing() add_subdirectory(test) endif() unset(GPUTILS_BUILD_TEST CACHE) From 6a771bf60970266ec4da9193410dec021e7cb9bb Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 22:05:18 +0000 Subject: [PATCH 04/23] more unit tests (corner cases) --- include/tensor.cuh | 192 ++++++++++++++++++++++++++++----------------- test/testTensor.cu | 170 ++++++++++++++++++++++++++++++++++++++- 2 files changed, 291 insertions(+), 71 deletions(-) diff --git a/include/tensor.cuh b/include/tensor.cuh index b73e280..e10c455 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -15,6 +15,11 @@ #include #include #include +#include +#include +#include +#include +#include /** @@ -131,12 +136,15 @@ static size_t s_numStreams = 1; * The cuSolver handle can be accessed anywhere by `Session::getInstance().cuSolverHandle()` */ class Session { + 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,6 +171,7 @@ private: gpuErrChk(cusolverDnCreate(&m_cusolverHandles[i])); gpuErrChk(cusolverDnSetStream(m_cusolverHandles[i], m_cublasStreams[i])); } + initialized() = true; } ~Session() { @@ -181,6 +190,8 @@ private: public: Session(Session const &) = delete; + size_t numStreams() const { return m_numCublasHandlesStreams; } + void operator=(Session const &) = delete; /** @@ -260,6 +271,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. @@ -549,7 +578,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 +589,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); @@ -662,7 +691,7 @@ public: template DTensor DTensor::setStreamIdx(size_t idx) { - if (idx >= s_numStreams) { + if (idx >= Session::getInstance().numStreams()) { throw std::invalid_argument("Invalid stream index; it exceeds the max allocated streams"); } m_idxStream = idx; @@ -689,12 +718,15 @@ void DTensor::initialisePointersToMatricesData() { 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 +749,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 +789,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 +835,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 +843,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 +867,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]; @@ -862,10 +901,10 @@ 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(); @@ -873,10 +912,12 @@ DTensor::DTensor(size_t m, size_t n, size_t k, bool zero) { 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 */ @@ -898,7 +939,9 @@ DTensor::DTensor(const DTensor &other) { 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 +949,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; @@ -1127,6 +1170,7 @@ 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 @@ -1277,6 +1321,8 @@ inline DTensor &DTensor::operator-=(const DTensor &rhs) template inline T DTensor::operator()(size_t i, size_t j, size_t k) const { + if (i >= m_numRows || j >= m_numCols || k >= m_numMats) + throw std::invalid_argument("tensor index out of range"); T hostDst; size_t offset = i + m_numRows * (j + m_numCols * k); gpuErrChk(cudaMemcpy(&hostDst, m_d_data + offset, sizeof(T), cudaMemcpyDeviceToHost)); @@ -1395,6 +1441,8 @@ inline void DTensor::leastSquaresBatched(DTensor &B) { template<> inline DTensor DTensor::getRows(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"); size_t rowsRangeLength = rowsTo - rowsFrom + 1; size_t n = numCols(), m = numRows(); DTensor rowsOnly(rowsRangeLength, numCols(), 1); @@ -1410,6 +1458,8 @@ inline DTensor DTensor::getRows(size_t rowsFrom, size_t rowsTo, template<> inline DTensor DTensor::getRows(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"); size_t rowsRangeLength = rowsTo - rowsFrom + 1; size_t n = numCols(), m = numRows(); DTensor rowsOnly(rowsRangeLength, numCols(), 1); @@ -1516,7 +1566,7 @@ private: * @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)"); } }; @@ -1701,6 +1751,7 @@ 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; @@ -1809,6 +1860,7 @@ 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; @@ -1930,9 +1982,9 @@ 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.); @@ -1964,9 +2016,9 @@ 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.); @@ -2045,7 +2097,7 @@ 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); auto aTranspose = a.tr(); @@ -2114,7 +2166,7 @@ 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_numRows = A.numRows(); m_numMats = A.numMats(); @@ -2234,7 +2286,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 +2300,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; diff --git a/test/testTensor.cu b/test/testTensor.cu index 691a107..8638b7b 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -25,6 +25,37 @@ 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 + * + * This test must remain before all tensor tests: Session is a singleton, so + * the number of streams has to be selected before its first use. + * --------------------------------------- */ + +TEST(SessionTest, streamSelectionAndExecution) { + EXPECT_THROW(Session::setStreams(0), std::invalid_argument); + Session::setStreams(2); + + 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); + Session::getInstance().synchronizeStream(1); + std::vector result; + C.download(result); + EXPECT_EQ((std::vector{23, 34, 31, 46}), result); +} + /* --------------------------------------- * Zero Tensor (Constructor) * --------------------------------------- */ @@ -359,7 +390,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()); @@ -1013,6 +1045,142 @@ TEST_F(TensorTest, tensorBytesAllocatedDeallocated) { tensorBytesAllocatedDeallocated(); } +/* --------------------------------------- + * 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); +} + /* ================================================================================================ From cf66da4786c940539b7dcced22761f2ca132002a Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 22:15:37 +0000 Subject: [PATCH 05/23] Avoid the host copy for column-major uploads. --- include/tensor.cuh | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/include/tensor.cuh b/include/tensor.cuh index e10c455..e9cf559 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -1175,15 +1175,12 @@ inline bool DTensor::upload(const std::vector &vec, StorageMode mode) { 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 (mode == StorageMode::rowMajor) { + std::vector vecCm(thisSize); rm2cm(vec, vecCm); + gpuErrChk(cudaMemcpy(m_d_data, vecCm.data(), size * sizeof(T), cudaMemcpyHostToDevice)); } else { - vecCm = vec; - } - if (size <= thisSize) { - size_t buffer_size = size * sizeof(T); - gpuErrChk(cudaMemcpy(m_d_data, vecCm.data(), buffer_size, cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpy(m_d_data, vec.data(), size * sizeof(T), cudaMemcpyHostToDevice)); } return true; } From a21a2481505e33e6ad26c009a307c973f790058a Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 22:35:27 +0000 Subject: [PATCH 06/23] benchmarks for getRows --- benchmark/CMakeLists.txt | 13 +++ benchmark/README.md | 73 ++++++++++++ benchmark/benchmark.cu | 238 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 324 insertions(+) create mode 100644 benchmark/CMakeLists.txt create mode 100644 benchmark/README.md create mode 100644 benchmark/benchmark.cu 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..60d1f06 --- /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`: current 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; } +} From 64abec7464fbae8c6851d0c7ebbe7f72de86d189 Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 22:47:57 +0000 Subject: [PATCH 07/23] Removed deep copy in setStreamIdx --- include/tensor.cuh | 45 ++++++++++++++++++++------------- test/testTensor.cu | 62 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 90 insertions(+), 17 deletions(-) diff --git a/include/tensor.cuh b/include/tensor.cuh index e9cf559..5149a05 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -20,6 +20,7 @@ #include #include #include +#include /** @@ -362,9 +363,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; } @@ -500,6 +503,9 @@ public: /** * 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) @@ -690,7 +696,7 @@ public: }; /* END OF DTENSOR */ template -DTensor DTensor::setStreamIdx(size_t idx) { +DTensor &DTensor::setStreamIdx(size_t idx) & { if (idx >= Session::getInstance().numStreams()) { throw std::invalid_argument("Invalid stream index; it exceeds the max allocated streams"); } @@ -698,6 +704,13 @@ DTensor DTensor::setStreamIdx(size_t 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 */ @@ -1443,13 +1456,12 @@ inline DTensor DTensor::getRows(size_t rowsFrom, size_t rowsTo, 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)); - } + rowsOnly.m_idxStream = m_idxStream; + cudaStream_t stream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(m_idxStream), &stream)); + gpuErrChk(cudaMemcpy2DAsync(rowsOnly.raw(), rowsRangeLength * sizeof(double), + raw() + rowsFrom + matIdx * n * m, m * sizeof(double), + rowsRangeLength * sizeof(double), n, cudaMemcpyDeviceToDevice, stream)); return rowsOnly; } @@ -1460,13 +1472,12 @@ inline DTensor DTensor::getRows(size_t rowsFrom, size_t rowsTo, si 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)); - } + rowsOnly.m_idxStream = m_idxStream; + cudaStream_t stream; + gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(m_idxStream), &stream)); + gpuErrChk(cudaMemcpy2DAsync(rowsOnly.raw(), rowsRangeLength * sizeof(float), + raw() + rowsFrom + matIdx * n * m, m * sizeof(float), + rowsRangeLength * sizeof(float), n, cudaMemcpyDeviceToDevice, stream)); return rowsOnly; } diff --git a/test/testTensor.cu b/test/testTensor.cu index 8638b7b..ca3ee75 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -50,12 +50,69 @@ TEST(SessionTest, streamSelectionAndExecution) { 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); } +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) * --------------------------------------- */ @@ -972,6 +1029,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) { From 0b268821ad557a5d65cbaeab0a546b2e493c7447 Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 22:51:56 +0000 Subject: [PATCH 08/23] Reusing least-squares workspace --- include/tensor.cuh | 27 +++++++++++++++++---- test/testTensor.cu | 59 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 82 insertions(+), 4 deletions(-) diff --git a/include/tensor.cuh b/include/tensor.cuh index 5149a05..e3724cb 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -610,6 +610,17 @@ 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. + */ + void leastSquaresBatched(DTensor &b, DTensor &infoArray); + /** * Batched `C <- bC + a*A*B`. * Performs the operation `Ci <- bCi + a*Ai*Bi` for each k-index `i`. @@ -1394,7 +1405,7 @@ inline void DTensor::addAB(const DTensor &A, const DTensor } template<> -inline void DTensor::leastSquaresBatched(DTensor &B) { +inline void DTensor::leastSquaresBatched(DTensor &B, DTensor &infoArray) { size_t batchSize = numMats(); size_t nColsB = B.numCols(); if (B.numRows() != m_numRows) @@ -1406,7 +1417,8 @@ 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, @@ -1422,7 +1434,7 @@ inline void DTensor::leastSquaresBatched(DTensor &B) { } template<> -inline void DTensor::leastSquaresBatched(DTensor &B) { +inline void DTensor::leastSquaresBatched(DTensor &B, DTensor &infoArray) { size_t batchSize = numMats(); size_t nColsB = B.numCols(); if (B.numRows() != m_numRows) @@ -1434,7 +1446,8 @@ 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, @@ -1449,6 +1462,12 @@ inline void DTensor::leastSquaresBatched(DTensor &B) { batchSize)); } +template +inline void DTensor::leastSquaresBatched(DTensor &B) { + DTensor infoArray(numMats()); + leastSquaresBatched(B, infoArray); +} + template<> inline DTensor DTensor::getRows(size_t rowsFrom, size_t rowsTo, size_t matIdx) const { if (rowsFrom > rowsTo || rowsTo >= numRows() || matIdx >= numMats()) diff --git a/test/testTensor.cu b/test/testTensor.cu index ca3ee75..b9b4ad6 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -59,6 +59,19 @@ TEST(SessionTest, streamSelectionAndExecution) { 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); } template @@ -1290,6 +1303,52 @@ TEST_F(LeastSquaresTest, tensorLS1) { tensorLeastSquares1(PRECISION_HIGH); } +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 From 12ac0630e2c23f5e04909c8b0eae4570ca36ca34 Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 23:08:05 +0000 Subject: [PATCH 09/23] Implemented getRowsInto --- README.md | 38 ++++++++++ benchmark/README.md | 2 +- include/tensor.cuh | 91 ++++++++++++++++-------- test/testTensor.cu | 169 +++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 265 insertions(+), 35 deletions(-) diff --git a/README.md b/README.md index 2728855..efe88af 100644 --- a/README.md +++ b/README.md @@ -8,6 +8,18 @@ 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. + +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. @@ -239,6 +251,25 @@ Then, we can solve the system by A.leastSquaresBatched(B); ``` +For repeated batched solves, allocate the per-matrix device status buffer once +during setup and pass it explicitly: + +```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 allocates its own +temporary status buffer. The supplied 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. + The `DTensor` `B` will be overwritten with the solution. > [!IMPORTANT] @@ -264,6 +295,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: @@ -409,6 +445,8 @@ 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. 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); diff --git a/benchmark/README.md b/benchmark/README.md index 60d1f06..70a373a 100644 --- a/benchmark/README.md +++ b/benchmark/README.md @@ -38,7 +38,7 @@ single column-major matrix. All options require positive integers. 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`: current per-row cuBLAS copy loop versus one asynchronous +- 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. diff --git a/include/tensor.cuh b/include/tensor.cuh index e3724cb..7b85fe6 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -513,6 +513,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. @@ -1468,19 +1478,33 @@ inline void DTensor::leastSquaresBatched(DTensor &B) { leastSquaresBatched(B, infoArray); } +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 { if (rowsFrom > rowsTo || rowsTo >= numRows() || matIdx >= numMats()) throw std::invalid_argument("invalid row or matrix range"); - size_t rowsRangeLength = rowsTo - rowsFrom + 1; - size_t n = numCols(), m = numRows(); - DTensor rowsOnly(rowsRangeLength, numCols(), 1); + DTensor rowsOnly(rowsTo - rowsFrom + 1, numCols(), 1); rowsOnly.m_idxStream = m_idxStream; - cudaStream_t stream; - gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(m_idxStream), &stream)); - gpuErrChk(cudaMemcpy2DAsync(rowsOnly.raw(), rowsRangeLength * sizeof(double), - raw() + rowsFrom + matIdx * n * m, m * sizeof(double), - rowsRangeLength * sizeof(double), n, cudaMemcpyDeviceToDevice, stream)); + getRowsInto(rowsOnly, rowsFrom, rowsTo, matIdx); return rowsOnly; } @@ -1488,15 +1512,9 @@ template<> inline DTensor DTensor::getRows(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"); - size_t rowsRangeLength = rowsTo - rowsFrom + 1; - size_t n = numCols(), m = numRows(); - DTensor rowsOnly(rowsRangeLength, numCols(), 1); + DTensor rowsOnly(rowsTo - rowsFrom + 1, numCols(), 1); rowsOnly.m_idxStream = m_idxStream; - cudaStream_t stream; - gpuErrChk(cublasGetStream(Session::getInstance().cuBlasHandle(m_idxStream), &stream)); - gpuErrChk(cudaMemcpy2DAsync(rowsOnly.raw(), rowsRangeLength * sizeof(float), - raw() + rowsFrom + matIdx * n * m, m * sizeof(float), - rowsRangeLength * sizeof(float), n, cudaMemcpyDeviceToDevice, stream)); + getRowsInto(rowsOnly, rowsFrom, rowsTo, matIdx); return rowsOnly; } @@ -1898,6 +1916,7 @@ public: /** * Factorise matrix. + * Work is enqueued on the bound matrix's selected stream. */ void factorise(); @@ -1906,6 +1925,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); @@ -1913,6 +1934,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) @@ -1926,7 +1948,7 @@ template<> inline void QRFactoriser::computeWorkspaceSize() { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); - gpuErrChk(cusolverDnDgeqrf_bufferSize(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnDgeqrf_bufferSize(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), m, n, nullptr, m, &m_workspaceSize)); @@ -1936,7 +1958,7 @@ template<> inline void QRFactoriser::computeWorkspaceSize() { size_t m = m_matrix->numRows(); size_t n = m_matrix->numCols(); - gpuErrChk(cusolverDnSgeqrf_bufferSize(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnSgeqrf_bufferSize(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), m, n, nullptr, m, &m_workspaceSize)); @@ -1946,7 +1968,7 @@ 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(), @@ -1959,7 +1981,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(), @@ -1972,7 +1994,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(), @@ -1991,7 +2013,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(), @@ -2022,13 +2044,15 @@ inline void QRFactoriser::getQR(DTensor &Q, DTensor &R) } Q.upload(vecQ, rowMajor); // 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++) { @@ -2056,13 +2080,15 @@ inline void QRFactoriser::getQR(DTensor &Q, DTensor &R) { } Q.upload(vecQ, rowMajor); // 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++) { @@ -2372,22 +2398,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/test/testTensor.cu b/test/testTensor.cu index b9b4ad6..168bf66 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -3,6 +3,8 @@ #include #include #include +#include +#include #define PRECISION_LOW 1e-4 #define PRECISION_HIGH 1e-10 @@ -28,13 +30,14 @@ protected: /* --------------------------------------- * Session streams * - * This test must remain before all tensor tests: Session is a singleton, so - * the number of streams has to be selected before its first use. + * 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 = [] { Session::setStreams(2); return true; }(); + TEST(SessionTest, streamSelectionAndExecution) { EXPECT_THROW(Session::setStreams(0), std::invalid_argument); - Session::setStreams(2); DTensor A({1, 2, 3, 4}, 2, 2); DTensor B({5, 6, 7, 8}, 2, 2); @@ -74,6 +77,103 @@ TEST(SessionTest, streamSelectionAndExecution) { 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 = Session::getInstance().cuBlasHandle(1); + cusolverDnHandle_t solver = Session::getInstance().cuSolverHandle(1); + cublasPointerMode_t oldMode; + NonblockingTestStream() { + 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 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 &>); @@ -1054,6 +1154,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 From a7d02147a56fa2ee2222df3356fa481a4d1c0593 Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 23:12:18 +0000 Subject: [PATCH 10/23] updated changelog --- CHANGELOG.md | 62 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5d853fe..345a6c2 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). + +## v1.10.0 - Unreleased + +### Added + +- `DTensor::getRowsInto(destination, rowsFrom, rowsTo, matIdx)` copies rows into a + preallocated tensor without allocating or resizing; validates indices, destination shape, + and non-overlapping storage, and enqueues the copy on the source tensor's stream. +- `DTensor::leastSquaresBatched(b, infoArray)` accepts a reusable device status buffer, + avoiding GPUtils's per-call workspace allocation while retaining the original overload. +- `Session::numStreams()` exposes the initialized stream count. +- Standalone Release-mode GPU benchmarks for pointer-based versus strided batched GEMM + and row-loop, 2D-copy, kernel, and actual `getRows` extraction paths. Includes configurable + dimensions, correctness checks, preallocated and allocation-inclusive measurements, + host/CUDA-event timings, CSV percentiles, and build/run documentation. +- Regression tests for setup validation, malformed files, dimension mismatches in checked + operations, non-default GEMM coefficients, stream ordering, move semantics, reusable + least-squares status buffers, and preallocated row extraction. + +### Changed + +- Column-major uploads copy directly from the supplied host vector; only row-major uploads + allocate a conversion buffer. +- `setStreamIdx` returns a reference for lvalues and moves rvalues into the returned tensor, + eliminating the previous implicit device allocation and deep copy. +- `getRows` uses one asynchronous 2D device copy instead of a cuBLAS call per row and + returns a tensor with the source's stream index. It shares its copy implementation with + `getRowsInto`; caller-supplied destinations retain their own stream index. +- Explicit tensor allocations require positive dimensions that fit CUDA's integer dimensions + and reject allocation-size overflow; default construction remains available for placeholders. + Vector constructors validate data length and storage mode before device allocation. +- Added validation for reshape sizes, slice axes/ranges, random bounds, storage modes, + row extraction, debug element access, empty factoriser inputs, and QR output batch shapes. +- Stream setup rejects zero streams and reconfiguration after session initialization; + stream selection validates against the actual initialized count. +- Documented unchecked numerical-operation preconditions, slice and asynchronous-buffer + lifetimes, stream dependencies, reusable workspace usage, and binary-format assumptions. +- Ignore generated CMake/CTest files, fetched dependencies, build products, and temporary + file-I/O test artifacts, including those produced by in-source builds. + +### Fixed + +- Enable testing at the project root so `ctest --test-dir build` discovers the test suite. +- Make the Python-format binary-reader test self-contained, removing its dependency on + separately generated Python test files and a particular working directory. +- Reject malformed text headers/data and incomplete or incorrectly sized binary files before + constructing tensors; reject extra payload data and report file-write failures. Handle short + file names safely and avoid writes into empty text-data vectors and uninitialized binary reads. +- Avoid null-pointer arithmetic when slicing a single-matrix tensor along the matrix axis. +- QR workspace queries, factorisation, and least-squares/Q extraction use the matrix-selected + cuSolver handle consistently with cuBLAS. The debugging `getQR` method synchronizes that + stream before reading the triangular factor on the host. +- Givens coefficient computation runs on the selected cuBLAS stream before the rotation, + checks kernel launch errors, and restores the caller's previous cuBLAS pointer mode. +- Configure test streams independently of test execution order and add nonblocking-stream + regressions with delayed inputs and downstream consumers for QR, Givens, and row extraction. +- Correct the row-slicing test to select a column before slicing its rows, rather than using + flat storage offsets as row indices. + + From 407cfde1d7aae81bf9b704b79f2b842c0601c6a7 Mon Sep 17 00:00:00 2001 From: alphaville Date: Tue, 3 Nov 2026 23:49:37 +0000 Subject: [PATCH 11/23] careful implementation of assignment op --- CHANGELOG.md | 8 +++ README.md | 31 ++++++++++- include/tensor.cuh | 78 +++++++++++++++++----------- test/testTensor.cu | 127 ++++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 211 insertions(+), 33 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 345a6c2..7f6e20d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +- **Breaking:** copy assignment now deep-copies tensor data and shape consistently with + copy construction, instead of creating a non-owning alias. Use explicit slices for + aliasing, and `deviceCopyTo` for copying into preallocated storage. Move assignment + transfers ownership and empties the source; moving a slice retains its borrowing semantics. - Column-major uploads copy directly from the supplied host vector; only row-major uploads allocate a conversion buffer. - `setStreamIdx` returns a reference for lvalues and moves rvalues into the returned tensor, @@ -49,6 +53,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Copy/move assignment releases the destination's previous owned storage, preserves + matrix-pointer and stream metadata, and handles self-assignment and empty tensors safely. + Copy construction and assignment complete copies on the source's selected stream. + Added regression coverage for ownership, allocation accounting, views, and stream ordering. - Enable testing at the project root so `ctest --test-dir build` discovers the test suite. - Make the Python-format binary-reader test self-contained, removing its dependency on separately generated Python test files and a particular working directory. diff --git a/README.md b/README.md index efe88af..c1ee358 100644 --- a/README.md +++ b/README.md @@ -94,8 +94,35 @@ DTensor y(10); x.deviceCopyTo(y); // x ---> 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 diff --git a/include/tensor.cuh b/include/tensor.cuh index 7b85fe6..a7b0abe 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -440,7 +440,20 @@ 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_idxStream, other.m_idxStream); + } /** * Slice constructor. @@ -674,8 +687,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); @@ -964,11 +981,21 @@ 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 @@ -1002,24 +1029,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 @@ -1293,12 +1304,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; } diff --git a/test/testTensor.cu b/test/testTensor.cu index 168bf66..89225a3 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -899,10 +899,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) { @@ -911,6 +918,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 * --------------------------------------- */ From 4e0922bef2aaf46dc0001184fae1f247c4908bb2 Mon Sep 17 00:00:00 2001 From: alphaville Date: Wed, 4 Nov 2026 00:06:33 +0000 Subject: [PATCH 12/23] fix Svd::rank() issue --- CHANGELOG.md | 4 +++ include/tensor.cuh | 23 +++++++++++++---- test/testTensor.cu | 64 ++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 86 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7f6e20d..542ce6e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -53,6 +53,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Repeated `Svd::rank()` calls reset device counters before counting, preventing + accumulation across calls and thresholds. SVD factorisation, rank reset/counting, + and result stream metadata consistently use the input matrix's selected stream. + Added repeated-call and nonblocking-stream regressions for float and double batches. - Copy/move assignment releases the destination's previous owned storage, preserves matrix-pointer and stream metadata, and handles self-assignment and empty tensors safely. Copy construction and assignment complete copies on the source's selected stream. diff --git a/include/tensor.cuh b/include/tensor.cuh index a7b0abe..bbbb0a8 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -1664,7 +1664,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()); } /** @@ -1707,16 +1711,25 @@ public: /** * 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; } @@ -1725,12 +1738,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)); } @@ -1747,7 +1760,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, @@ -1774,7 +1787,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, diff --git a/test/testTensor.cu b/test/testTensor.cu index 89225a3..ab3a9b2 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -1797,6 +1797,70 @@ 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 matrices(std::vector{5,0,0,2, 4,0,0,T(0.05)},2,2,2); + 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)); + 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(); +} + /* ================================================================================================ * CHOLESKY TESTS * ================================================================================================ */ From d41bd6f5a49c9b83c81845b7e41bcf089fc14d41 Mon Sep 17 00:00:00 2001 From: alphaville Date: Wed, 4 Nov 2026 00:21:54 +0000 Subject: [PATCH 13/23] fixed issues with tensor reshaping --- CHANGELOG.md | 15 +++++ README.md | 52 +++++++++++++++++ include/tensor.cuh | 119 ++++++++++++++++++++++++++------------ scripts/build_and_test.sh | 31 ++++++++++ scripts/check_memory.sh | 45 ++++++++++++++ test/testTensor.cu | 91 +++++++++++++++++++++++++++++ 6 files changed, 315 insertions(+), 38 deletions(-) create mode 100755 scripts/build_and_test.sh create mode 100755 scripts/check_memory.sh diff --git a/CHANGELOG.md b/CHANGELOG.md index 542ce6e..c7f063c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- `scripts/build_and_test.sh` builds the main executable and tests and runs all + C++/Python tests; `scripts/check_memory.sh` additionally runs full CUDA memcheck + with leak detection, shuffled repeated tests, main-executable checks, and saved logs. +- Size-safe manual accounting methods `Session::recordAllocation(size_t)` and + `recordDeallocation(size_t)`, plus a counted `cudaRelease` wrapper. - `DTensor::getRowsInto(destination, rowsFrom, rowsTo, matIdx)` copies rows into a preallocated tensor without allocating or resizing; validates indices, destination shape, and non-overlapping storage, and enqueues the copy on the source tensor's stream. @@ -28,6 +33,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +- **Breaking:** replace `incrementAllocatedBytes(int)` with separate size-safe + allocation/deallocation accounting methods. Reshape retains owned pointer-array + capacity when shrinking and reuses it when growing within capacity. - **Breaking:** copy assignment now deep-copies tensor data and shape consistently with copy construction, instead of creating a non-owning alias. Use explicit slices for aliasing, and `deviceCopyTo` for copying into preallocated storage. Move assignment @@ -53,6 +61,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Keep QR debugging-helper host uploads on the selected solver stream and complete + them before returning, fixing an intermittent nonblocking-stream failure exposed + by the shuffled full-suite sanitizer run. +- Track actual owned element bytes and matrix-pointer capacity independently of shape, + update counters only after successful allocation/free calls, and transfer sizes with + moves/swaps. Reshaped views own independent pointer arrays without freeing borrowed + data. Added reshape/reallocation/ownership and large-counter arithmetic regressions. - Repeated `Svd::rank()` calls reset device counters before counting, preventing accumulation across calls and thresholds. SVD factorisation, rank reset/counting, and result stream metadata consistently use the input matrix's selected stream. diff --git a/README.md b/README.md index c1ee358..7ad8b91 100644 --- a/README.md +++ b/README.md @@ -466,6 +466,19 @@ 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++ @@ -482,4 +495,43 @@ auto b = DTensor::createRandomTensor(3, 6, 4, -1, 1).setStreamIdx(1); Session::getInstance().synchronizeAllStreams(); ``` +## 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/include/tensor.cuh b/include/tensor.cuh index bbbb0a8..6a18a75 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -128,6 +128,25 @@ inline void gpuAssert(T code, const char *file, int line, bool abort = true) { */ 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. * Cuda library functions require a handle. @@ -185,7 +204,7 @@ private: 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: @@ -219,9 +238,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; } @@ -229,13 +249,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 @@ -308,19 +337,21 @@ 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) 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; } /** @@ -452,6 +483,8 @@ public: 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); } @@ -659,9 +692,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 * @@ -928,20 +962,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; @@ -1199,15 +1231,20 @@ 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)); + m_dataBytes = buffer_size; if (zero) gpuErrChk(cudaMemset(m_d_data, 0, buffer_size)); // set to zero all elements 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; } @@ -2070,10 +2107,12 @@ inline void QRFactoriser::getQR(DTensor &Q, DTensor &R) 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(m_matrix->streamIdx()), CUBLAS_SIDE_LEFT, CUBLAS_OP_N, m, n, n, @@ -2088,10 +2127,11 @@ inline void QRFactoriser::getQR(DTensor &Q, DTensor &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<> @@ -2106,10 +2146,12 @@ inline void QRFactoriser::getQR(DTensor &Q, DTensor &R) { 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(m_matrix->streamIdx()), CUBLAS_SIDE_LEFT, CUBLAS_OP_N, m, n, n, @@ -2124,10 +2166,11 @@ inline void QRFactoriser::getQR(DTensor &Q, DTensor &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. } diff --git a/scripts/build_and_test.sh b/scripts/build_and_test.sh new file mode 100755 index 0000000..e2f250d --- /dev/null +++ b/scripts/build_and_test.sh @@ -0,0 +1,31 @@ +#!/usr/bin/env bash +# Configure/build GPUtils and run every C++ and Python unit test. +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 +jobs="${JOBS:-4}" + +if ! "$python" -c 'import numpy' >/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}" \ + "$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/test/testTensor.cu b/test/testTensor.cu index ab3a9b2..98a5a07 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -1408,6 +1408,97 @@ 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 * --------------------------------------- */ From 1480b2fe8c645973a6c2b958ea4f4d6b644d37ea Mon Sep 17 00:00:00 2001 From: alphaville Date: Wed, 4 Nov 2026 00:35:25 +0000 Subject: [PATCH 14/23] numerically stable rank --- CHANGELOG.md | 12 +++ README.md | 49 ++++++++++ include/tensor.cuh | 129 +++++++++++++++++++++------ test/testTensor.cu | 218 ++++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 380 insertions(+), 28 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c7f063c..0942eff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Explicit scale-aware `Svd::rankRelative(rtol, atol = 0)` counts singular values + above `max(atol, rtol * sigmaMax)` using each matrix's own scale on the GPU. + Validates finite, nonnegative tolerances, shares the existing reusable rank buffer, + and preserves the absolute `rank()` API and nullspace policy. Added rescaling, + per-batch, zero-matrix, boundary, repeated-call, and nonblocking-stream tests. - `scripts/build_and_test.sh` builds the main executable and tests and runs all C++/Python tests; `scripts/check_memory.sh` additionally runs full CUDA memcheck with leak detection, shuffled repeated tests, main-executable checks, and saved logs. @@ -61,6 +66,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Complete selected-stream use for single and batched Cholesky, including workspace + queries and solves; solver status tensors follow the matrix stream. Host transfers, + debug scalar reads, initialization, and pointer-array uploads complete on their selected + stream. Transpose results inherit the source stream. Device copies order prior destination + work and complete on the source stream before returning. Added nonblocking-stream + regressions for Cholesky, transfers, and delayed SVD inputs, and documented cross-stream + operand dependencies and blocking transfer semantics. - Keep QR debugging-helper host uploads on the selected solver stream and complete them before returning, fixing an intermittent nonblocking-stream failure exposed by the shuffled full-suite sanitizer run. diff --git a/README.md b/README.md index 7ad8b91..4516c62 100644 --- a/README.md +++ b/README.md @@ -434,6 +434,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. @@ -495,6 +524,26 @@ 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 diff --git a/include/tensor.cuh b/include/tensor.cuh index 6a18a75..422c080 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -341,6 +341,20 @@ private: 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) { gpuErrChk(Session::getInstance().cudaRelease(m_d_data, m_dataBytes)); @@ -543,6 +557,8 @@ 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; @@ -798,7 +814,9 @@ 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 @@ -1154,7 +1172,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; } @@ -1165,7 +1183,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; } @@ -1176,7 +1194,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; } @@ -1187,7 +1205,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; } @@ -1232,7 +1250,11 @@ inline void DTensor::allocateOnDevice(size_t size, bool zero) { size_t buffer_size = size * sizeof(T); gpuErrChk(Session::getInstance().cudaAllocate((void**)&m_d_data, buffer_size)); m_dataBytes = buffer_size; - if (zero) gpuErrChk(cudaMemset(m_d_data, 0, buffer_size)); // set to zero all elements + if (zero) { + auto stream = selectedStream(); + gpuErrChk(cudaMemsetAsync(m_d_data, 0, buffer_size, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); + } if (numMats() > 1) { m_doDestroyPtrMatrices = true; @@ -1257,12 +1279,16 @@ inline bool DTensor::upload(const std::vector &vec, StorageMode mode) { 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"); + if (size == 0) return true; + auto stream = selectedStream(); if (mode == StorageMode::rowMajor) { std::vector vecCm(thisSize); rm2cm(vec, vecCm); - gpuErrChk(cudaMemcpy(m_d_data, vecCm.data(), size * sizeof(T), cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpyAsync(m_d_data, vecCm.data(), size * sizeof(T), cudaMemcpyHostToDevice, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); // Keep the conversion buffer alive. } else { - gpuErrChk(cudaMemcpy(m_d_data, vec.data(), size * sizeof(T), cudaMemcpyHostToDevice)); + gpuErrChk(cudaMemcpyAsync(m_d_data, vec.data(), size * sizeof(T), cudaMemcpyHostToDevice, stream)); + gpuErrChk(cudaStreamSynchronize(stream)); } return true; } @@ -1270,10 +1296,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 @@ -1290,6 +1319,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++) { @@ -1306,6 +1336,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++) { @@ -1324,10 +1355,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<> @@ -1409,10 +1445,8 @@ template inline T DTensor::operator()(size_t i, size_t j, size_t k) const { if (i >= m_numRows || j >= m_numCols || k >= m_numMats) throw std::invalid_argument("tensor index out of range"); - T hostDst; 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<> @@ -1640,6 +1674,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. @@ -1691,6 +1736,7 @@ public: checkMatrix(mat); m_destroyMatrix = destroyMatrix; m_tensor = (destroyMatrix) ? &mat : new DTensor(mat); + m_info->setStreamIdx(mat.streamIdx()); m_computeU = computeU; size_t m = mat.numRows(); size_t n = mat.numCols(); @@ -1770,6 +1816,32 @@ public: } 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; + } }; @@ -1868,6 +1940,7 @@ public: 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); } @@ -1890,7 +1963,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)); } @@ -1898,7 +1971,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)); } @@ -1906,7 +1979,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, @@ -1917,7 +1990,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, @@ -1927,7 +2000,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, @@ -1938,7 +2011,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, @@ -1977,6 +2050,7 @@ public: 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()); computeWorkspaceSize(); m_workspace = std::make_unique >(m_workspaceSize); m_householder = std::make_unique >(m_matrix->numCols()); @@ -2295,6 +2369,7 @@ public: m_factorisationDone(factorised) { 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(); } @@ -2314,7 +2389,7 @@ public: template<> inline void CholeskyBatchFactoriser::factorise() { if (m_factorisationDone) return; - gpuErrChk(cusolverDnDpotrfBatched(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnDpotrfBatched(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, m_numRows, m_matrix->ptrMatrices(), @@ -2327,7 +2402,7 @@ inline void CholeskyBatchFactoriser::factorise() { template<> inline void CholeskyBatchFactoriser::factorise() { if (m_factorisationDone) return; - gpuErrChk(cusolverDnSpotrfBatched(Session::getInstance().cuSolverHandle(), + gpuErrChk(cusolverDnSpotrfBatched(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, m_numRows, m_matrix->ptrMatrices(), @@ -2344,7 +2419,7 @@ 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(), + gpuErrChk(cusolverDnDpotrsBatched(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, m_numRows, 1, ///< only supports rhs = 1 @@ -2363,7 +2438,7 @@ 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(), + gpuErrChk(cusolverDnSpotrsBatched(Session::getInstance().cuSolverHandle(m_matrix->streamIdx()), CUBLAS_FILL_MODE_LOWER, m_numRows, 1, ///< only supports rhs = 1 diff --git a/test/testTensor.cu b/test/testTensor.cu index 98a5a07..7ddc1bf 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -107,6 +107,124 @@ static void CUDART_CB delayTestInput(void *) { 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; @@ -1922,11 +2040,13 @@ TEST_F(SvdTest, rankRepeatedCallsRecomputeCounts) { template void rankCallsAreOrderedOnSelectedStream() { NonblockingTestStream execution; - DTensor matrices(std::vector{5,0,0,2, 4,0,0,T(0.05)},2,2,2); + 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()); @@ -1952,6 +2072,102 @@ TEST_F(SvdTest, rankCallsUseSelectedStream) { 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 * ================================================================================================ */ From b3549c3c889a16e6ddefe12f6086579dd913c279 Mon Sep 17 00:00:00 2001 From: alphaville Date: Wed, 4 Nov 2026 00:41:25 +0000 Subject: [PATCH 15/23] added cudaStreamDestroy --- CHANGELOG.md | 3 +++ include/tensor.cuh | 3 +++ 2 files changed, 6 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 0942eff..b2c9240 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -66,6 +66,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Session shutdown completes work on all owned CUDA streams, releases library + handles, and explicitly destroys each stream instead of retaining stream resources + until CUDA context teardown. - Complete selected-stream use for single and batched Cholesky, including workspace queries and solves; solver status tensors follow the matrix stream. Host transfers, debug scalar reads, initialization, and pointer-array uploads complete on their selected diff --git a/include/tensor.cuh b/include/tensor.cuh index 422c080..7594ab0 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -195,9 +195,12 @@ private: } ~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])); } } From 6d01ba3af36bca5a3dbf4f3be9b5ee0cb1a2ba38 Mon Sep 17 00:00:00 2001 From: alphaville Date: Wed, 4 Nov 2026 00:47:47 +0000 Subject: [PATCH 16/23] memory leak revealed --- test/CMakeLists.txt | 1 + test/hostAllocationFailure.cpp | 36 +++++++++++++++++++ test/hostAllocationFailure.h | 22 ++++++++++++ test/testTensor.cu | 66 ++++++++++++++++++++++++++++++++++ 4 files changed, 125 insertions(+) create mode 100644 test/hostAllocationFailure.cpp create mode 100644 test/hostAllocationFailure.h diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4abfa1c..e768273 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -13,6 +13,7 @@ add_executable(gputils_test) target_sources(gputils_test # add files PRIVATE testTensor.cu + hostAllocationFailure.cpp ) target_link_libraries(gputils_test PRIVATE 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/testTensor.cu b/test/testTensor.cu index 7ddc1bf..0d6a00e 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -5,6 +5,7 @@ #include #include #include +#include "hostAllocationFailure.h" #define PRECISION_LOW 1e-4 #define PRECISION_HIGH 1e-10 @@ -366,6 +367,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 * --------------------------------------- */ From 17180cc37e179fb032503095b4bc37897707e764 Mon Sep 17 00:00:00 2001 From: alphaville Date: Wed, 4 Nov 2026 00:52:27 +0000 Subject: [PATCH 17/23] fixed mem mgmt issue --- CHANGELOG.md | 5 +++++ README.md | 6 ++++++ include/tensor.cuh | 22 +++++++++++++++------- 3 files changed, 26 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b2c9240..b776c5f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -66,6 +66,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Allocating tensor constructors release acquired device storage and preserve the + original exception when host allocation fails during row-major conversion or + matrix-pointer setup. Regression tests inject actual test-only host `operator new` + failures after GPU allocation, verify accounting recovery, and exercise successful + construction afterward for float, double, and int. - Session shutdown completes work on all owned CUDA streams, releases library handles, and explicitly destroys each stream instead of retaining stream resources until CUDA context teardown. diff --git a/README.md b/README.md index 4516c62..57baa93 100644 --- a/README.md +++ b/README.md @@ -14,6 +14,12 @@ 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 diff --git a/include/tensor.cuh b/include/tensor.cuh index 7594ab0..fbcebf2 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -1009,9 +1009,13 @@ DTensor::DTensor(size_t m, size_t n, size_t k, bool zero) { m_numRows = m; m_numCols = n; m_numMats = k; - allocateOnDevice(size, zero); - /* Initialise m_d_ptrMatrices */ - initialisePointersToMatricesData(); + try { + allocateOnDevice(size, zero); + initialisePointersToMatricesData(); + } catch (...) { + destroy(); + throw; + } } template @@ -1022,10 +1026,14 @@ DTensor::DTensor(const std::vector &data, size_t m, size_t n, size_t k, St m_numRows = m; m_numCols = n; m_numMats = k; - allocateOnDevice(size); - upload(data, mode); - /* Initialise m_d_ptrMatrices */ - initialisePointersToMatricesData(); + try { + allocateOnDevice(size); + upload(data, mode); + initialisePointersToMatricesData(); + } catch (...) { + destroy(); + throw; + } } template From 550046d4cab8efde3a7ad49975d34c267bf9b4a7 Mon Sep 17 00:00:00 2001 From: alphaville Date: Wed, 4 Nov 2026 01:02:52 +0000 Subject: [PATCH 18/23] using unique_ptr in Svd --- CHANGELOG.md | 4 +++ README.md | 6 ++++ include/tensor.cuh | 21 +++++++------- test/testTensor.cu | 68 ++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 88 insertions(+), 11 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b776c5f..66f1ae8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -66,6 +66,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- SVD uses a `unique_ptr` owner for its optional input copy and a separate raw + observer for numerical operations, removing manual deletion and preventing leaks + if later workspace/output construction throws. Fault-injection regressions verify + cleanup, input preservation, and borrowed-input lifetime for float and double. - Allocating tensor constructors release acquired device storage and preserve the original exception when host allocation fails during row-major conversion or matrix-pointer setup. Regression tests inject actual test-only host `operator new` diff --git a/README.md b/README.md index 57baa93..9976271 100644 --- a/README.md +++ b/README.md @@ -421,6 +421,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: diff --git a/include/tensor.cuh b/include/tensor.cuh index fbcebf2..5f817d3 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -1708,14 +1708,14 @@ 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. @@ -1745,8 +1745,12 @@ 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(); @@ -1793,13 +1797,8 @@ 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; /** diff --git a/test/testTensor.cu b/test/testTensor.cu index 0d6a00e..26edeb6 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -1925,6 +1925,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 From bd6a71c130905343b555f97c4ef708318e1cac00 Mon Sep 17 00:00:00 2001 From: alphaville Date: Wed, 4 Nov 2026 01:11:12 +0000 Subject: [PATCH 19/23] s_numStreams became a Session static member --- CHANGELOG.md | 4 ++++ README.md | 3 +++ include/tensor.cuh | 8 ++------ test/CMakeLists.txt | 1 + test/sessionAcrossTranslationUnits.cu | 19 +++++++++++++++++++ test/sessionAcrossTranslationUnits.h | 13 +++++++++++++ test/testTensor.cu | 23 ++++++++++++++++++++++- 7 files changed, 64 insertions(+), 7 deletions(-) create mode 100644 test/sessionAcrossTranslationUnits.cu create mode 100644 test/sessionAcrossTranslationUnits.h diff --git a/CHANGELOG.md b/CHANGELOG.md index 66f1ae8..d337964 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -66,6 +66,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Fixed +- Store the configured stream count in a C++17 inline static `Session` member, + eliminating per-translation-unit copies. A two-CUDA-translation-unit regression + verifies cross-file configuration, singleton identity, second-stream execution, + allocation accounting, and rejection of reconfiguration after initialization. - SVD uses a `unique_ptr` owner for its optional input copy and a separate raw observer for numerical operations, removing manual deletion and preventing leaks if later workspace/output construction throws. Fault-injection regressions verify diff --git a/README.md b/README.md index 9976271..fc17797 100644 --- a/README.md +++ b/README.md @@ -528,6 +528,9 @@ 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); diff --git a/include/tensor.cuh b/include/tensor.cuh index 5f817d3..e8777a5 100644 --- a/include/tensor.cuh +++ b/include/tensor.cuh @@ -122,12 +122,6 @@ 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. @@ -156,6 +150,8 @@ public: * 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: /** diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e768273..a4e0cdf 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -14,6 +14,7 @@ target_sources(gputils_test # add files PRIVATE testTensor.cu hostAllocationFailure.cpp + sessionAcrossTranslationUnits.cu ) target_link_libraries(gputils_test PRIVATE 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/testTensor.cu b/test/testTensor.cu index 26edeb6..ba7c252 100644 --- a/test/testTensor.cu +++ b/test/testTensor.cu @@ -6,6 +6,7 @@ #include #include #include "hostAllocationFailure.h" +#include "sessionAcrossTranslationUnits.h" #define PRECISION_LOW 1e-4 #define PRECISION_HIGH 1e-10 @@ -35,7 +36,27 @@ protected: * have the same stream setup and do not depend on test execution order. * --------------------------------------- */ -static const bool testStreamsConfigured = [] { Session::setStreams(2); return true; }(); +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); From 126f22021969fa6705428210469bee2ab2983bd7 Mon Sep 17 00:00:00 2001 From: alphaville Date: Wed, 4 Nov 2026 01:52:06 +0000 Subject: [PATCH 20/23] bug fix: QR memory allocation --- CHANGELOG.md | 153 +++++++++++++-------------------------------- include/tensor.cuh | 30 ++++++--- test/testTensor.cu | 58 +++++++++++++++++ 3 files changed, 120 insertions(+), 121 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d337964..4e65f0b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,114 +12,51 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added -- Explicit scale-aware `Svd::rankRelative(rtol, atol = 0)` counts singular values - above `max(atol, rtol * sigmaMax)` using each matrix's own scale on the GPU. - Validates finite, nonnegative tolerances, shares the existing reusable rank buffer, - and preserves the absolute `rank()` API and nullspace policy. Added rescaling, - per-batch, zero-matrix, boundary, repeated-call, and nonblocking-stream tests. -- `scripts/build_and_test.sh` builds the main executable and tests and runs all - C++/Python tests; `scripts/check_memory.sh` additionally runs full CUDA memcheck - with leak detection, shuffled repeated tests, main-executable checks, and saved logs. -- Size-safe manual accounting methods `Session::recordAllocation(size_t)` and - `recordDeallocation(size_t)`, plus a counted `cudaRelease` wrapper. -- `DTensor::getRowsInto(destination, rowsFrom, rowsTo, matIdx)` copies rows into a - preallocated tensor without allocating or resizing; validates indices, destination shape, - and non-overlapping storage, and enqueues the copy on the source tensor's stream. -- `DTensor::leastSquaresBatched(b, infoArray)` accepts a reusable device status buffer, - avoiding GPUtils's per-call workspace allocation while retaining the original overload. +- `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. -- Standalone Release-mode GPU benchmarks for pointer-based versus strided batched GEMM - and row-loop, 2D-copy, kernel, and actual `getRows` extraction paths. Includes configurable - dimensions, correctness checks, preallocated and allocation-inclusive measurements, - host/CUDA-event timings, CSV percentiles, and build/run documentation. -- Regression tests for setup validation, malformed files, dimension mismatches in checked - operations, non-default GEMM coefficients, stream ordering, move semantics, reusable - least-squares status buffers, and preallocated row extraction. +- 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 accounting methods. Reshape retains owned pointer-array - capacity when shrinking and reuses it when growing within capacity. -- **Breaking:** copy assignment now deep-copies tensor data and shape consistently with - copy construction, instead of creating a non-owning alias. Use explicit slices for - aliasing, and `deviceCopyTo` for copying into preallocated storage. Move assignment - transfers ownership and empties the source; moving a slice retains its borrowing semantics. -- Column-major uploads copy directly from the supplied host vector; only row-major uploads - allocate a conversion buffer. -- `setStreamIdx` returns a reference for lvalues and moves rvalues into the returned tensor, - eliminating the previous implicit device allocation and deep copy. -- `getRows` uses one asynchronous 2D device copy instead of a cuBLAS call per row and - returns a tensor with the source's stream index. It shares its copy implementation with - `getRowsInto`; caller-supplied destinations retain their own stream index. -- Explicit tensor allocations require positive dimensions that fit CUDA's integer dimensions - and reject allocation-size overflow; default construction remains available for placeholders. - Vector constructors validate data length and storage mode before device allocation. -- Added validation for reshape sizes, slice axes/ranges, random bounds, storage modes, - row extraction, debug element access, empty factoriser inputs, and QR output batch shapes. -- Stream setup rejects zero streams and reconfiguration after session initialization; - stream selection validates against the actual initialized count. -- Documented unchecked numerical-operation preconditions, slice and asynchronous-buffer - lifetimes, stream dependencies, reusable workspace usage, and binary-format assumptions. -- Ignore generated CMake/CTest files, fetched dependencies, build products, and temporary - file-I/O test artifacts, including those produced by in-source builds. +- **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 -- Store the configured stream count in a C++17 inline static `Session` member, - eliminating per-translation-unit copies. A two-CUDA-translation-unit regression - verifies cross-file configuration, singleton identity, second-stream execution, - allocation accounting, and rejection of reconfiguration after initialization. -- SVD uses a `unique_ptr` owner for its optional input copy and a separate raw - observer for numerical operations, removing manual deletion and preventing leaks - if later workspace/output construction throws. Fault-injection regressions verify - cleanup, input preservation, and borrowed-input lifetime for float and double. -- Allocating tensor constructors release acquired device storage and preserve the - original exception when host allocation fails during row-major conversion or - matrix-pointer setup. Regression tests inject actual test-only host `operator new` - failures after GPU allocation, verify accounting recovery, and exercise successful - construction afterward for float, double, and int. -- Session shutdown completes work on all owned CUDA streams, releases library - handles, and explicitly destroys each stream instead of retaining stream resources - until CUDA context teardown. -- Complete selected-stream use for single and batched Cholesky, including workspace - queries and solves; solver status tensors follow the matrix stream. Host transfers, - debug scalar reads, initialization, and pointer-array uploads complete on their selected - stream. Transpose results inherit the source stream. Device copies order prior destination - work and complete on the source stream before returning. Added nonblocking-stream - regressions for Cholesky, transfers, and delayed SVD inputs, and documented cross-stream - operand dependencies and blocking transfer semantics. -- Keep QR debugging-helper host uploads on the selected solver stream and complete - them before returning, fixing an intermittent nonblocking-stream failure exposed - by the shuffled full-suite sanitizer run. -- Track actual owned element bytes and matrix-pointer capacity independently of shape, - update counters only after successful allocation/free calls, and transfer sizes with - moves/swaps. Reshaped views own independent pointer arrays without freeing borrowed - data. Added reshape/reallocation/ownership and large-counter arithmetic regressions. -- Repeated `Svd::rank()` calls reset device counters before counting, preventing - accumulation across calls and thresholds. SVD factorisation, rank reset/counting, - and result stream metadata consistently use the input matrix's selected stream. - Added repeated-call and nonblocking-stream regressions for float and double batches. -- Copy/move assignment releases the destination's previous owned storage, preserves - matrix-pointer and stream metadata, and handles self-assignment and empty tensors safely. - Copy construction and assignment complete copies on the source's selected stream. - Added regression coverage for ownership, allocation accounting, views, and stream ordering. +- 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 self-contained, removing its dependency on - separately generated Python test files and a particular working directory. -- Reject malformed text headers/data and incomplete or incorrectly sized binary files before - constructing tensors; reject extra payload data and report file-write failures. Handle short - file names safely and avoid writes into empty text-data vectors and uninitialized binary reads. +- 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 workspace queries, factorisation, and least-squares/Q extraction use the matrix-selected - cuSolver handle consistently with cuBLAS. The debugging `getQR` method synchronizes that - stream before reading the triangular factor on the host. -- Givens coefficient computation runs on the selected cuBLAS stream before the rotation, - checks kernel launch errors, and restores the caller's previous cuBLAS pointer mode. -- Configure test streams independently of test execution order and add nonblocking-stream - regressions with delayed inputs and downstream consumers for QR, Givens, and row extraction. -- Correct the row-slicing test to select a column before slicing its rows, rather than using - flat storage offsets as row indices. +- 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. -## v1.10.0 - Unreleased +## v2.0.0 - 04-11-2026 ### Added