CUDA：在内核中使用设备函子

Question

我试图制作一个设备仿函数，它基本上像这样执行（未优化的）矩阵向量乘法

namespace cusolve
{

template <class value_type,
          class matrix_type = value_type*,
          class vector_type = value_type*>
struct linear_operator 
{
    const matrix_type matrix;
    const size_t width;
    __device__
    linear_operator(const matrix_type matrix, size_t width)
        : matrix(matrix), width(width) { }

    __device__
    void operator()(const vector_type x, vector_type x_out)
    {
        auto col = blockIdx.x * blockDim.x + threadIdx.x;
        auto row = blockIdx.y * blockDim.y + threadIdx.y;
        x_out[row] = 0;
        if (row < width)
        {
            for (size_t i = 0; i < width; i++)
            {
                x_out[row] += matrix[row*width + i] * x[i];
            }
        }
        return;              
    }
};

因此，这假定 matrix、x 和 x_out 是设备指针。所以，为了测试它，我尝试从一个简单的内核调用它

__global__
void
operateKernel(double *d_matrix,
        double *d_vector, double *d_vector_out,
        size_t width)
{
    cusolve::linear_operator<double> matmul(d_matrix, width);
    matmul(d_vector, d_vector_out);
}


void
operate(double *matrix, double *vector, double *vector_out, size_t width)
{
    const dim3 blockConfig(16, 16);
    const size_t gridWidth = (size_t) ((double) width) / 16.0l;
    const dim3 gridConfig(gridWidth, gridWidth);

    double *d_matrix, *d_vector, *d_vector_out;
    auto mem_vector = width * sizeof(double);
    auto mem_matrix = mem_vector * width;
    cudaMalloc((void **) &d_matrix, mem_matrix);
    cudaMalloc((void **) &d_vector, mem_vector);
    cudaMalloc((void **) &d_vector_out, mem_vector);

    cudaMemcpy(d_matrix, matrix, mem_matrix, cudaMemcpyHostToDevice);
    cudaMemcpy(d_vector, vector, mem_vector, cudaMemcpyHostToDevice);

    operateKernel<<<gridConfig, blockConfig>>>(d_matrix, d_vector, d_vector_out, width);
    cudaMemcpy(vector_out, d_vector_out, mem_vector, cudaMemcpyDeviceToHost);
    cudaFree(d_vector);
    cudaFree(d_matrix);
    cudaFree(d_vector_out);
}

但是，当我尝试使用分配并初始化为非空向量和矩阵从 main() 调用 operate() 时，输出全为零。我已经为此苦恼了一段时间，但一直无法弄清楚我做错了什么。 P.S：我故意尝试这样做，而不是作为学习练习。

Answer 1

在计算网格尺寸时忘记使用ceil。

const size_t gridWidth = ceil( ((double) width) / 16.0l );