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