I've been learning CUDA on my own, and one project I was trying to do is implementing a linear algebra library. I created a mathematical vector class in C++ that computes the dot product via a GPU. I'm running into issues when it comes to freeing memory on the device via my class destructor.
Specifically I'm running into this custom error Vec destructor failed to delete device pointer: an illegal memory access was encountered, which happens when the Vec destructor is called.
// vec.hpp
#pragma once
#include <vector>
template <typename T>
class Vec {
private:
size_t size;
// points to a memory location on device NOT host
T* data_ptr;
public:
static_assert(std::is_same<T, float>::value || std::is_same<T, double>::value,
"Vec<T> can only be instantiated with float or double.");
// initialize vector of zeros of dimension size
Vec(size_t size);
// initialize from vector of values
Vec(const std::vector<T>& vec);
// destructor for device memory
~Vec();
// dot product [TODO: Add an exterior product]
template <typename U>
friend U operator%(const Vec<U>& A, const Vec<U>& B);
};
extern template class Vec<float>;
extern template class Vec<double>;
I've implemented these methods in this file:
#include "vec.hpp"
#include <stdexcept>
#include <iostream>
constexpr int BLOCK_SIZE = 512;
// dot product kernel using shared memory
template <typename T>
__global__ void vec_dot(const T* vecA, const T* vecB, int size, T* out) {
static_assert(std::is_arithmetic<T>::value, "T must be a numeric type");
__shared__ T sdata[BLOCK_SIZE];
int thread_index = threadIdx.x;
sdata[threadIdx.x] = 0;
size_t idx = threadIdx.x + blockDim.x*blockIdx.x;
while (idx < size) {
sdata[thread_index] += vecA[idx] * vecB[idx];
idx += blockDim.x * gridDim.x;
}
for (unsigned int s=blockDim.x/2; s>0; s>>=1) {
__syncthreads();
if (thread_index < s) sdata[thread_index] += vecA[thread_index + s] * vecB[thread_index + s];
}
if (thread_index == 0) atomicAdd(out, sdata[0]);
}
template <typename T>
T operator%(const Vec<T>& A, const Vec<T>& B) {
if (A.size != B.size) {
throw std::invalid_argument("Vectors are not the same size.");
}
T* out_device;
cudaMalloc(&out_device, sizeof(T));
vec_dot<<<100,1024>>>(A.data_ptr, B.data_ptr, A.size, out_device);
cudaDeviceSynchronize();
T out_host;
cudaMemcpy(&out_host, out_device, sizeof(T), cudaMemcpyDeviceToHost);
cudaFree(out_device);
return out_host;
}
template <typename T>
Vec<T>::~Vec() {
cudaError_t err = cudaFree(data_ptr);
if (err != cudaSuccess) {
std::cerr << "Vec destructor failed to delete device pointer: " << cudaGetErrorString(err) << std::endl;
}
}
template <typename T>
Vec<T>::Vec(size_t size) : size(size), data_ptr(nullptr) {
cudaError_t err = cudaMalloc(&data_ptr, size * sizeof(T));
if (err != cudaSuccess) {
std::cerr << "Vec malloc failed: " << cudaGetErrorString(err) << std::endl;
std::terminate();
}
cudaError_t memset_err = cudaMemset(data_ptr, 0, size * sizeof(T));
if (memset_err != cudaSuccess) {
std::cerr << "Setting memory to zeros failed: " << cudaGetErrorString(err) << std::endl;
std::terminate();
}
}
template <typename T>
Vec<T>::Vec(const std::vector<T>& vec) : size(vec.size()), data_ptr(nullptr) {
cudaError_t err = cudaMalloc(&data_ptr, size * sizeof(T));
if (err != cudaSuccess) {
std::cerr << "Vec malloc failed: " << cudaGetErrorString(err) << std::endl;
std::terminate();
}
cudaMemcpy(data_ptr, vec.data(), size * sizeof(T), cudaMemcpyHostToDevice);
}
template class Vec<float>;
template class Vec<double>;
template float operator%<float>(const Vec<float>&, const Vec<float>&);
template double operator%<double>(const Vec<double>&, const Vec<double>&);
My main.cpp is:
#include <iostream>
#include "vec.hpp"
#include <vector>
int main()
{
std::vector<double> vec;
for (int i=0; i<10000; i++) {
vec.push_back(1.0f);
}
Vec<double> vector(vec);
double product = vector%vector;
std::cout << product << std::endl;
return 0;
}
I tested this a little bit by checking that the destructor was not called twice (just added a print statement to verify it was called only once), but I'm not really sure why my code doesn't work since I allocate device memory on vector creation and the only time I wipe that memory is on destruction.
The problem is that you've set your BLOCK_SIZE at 512, but call your block with 1024 threads. You should use the same BLOCK_SIZE everywhere to make sure you do not overflow the shared buffer.
Here's the corrected code (also on Godbolt: https://cuda.godbolt.org/z/or38M3hK3):
#include <vector>
#include <iostream>
#include <vector>
#include <stdexcept>
#include <iostream>
#include <cuda.h>
#include <cassert>
void check(const auto err, const int line) {
if (err != cudaSuccess) {
std::cerr << "error: " << cudaGetErrorString(err) << "in line: " << line << "\n";
}
}
#define CUDA(err) check(err, __LINE__)
template <typename T>
class Vec {
size_t size;
// points to a memory location on device NOT host
T* data_ptr;
public:
static_assert(std::is_same<T, float>::value || std::is_same<T, double>::value,
"Vec<T> can only be instantiated with float or double.");
// initialize vector of zeros of dimension size
Vec(size_t size);
// initialize from vector of values
Vec(const std::vector<T>& vec);
// destructor for device memory
~Vec();
// dot product [TODO: Add an exterior product]
template <typename U>
friend U dot_product(const Vec<U>& A, const Vec<U>& B);
};
extern template class Vec<float>;
extern template class Vec<double>;
constexpr auto BLOCK_SIZE = 512u;
static_assert(BLOCK_SIZE <= 1024);
// dot product kernel using shared memory
template <typename T>
__global__ void vec_dot(const T* vecA, const T* vecB, const int size, T* out) {
static_assert(std::is_arithmetic<T>::value, "T must be a numeric type");
assert(BLOCK_SIZE == blockDim.x);
__shared__ T sdata[BLOCK_SIZE];
const auto tid = threadIdx.x;
sdata[tid] = 0;
size_t idx = threadIdx.x + blockDim.x * blockIdx.x;
while (idx < size) {
sdata[tid] += vecA[idx] * vecB[idx];
idx += blockDim.x * gridDim.x;
}
for (auto s = blockDim.x/2; s > 0; s >>= 1) {
__syncthreads();
//every block in {}, even if it's a single statement
if (tid < s) { sdata[tid] += vecA[tid + s] * vecB[tid + s]; }
}
__syncthreads();
if (tid == 0) { atomicAdd(out, sdata[0]); }
}
template <typename T>
T dot_product(const Vec<T>& A, const Vec<T>& B) {
if (A.size != B.size) {
throw std::invalid_argument("Vectors are not the same size.");
}
T* out_device;
CUDA(cudaMalloc(&out_device, sizeof(T)));
vec_dot<<<100, BLOCK_SIZE>>>(A.data_ptr, B.data_ptr, A.size, out_device);
CUDA(cudaDeviceSynchronize());
T out_host;
CUDA(cudaMemcpy(&out_host, out_device, sizeof(T), cudaMemcpyDeviceToHost));
CUDA(cudaFree(out_device));
return out_host;
}
template <typename T>
Vec<T>::~Vec() {
CUDA(cudaFree(data_ptr));
}
template <typename T>
Vec<T>::Vec(size_t size) : size(size), data_ptr(nullptr) {
CUDA(cudaMalloc(&data_ptr, size * sizeof(T)));
CUDA(cudaMemset(data_ptr, 0, size * sizeof(T)));
}
template <typename T>
Vec<T>::Vec(const std::vector<T>& vec) : size(vec.size()) {
CUDA(cudaMalloc(&data_ptr /*nullptr*/, size * sizeof(T)));
CUDA(cudaMemcpy(data_ptr, vec.data(), size * sizeof(T), cudaMemcpyHostToDevice));
}
template class Vec<float>;
template class Vec<double>;
template float dot_product<float>(const Vec<float>&, const Vec<float>&);
template double dot_product<double>(const Vec<double>&, const Vec<double>&);
int main()
{
std::vector<float> vec;
for (auto i = 0; i < 10000; i++) {
vec.push_back(1.0f);
}
Vec<float> vector(vec);
auto product = dot_product(vector, vector);
std::cout << "dot product = " << product << "\n";
return 0;
}
Note the error checking macro, you should bracket every call to a cudaXXX function with error checking. That way you'll get the error close to the source.
See: What is the canonical way to check for errors using the CUDA runtime API?
This is important, because error reporting in CUDA is sticky. If an error occurs somewhere, it will keep reporting that error in every subsequent cuda call.
Also as a matter of C++ style, it's best not to use non-standard operator overloads. If no canonical operator is available, then just use a named function. That way it's clear that you're doing a dot product and not a modulus operation.
This version outputs:
Program returned: 0
Program stdout
dot product = 920
I don't know if that's the correct result (and I'm too lazy to check). It's best to add an assert to check that your code does indeed calculate correctly.