用 cuda C++ 对向量值求和

Question

我尝试使用 CUDA c++ 对许多向量值求和。我找到了两个向量的一些解决方案。如您所见，可以添加两个向量，但我想动态生成具有相同长度的向量。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// CUDA kernel. Each thread takes care of one element of c
 __global__ void vecAdd(double *a, double *b, double *c, int n)
{
// Get our global thread ID
int id = blockIdx.x*blockDim.x+threadIdx.x;
// Make sure we do not go out of bounds
if (id < n)
    c[id] = a[id] + b[id];
}

int main( int argc, char* argv[] )
{
// Size of vectors
int n = 100000;

// Host input vectors
double *h_a;
double *h_b;
//Host output vector
double *h_c;

// Device input vectors
double *d_a;
double *d_b;
//Device output vector
double *d_c;

// Size, in bytes, of each vector
size_t bytes = n*sizeof(double);

// Allocate memory for each vector on host
h_a = (double*)malloc(bytes);
h_b = (double*)malloc(bytes);
h_c = (double*)malloc(bytes);

// Allocate memory for each vector on GPU
cudaMalloc(&d_a, bytes);
cudaMalloc(&d_b, bytes);
cudaMalloc(&d_c, bytes);

int i;
// Initialize vectors on host
for( i = 0; i < n; i++ ) {
    h_a[i] = sin(i)*sin(i);
    h_b[i] = cos(i)*cos(i);
}

// Copy host vectors to device
cudaMemcpy( d_a, h_a, bytes, cudaMemcpyHostToDevice);
cudaMemcpy( d_b, h_b, bytes, cudaMemcpyHostToDevice);

int blockSize, gridSize;

// Number of threads in each thread block
blockSize = 1024;

// Number of thread blocks in grid
gridSize = (int)ceil((float)n/blockSize);

// Execute the kernel
vecAdd<<<gridSize, blockSize>>>(d_a, d_b, d_c, n);

// Copy array back to host
cudaMemcpy( h_c, d_c, bytes, cudaMemcpyDeviceToHost );

// Sum up vector c and the print result divided by n, this should equal 1 
within error
double sum = 0;
for(i=0; i<n; i++)
    sum += h_c[i];
printf("final result: %f\n", sum/n);

// Release device memory
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);

// Release host memory
free(h_a);
free(h_b);
free(h_c);

return 0;
}

有没有办法为许多向量做到这一点？我的向量大小是：

#vector length
N = 1000 
#number of vectors
i = 300000
v[i] = [1,2,..., N]

结果我需要得到：

out[i]= [sum(v[1]), sum(v[2]),..., sum(v[i])]

感谢您的建议。

Answer 1

以类似于您显示的代码的方式将多个向量相加（即生成元素总和）相当于对矩阵的列求和。而这个想法代表了实现解决方案的明智方法。

我们会将您的向量视为一个矩阵，其中每个向量都是矩阵中的一行。 CUDA 内核将为每一列分配一个线程，并将该列的元素相加，产生一个数字结果。该单个数字结果将成为整个问题向量结果的一个元素。

这是一个完整的示例，展示了一种可能的方法：

$ cat t2.cu
#include <iostream>

typedef double mt;
const int nTPB = 64;

template <typename T>
__global__ void column_sum(T *matrix, T *sums, unsigned n_vectors, unsigned vector_length){

  unsigned idx = threadIdx.x+blockDim.x*blockIdx.x;
  if (idx < vector_length){
    T temp = 0;
    for (unsigned i = 0; i < n_vectors; i++)
      temp += matrix[i*vector_length+idx];
    sums[idx] = temp;}
}

int main(){
  const unsigned vlen = 1000;
  const unsigned nvec = 300000;
  mt *h_matrix, *d_matrix, *h_sums, *d_sums;
  // create the desired number of vectors as a single matrix
  h_sums = new mt[vlen];
  h_matrix = new mt[vlen*nvec];
  cudaMalloc(&d_matrix, vlen*nvec*sizeof(mt));
  cudaMalloc(&d_sums, vlen*sizeof(mt));
  size_t count = 0;
  for (unsigned i = 0; i < nvec; i++)
    for (unsigned j = 0; j < vlen; j++)
      h_matrix[count++] = j;
  cudaMemcpy(d_matrix, h_matrix, vlen*nvec*sizeof(mt), cudaMemcpyHostToDevice);
  column_sum<<<(vlen+nTPB-1)/nTPB,nTPB>>>(d_matrix, d_sums, nvec, vlen);
  cudaMemcpy(h_sums, d_sums, vlen*sizeof(mt), cudaMemcpyDeviceToHost);
  for (unsigned i = 0; i < vlen; i++) if (h_sums[i] != ((mt)nvec)*i) {std::cout << " mismatch at " << i << " was: " << h_sums[i] << " should be: " << ((mt)nvec)*i << std::endl; return -1;}
  std::cout << cudaGetErrorString(cudaGetLastError()) << std::endl;
}
$ nvcc -o t2 t2.cu
$ cuda-memcheck ./t2
========= CUDA-MEMCHECK
no error
========= ERROR SUMMARY: 0 errors
$

请注意，此方法仅在 GPU 上创建与矢量元素一样多的线程（上例中为 1000 个）。 1000 个线程足以让最小的 GPU 保持忙碌。但是，如果您的向量长度为​​ 10,000 或更长，则此算法在大多数 GPU 上将是有效的。如果您想探索为小问题创建更有效的算法，您可以研究classical parallel reduction 的想法。