CUDA：如何使用推力进行矩阵乘法？

Question

我是 CUDA 和 Thrust 的新手，我正在尝试实现矩阵乘法，我想只使用推力算法来实现这一点，因为我想避免手动调用内核。

有没有办法可以有效地实现这一目标？ （至少不使用 2 个嵌套的 for 循环）

还是我必须辞职并调用 CUDA 内核？

//My data
thrust::device_vector<float> data(n*m);
thrust::device_vector<float> other(m*r);
thrust::device_vector<float> result(n*r);

// To make indexing faster, not really needed
transpose(other);

// My current approach
for (int i = 0; i < n; ++i)
{
   for (int j = 0; j < r;++j)
   {
       result[i*r+ j] = thrust::inner_product(data.begin()+(i*m), data.begin()+((i+1)*m),other+(j*m), 0.0f);
   }
}

Answer 1

如果您对性能感兴趣（通常是人们使用 GPU 执行计算任务的原因），则不应使用推力，也不应调用或编写自己的 CUDA 内核。您应该使用 CUBLAS 库。学习练习，如果想学习自己的CUDA内核，可以参考a first-level-optimized CUDA version in the CUDA programming guide in the shared memory section。如果你真的想在一次推力调用中使用推力，这是可能的。

基本思想是使用元素方面的操作，如 here 所述的推力::变换。每个输出数组元素的点积是使用由循环组成的函子计算的。

这是一个考虑 3 种方法的工作示例。您原来的双嵌套循环方法（相对较慢）、单推力调用方法（更快）和 cublas 方法（最快，当然对于较大的矩阵大小）。下面的代码只对 200 或更小的矩阵边尺寸运行方法 1，因为它太慢了。以下是 Tesla P100 的示例：

$ cat t463.cu
#include <thrust/device_vector.h>
#include <thrust/transform.h>
#include <thrust/inner_product.h>
#include <thrust/execution_policy.h>
#include <thrust/equal.h>
#include <thrust/iterator/constant_iterator.h>
#include <cublas_v2.h>
#include <iostream>
#include <time.h>
#include <sys/time.h>
#include <cstdlib>
#define USECPSEC 1000000ULL

long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv, 0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}

struct dp
{
  float *A, *B;
  int m,n,r;
  dp(float *_A, float *_B, int _m, int _n, int _r): A(_A), B(_B), m(_m), n(_n), r(_r) {};
  __host__ __device__
  float operator()(size_t idx){
    float sum = 0.0f;
    int row = idx/r;
    int col = idx - (row*r); // cheaper modulo
    for (int i = 0; i < m; i++)
      sum += A[col + row*i] * B[col + row*i];
    return sum;}
};

const int dsd = 200;
int main(int argc, char *argv[]){
  int ds = dsd;
  if (argc > 1) ds = atoi(argv[1]);
  const int n = ds;
  const int m = ds;
  const int r = ds;
  // data setup
  thrust::device_vector<float> data(n*m,1);
  thrust::device_vector<float> other(m*r,1);
  thrust::device_vector<float> result(n*r,0);
  // method 1
  //let's pretend that other is (already) transposed for efficient memory access by thrust
  // therefore each dot-product is formed using a row of data and a row of other
  long long dt = dtime_usec(0);
    if (ds < 201){
    for (int i = 0; i < n; ++i)
    {
      for (int j = 0; j < r;++j)
      {
         result[i*r+ j] = thrust::inner_product(data.begin()+(i*m), data.begin()+((i+1)*m),other.begin()+(j*m), 0.0f);
      }
    }
    cudaDeviceSynchronize();
    dt = dtime_usec(dt);
    if (thrust::equal(result.begin(), result.end(), thrust::constant_iterator<float>(m)))
      std::cout << "method 1 time: " << dt/(float)USECPSEC << "s" << std::endl;
    else
      std::cout << "method 1 failure" << std::endl;
    }
  thrust::fill(result.begin(), result.end(), 0);
  cudaDeviceSynchronize();
// method 2
  //let's pretend that data is (already) transposed for efficient memory access by thrust
  // therefore each dot-product is formed using a column of data and a column of other
  dt = dtime_usec(0);
  thrust::transform(thrust::counting_iterator<int>(0), thrust::counting_iterator<int>(n*r), result.begin(), dp(thrust::raw_pointer_cast(data.data()), thrust::raw_pointer_cast(other.data()), m, n, r));
  cudaDeviceSynchronize();
  dt = dtime_usec(dt);
  if (thrust::equal(result.begin(), result.end(), thrust::constant_iterator<float>(m)))
    std::cout << "method 2 time: " << dt/(float)USECPSEC << "s" << std::endl;
  else
    std::cout << "method 2 failure" << std::endl;
// method 3
  // once again, let's pretend the data is ready to go for CUBLAS
  cublasHandle_t h;
  cublasCreate(&h);
  thrust::fill(result.begin(), result.end(), 0);
  float alpha = 1.0f;
  float beta = 0.0f;
  cudaDeviceSynchronize();
  dt = dtime_usec(0);
  cublasSgemm(h, CUBLAS_OP_T, CUBLAS_OP_T, n, r, m, &alpha, thrust::raw_pointer_cast(data.data()), n, thrust::raw_pointer_cast(other.data()), m, &beta, thrust::raw_pointer_cast(result.data()), n);
  cudaDeviceSynchronize();
  dt = dtime_usec(dt);
  if (thrust::equal(result.begin(), result.end(), thrust::constant_iterator<float>(m)))
    std::cout << "method 3 time: " << dt/(float)USECPSEC << "s" << std::endl;
  else
    std::cout << "method 3 failure" << std::endl;
}
$ nvcc -o t463 t463.cu -lcublas
$ ./t463
method 1 time: 20.1648s
method 2 time: 6.3e-05s
method 3 time: 5.7e-05s
$ ./t463 1024
method 2 time: 0.008063s
method 3 time: 0.000458s
$

对于默认维度 200 的情况，单推力调用和 cublas 方法相当接近，但比循环方法快得多。对于 1024 的边尺寸，cublas 方法比单推力调用方法快近 20 倍。

请注意，我为所有 3 种方法选择了“最佳”转置配置。对于方法 1，最佳情况下的时间是 inner_product 使用每个输入矩阵中的“行”（实际上是第二个输入矩阵的转置）。对于方法 2，最佳情况下的时间是当函子从每个输入矩阵（实际上是第一个输入矩阵的转置）中遍历一个“列”时。对于方法 3，为两个输入矩阵选择 CUBLAS_OP_T 似乎是最快的。实际上，只有 cublas 方法具有灵活性，可用于各种输入情况并具有良好的性能。