A CUDA compatible Vector class

Question

Please review the following CUDA-compatible Vector class in C++.

#include <cstdlib>
#include <ctime>
#include <iostream>

class Vector {
private:
    int length;
    float * elements;
    bool is_column; // column-vector = single column, multiple rows.
    // row-vector = single row, multiple columns.
public:
    Vector()
    {
        this->length = 0;
        this->elements = nullptr;
        this->is_column = true;
    }
    Vector(int length)
    {
        this->set_length(length);
        is_column = true;
    }
    void set_length(int len)
    {
        this->length = len;
        if(elements!= nullptr)
        {
            cudaFree(elements);
        }
        cudaMallocManaged(&elements, len * sizeof(float));
    }
    void set_row_vector(bool is_row)
    {
        is_column = !is_row;
    }
    int get_length() const
    {
        return length;
    }
    // Copy constructor
    Vector(const Vector& other)
    {
        this->length = other.length;
        cudaMallocManaged(&elements, length * sizeof(float));
        for (int i = 0; i < length; i++) {
            this->elements[i] = other.elements[i];
        }
        this->is_column = true;
    }
    // Destructor
    ~Vector()
    {
        if(this->elements!=nullptr) {
            this->length = 0;
            cudaFree(this->elements);
        }
    }

    float dot_product(const Vector& other) const // inner product
    {
        if (!(!this->is_column && other.is_column_vector())) {
            throw std::runtime_error("Error: lhs must be row vector, rhs must be column vector");
        }

        if (this->length != other.length) {
            throw std::runtime_error("Error: Vectors must have the same length to perform Dot product");
        }
        float result = 0.0;
        for (int i = 0; i < this->length; i++) {
            result += this->elements[i] * other.elements[i];
        }
        return result;
    }

    void Vector::hadamard_product(const Vector& rhs, Vector** returns)
    {
        if (!(this->is_column && !rhs.is_column_vector())) {
            throw std::runtime_error("Error: rhs must be a row-vector");
        }
        *returns = new Vector[this->length];
        int block_size = 256;
        int num_blocks = (this->length + block_size - 1) / block_size;
        float* result;
        cudaMallocManaged(&result, this->length * sizeof(float));
        hadamard_product_kernel<<<num_blocks, block_size>>>(this->elements, rhs.elements, result, this->length);
        cudaDeviceSynchronize();
        int index = 0;
        for (int i = 0; i < this->length; i++) {
            Vector temp = rhs.clone();
            temp.scalar_product(result[i]);
            (*returns)[index] = temp;
            index++;
        }
        cudaFree(result);
    }
    
    void scalar_product(float val)
    {
        for (int i = 0; i < length; i++) {
            elements[i] *= val;
        }
    }
    
    // Overloaded assignment operator
    Vector& operator=(const Vector& other)
    {
        if (this != &other)
        {
            this->length = other.length;
            if(this->elements != nullptr)
            {
                cudaFree(this->elements);
            }
            cudaMallocManaged(&elements, length * sizeof(float));
            for (int i = 0; i < length; i++) {
                this->elements[i] = other.elements[i];
            }
        }
        return *this;
    }
    float& operator[](int index) // will modify the state of the object
    // when assigned a value
    {
        return elements[index];
    }
    bool operator==(const Vector& other) const
    {
        if (this->length != other.length) {
            return false;
        }
        for (int i = 0; i < this->length; i++) {
            if (this->elements[i] != other.elements[i]) {
                return false;
            }
        }
        return true;
    }
    bool is_column_vector() const
    {
        return is_column;
    }
};

__global__ void hadamard_product_kernel(float* lhs, float* rhs, float* result, int length)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < length) {
        result[idx] = lhs[idx] * rhs[idx];
    }
}
```

Answer 1

Separate concerns

Your class has too many responsibilities. It is:

1. Being a container like std::vector.
2. Doing CUDA memory management.
3. Adding row/column vector product functionality.

It might look like reasonable code now, but what if you add support for more mathematical operations? What if you want to support OpenCL? What if you want it to work with doubles instead of floats? It's best to separate concerns to keep everything managable.

First, the standard library already has a variable-length vector class: std::vector. It would be great if we could reuse that. By default, it uses the default memory allocator under the hood to allocate memory. However, you want to use CUDA functions to allocate memory. Luckily, you can provide a custom memory allocator to most STL contains:

template<typename T>
class CudaAllocator {
public:
    using value_type = T;

    T* allocate(std::size_t n) {
        T* result;
        cudaMallocManaged(&result, n * sizeof(T));
        return result;
    }

    void deallocate(T* p, std::size_t) {
        cudaFree(p);
    }
};

Now you can create a vector with memory managed by CUDA like so:

std::vector<float, CudaAllocator<float>> cudavec;

You can make a type alias for it to save some typing:

template<typename T = float>
using CudaVector = std::vector<T, CudaAllocator<T>>;

Instead of making the math functions member functions of a class, you can write out-of-class functions taking two parameters. For example:

template<typename T>
T dot_product(const CudaVector<T>& lhs, const CudaVector<T>& rhs)
{
    …
    if (lhs.size() != rhs.size()) {
        throw std::runtime_error("Error: Vectors must have the same length to perform Dot product");
    }

    T result{};

    for (std::size_t i = 0; i < lhs.size(); ++i) {
        result += lhs[i] * rhs[i];
    }

    return result;
}

The remaining issue is that you made a distinction between row and column vectors. You can create new types (not just aliases) for these:

tempate<typename T = float>
class CudaRowVector: public CudaVector<T> {
public:
    using CudaVector<T>::CudaVector; // expose the constructor of CudaVector<>
};

tempate<typename T = float>
class CudaColumnVector: public CudaVector<T> {
public:
    using CudaVector<T>::CudaVector;
};

And then the declaration of the dot product operation becomes:

template<typename T>
T dot_product(const CudaRowVector<T>& lhs, const CudaColumnVector<T>& rhs)
{
    …
}

This has the benefit that row/column errors are now caught at compile-time instead of at run-time.

Make use of standard library algorithms

The standard library comes with lots of algorithms that you can use, if your class works like other STL containers, for example by providing begin() and end() member functions so it can be iterated over. Since the CudaVector I've shown is just a std::vector, it means the algorithms will work on it as well. So you can simplify dot_product() a lot by making use of std::inner_product():

template<typename T>
T dot_product(const CudaRowVector<T>& lhs, const CudaColumnVector<T>& rhs)
{
    if (lhs.size() != rhs.size()) {
        throw std::runtime_error("Error: Vectors must have the same length to perform Dot product");
    }

    return std::inner_product(lhs.begin(), lhs.end(), rhs.begin(), T{});
}

I am sure there are CUDA equivalents of library functions for vector and matrix manipulations, so you could use those instead of writing your own kernels.