如何有效地将向量重复到cuda中的矩阵？

Question

我想在 cuda 中重复一个向量来形成一个矩阵，避免过多的 memcopy。向量和矩阵都分配在 GPU 上。

例如：

我有一个向量：

a = [1 2 3 4]

将其展开为矩阵：

b = [1 2 3 4;
     1 2 3 4;
     .......
     1 2 3 4]

我尝试的是分配 b 的每个元素。但这涉及到大量 GPU 内存到 GPU 内存的复制。

我知道这在 matlab 中很容易（使用 repmat），但是如何在 cuda 中有效地做到这一点？我没有在cublas中找到任何例程。

Answer 1

EDIT 基于 cmets，我已将代码更新为可处理行优先或列优先底层存储的版本。

这样的事情应该相当快：

// for row_major, blocks*threads should be a multiple of vlen
// for column_major, blocks should be equal to vlen
template <typename T>
__global__ void expand_kernel(const T* vector, const unsigned vlen, T* matrix, const unsigned mdim, const unsigned col_major=0){
  if (col_major){
    int idx = threadIdx.x+blockIdx.x*mdim;
    T myval = vector[blockIdx.x];
    while (idx < ((blockIdx.x+1)*mdim)){
      matrix[idx] = myval;
      idx += blockDim.x;
      }
    }
  else{
    int idx = threadIdx.x + blockDim.x * blockIdx.x;
    T myval = vector[idx%vlen];
    while (idx < mdim*vlen){
      matrix[idx] = myval;
      idx += gridDim.x*blockDim.x;
      }
    }
}

这假设您的矩阵的维度为 mdim 行 x vlen 列（似乎是您在问题中概述的内容。）

您可以调整网格和块的尺寸，以找出最适合您的特定 GPU 的方法。对于 row-major 情况，从每个块 256 或 512 个线程开始，并将块数设置为等于或大于 GPU 中 SM 数量的 4 倍。选择网格和块尺寸的乘积等于矢量长度vlen 的整数倍。如果这很困难，那么选择一个任意但“大”的线程块大小（例如 250 或 500）应该不会导致效率损失太大。

对于列优先的情况，每个块选择 256 或 512 个线程，并选择等于向量长度 vlen 的块数。如果vlen > 65535，则需要编译它以获得计算能力 3.0 或更高版本。如果vlen 很小，可能小于32，这种方法的效率可能会大大降低。如果您将每个块的线程数增加到 GPU 的最大值（512 或 1024），将会发现一些缓解措施。可能还有其他“扩展”实现可能更适合列主要“窄”矩阵的情况。例如，对列优先代码的直接修改将允许每个向量元素两个块，或每个向量元素四个块，然后启动的块总数将是 2*vlen 或 4*vlen，例如。

这是一个完整的示例，以及运行带宽测试，以证明上述代码实现了 bandwidthTest 指示的吞吐量的约 90%：

$ cat t546.cu
#include <stdio.h>

#define W 512
#define H (512*1024)
// for row_major, blocks*threads should be a multiple of vlen
// for column_major, blocks should be equal to vlen
template <typename T>
__global__ void expand_kernel(const T* vector, const unsigned vlen, T* matrix, const unsigned mdim, const unsigned col_major=0){
  if (col_major){
    int idx = threadIdx.x+blockIdx.x*mdim;
    T myval = vector[blockIdx.x];
    while (idx < ((blockIdx.x+1)*mdim)){
      matrix[idx] = myval;
      idx += blockDim.x;
      }
    }
  else{
    int idx = threadIdx.x + blockDim.x * blockIdx.x;
    T myval = vector[idx%vlen];
    while (idx < mdim*vlen){
      matrix[idx] = myval;
      idx += gridDim.x*blockDim.x;
      }
    }
}

template <typename T>
__global__ void check_kernel(const T* vector, const unsigned vlen, T* matrix, const unsigned mdim, const unsigned col_major=0){
  unsigned i = 0;
  while (i<(vlen*mdim)){
    unsigned idx = (col_major)?(i/mdim):(i%vlen);
    if (matrix[i] != vector[idx]) {printf("mismatch at offset %d\n",i); return;}
    i++;}
}

int main(){

  int *v, *m;
  cudaMalloc(&v, W*sizeof(int));
  cudaMalloc(&m, W*H*sizeof(int));
  int *h_v = (int *)malloc(W*sizeof(int));
  for (int i = 0; i < W; i++)
    h_v[i] = i;
  cudaMemcpy(v, h_v, W*sizeof(int), cudaMemcpyHostToDevice);

  // test row-major

  cudaEvent_t start, stop;
  cudaEventCreate(&start);
  cudaEventCreate(&stop);
  cudaEventRecord(start);
  expand_kernel<<<44, W>>>(v, W, m, H);
  cudaEventRecord(stop);
  float et;
  cudaEventSynchronize(stop);
  cudaEventElapsedTime(&et, start, stop);
  printf("row-majortime: %fms, bandwidth: %.0fMB/s\n", et, W*H*sizeof(int)/(1024*et));
  check_kernel<<<1,1>>>(v, W, m, H);
  cudaDeviceSynchronize();
  // test col-major

  cudaEventRecord(start);
  expand_kernel<<<W, 256>>>(v, W, m, H, 1);
  cudaEventRecord(stop);
  cudaEventSynchronize(stop);
  cudaEventElapsedTime(&et, start, stop);
  printf("col-majortime: %fms, bandwidth: %.0fMB/s\n", et, W*H*sizeof(int)/(1024*et));
  check_kernel<<<1,1>>>(v, W, m, H, 1);
  cudaDeviceSynchronize();
  return 0;
}

$ nvcc -arch=sm_20 -o t546 t546.cu
$ ./t546
row-majortime: 13.066944ms, bandwidth: 80246MB/s
col-majortime: 12.806720ms, bandwidth: 81877MB/s
$ /usr/local/cuda/samples/bin/x86_64/linux/release/bandwidthTest
[CUDA Bandwidth Test] - Starting...
Running on...

 Device 0: Quadro 5000
 Quick Mode

 Host to Device Bandwidth, 1 Device(s)
 PINNED Memory Transfers
   Transfer Size (Bytes)        Bandwidth(MB/s)
   33554432                     5864.2

 Device to Host Bandwidth, 1 Device(s)
 PINNED Memory Transfers
   Transfer Size (Bytes)        Bandwidth(MB/s)
   33554432                     6333.1

 Device to Device Bandwidth, 1 Device(s)
 PINNED Memory Transfers
   Transfer Size (Bytes)        Bandwidth(MB/s)
   33554432                     88178.6

Result = PASS
$

CUDA 6.5、RHEL 5.5

这也可以使用CUBLAS Rank-1 update function 来实现，但它会比上述方法慢很多。