c++ - 矩阵 vector 乘法 (cublasDgemv) 返回零

Question

对于我第一次接触 CUDA/cuBLAS，我尝试编写一个简单的函数，将 MxN 矩阵(用 vector 的 vector 表示，std::vector)与 Nx1 相乘“ones” vector ，以便获得矩阵的行(？)总和。这将利用 cublas_gemv() 加上其他基本 CUDA 操作，我认为这是一个很好的起点。

在处理设置问题和阅读/复制示例代码后，这是我所拥有的:

std::vector<double> test(std::vector<std::vector<double>> in)
{
    std::vector<double> out;
    long in_m = in.size();
    long in_n = in[0].size();
    cudaError_t cudaStat;
    cublasStatus_t stat;
    cublasHandle_t handle;
    // This just converts a vector-of-vectors into a col-first array
    double* p_in = vec2d_to_colfirst_array(in);
    double* p_ones = new double[in_n];
    double* p_out = new double[in_m];
    std::fill(p_ones, p_ones + in_n, 1.0);
    double* dev_in;
    double* dev_ones;
    double* dev_out;
    cudaStat = cudaMalloc((void**)&dev_in, in_m * in_n * sizeof(double));
    cudaStat = cudaMalloc((void**)&dev_ones, in_n * sizeof(double));
    cudaStat = cudaMalloc((void**)&dev_out, in_m * sizeof(double));
    stat = cublasCreate(&handle);
    cudaStat = cudaMemcpy(dev_in, p_in, in_m*in_n * sizeof(double), cudaMemcpyHostToDevice);
    cudaStat = cudaMemcpy(dev_ones, p_ones, in_n * sizeof(double), cudaMemcpyHostToDevice);
    double alpha = 1.0;
    double beta = 0.0;
    stat = cublasDgemv(handle, CUBLAS_OP_N, in_m, in_n, &alpha, dev_in, in_m, dev_ones, 1, &beta, dev_ones, 1);
    cudaStat = cudaMemcpy(p_out, dev_out, in_m * sizeof(double), cudaMemcpyDeviceToHost);
    out.assign(p_out, p_out + in_m);
    cudaFree(dev_in);
    cudaFree(dev_ones);
    cudaFree(dev_out);
    cublasDestroy(handle);
    free(p_in);
    free(p_ones);
    free(p_out);
    return out;
}

它看起来与我阅读的样本没有太大区别，所以我希望它“能正常工作”。但是，当我检查 p_out 时，它全为零。我肯定没有在 矩阵中输入零。

我验证了 vec2d_to_colfirst_array() 完成了它的工作，而且 dev_in/dev_ones 通过从设备返回主机，然后读取。也许问题出在对 cublasDgemv() 的调用中，但由于我是新手(而且因为与 Eigen 等相比，BLAS 语法更不直观)，在经历了很多挫折之后我就是不能看看有什么问题。

感谢任何帮助!

Answer 1

错误看起来相当简单。您希望从 dev_out 复制结果:

cudaStat = cudaMemcpy(p_out, dev_out, in_m * sizeof(double), cudaMemcpyDeviceToHost);

但是你永远不会在你的 cublas 调用中使用 dev_out:

stat = cublasDgemv(handle, CUBLAS_OP_N, in_m, in_n, &alpha, dev_in, in_m, dev_ones, 1, &beta, dev_ones, 1);

这似乎只是一个复制粘贴错误。如果您将 cublas 调用中的最后一个 dev_ones 实例替换为 dev_out，您的代码对我有用:

stat = cublasDgemv(handle, CUBLAS_OP_N, in_m, in_n, &alpha, dev_in, in_m, dev_ones, 1, &beta, dev_out, 1);

这是一个完整的示例，其中包含该更改:

$ cat t315.cu
#include <vector>
#include <cublas_v2.h>
#include <iostream>

const long idim1 = 8;
const long idim2 = 8;

double* vec2d_to_colfirst_array(std::vector<std::vector<double>> in){
    long dim1 = in.size();
    long dim2 = in[0].size();
    long k = 0;
    double *res = new double[dim1*dim2];
    for (int i = 0; i < dim1; i++)
      for (int j = 0; j < dim2; j++) res[k++] = in[i][j];
    return res;
}


std::vector<double> test(std::vector<std::vector<double>> in)
{
    std::vector<double> out;
    long in_m = in.size();
    long in_n = in[0].size();
    cudaError_t cudaStat;
    cublasStatus_t stat;
    cublasHandle_t handle;
    // This just converts a vector-of-vectors into a col-first array
    double* p_in = vec2d_to_colfirst_array(in);
    double* p_ones = new double[in_n];
    double* p_out = new double[in_m];
    std::fill(p_ones, p_ones + in_n, 1.0);
    double* dev_in;
    double* dev_ones;
    double* dev_out;
    cudaStat = cudaMalloc((void**)&dev_in, in_m * in_n * sizeof(double));
    cudaStat = cudaMalloc((void**)&dev_ones, in_n * sizeof(double));
    cudaStat = cudaMalloc((void**)&dev_out, in_m * sizeof(double));
    stat = cublasCreate(&handle);
    cudaStat = cudaMemcpy(dev_in, p_in, in_m*in_n * sizeof(double), cudaMemcpyHostToDevice);
    cudaStat = cudaMemcpy(dev_ones, p_ones, in_n * sizeof(double), cudaMemcpyHostToDevice);
    double alpha = 1.0;
    double beta = 0.0;
    stat = cublasDgemv(handle, CUBLAS_OP_N, in_m, in_n, &alpha, dev_in, in_m, dev_ones, 1, &beta, dev_out, 1);
    cudaStat = cudaMemcpy(p_out, dev_out, in_m * sizeof(double), cudaMemcpyDeviceToHost);
    out.assign(p_out, p_out + in_m);
    cudaFree(dev_in);
    cudaFree(dev_ones);
    cudaFree(dev_out);
    cublasDestroy(handle);

    free(p_in);
    free(p_ones);
    free(p_out);
    return out;
}

int main(){

  std::vector<double> a(idim2, 1.0);
  std::vector<std::vector<double>> b;
  for (int i = 0; i <  idim1; i++) b.push_back(a);
  std::vector<double> c = test(b);
  for (int i = 0; i < c.size(); i++) std::cout << c[i] << ",";
  std::cout << std::endl;
}

$ nvcc -std=c++11 -o t315 t315.cu -lcublas
t315.cu(24): warning: variable "cudaStat" was set but never used

t315.cu(25): warning: variable "stat" was set but never used

$ cuda-memcheck ./t315
========= CUDA-MEMCHECK
8,8,8,8,8,8,8,8,
========= ERROR SUMMARY: 0 errors
$

请注意，我不认为 free() 是与 new 一起使用的正确 API，但这似乎不是您问题的症结所在.