c++ 在GPU上使用popcnt

uqcuzwp8  于 2023-04-08  发布在  其他
关注(0)|答案(2)|浏览(303)

我需要计算
(a & b).count()
在大集合(〉10000)位向量(std::bitset<N>)上,其中N是从2^10到2^16的任何值。

const size_t N = 2048;
std::vector<std::vector<char>> distances;
std::vector<std::bitset<N>> bits(100000);
load_from_file(bits);
for(int i = 0; i < bits.size(); i++){
    for(int j = 0; j < bits.size(); j++){
        distance[i][j] = (bits[i] & bits[j]).count();
    }
}

目前我依靠分块多线程和SSE/AVX来计算distances。幸运的是,我可以使用AVX的vpand来计算&,但我的代码仍然使用popcnt (%rax)和一个循环来计算位数。
有没有一种方法可以在我的GPU(nVidia 760 m)上计算(a & b).count()函数?理想情况下,我只需要传递2块N位的内存。我正在寻找使用推力,但我找不到popcnt函数。

编辑:

当前CPU实现。

double validate_pooled(const size_t K) const{                           
    int right = 0;                                                          
    const size_t num_examples = labels.size();                              
    threadpool tp;                                                          
    std::vector<std::future<bool>> futs;                                    
    for(size_t i = 0; i < num_examples; i++){                               
        futs.push_back(tp.enqueue(&kNN<N>::validate_N, this, i, K));       
    }                                                                       
    for(auto& fut : futs)                                                   
        if(fut.get()) right++;                                              

    return right / (double) num_examples;                                   
}      

bool validate_N(const size_t cmp, const size_t n) const{                    
    const size_t num_examples = labels.size();                              
    std::vector<char> dists(num_examples, -1);                              
    for(size_t i = 0; i < num_examples; i++){                               
        if(i == cmp) continue;                                              
        dists[i] = (bits[cmp] & bits[i]).count();                           

    }                                                                       
    typedef std::unordered_map<std::string,size_t> counter;                 
    counter counts;                                                         
    for(size_t i = 0; i < n; i++){                                          
        auto iter = std::max_element(dists.cbegin(), dists.cend());         
        size_t idx = std::distance(dists.cbegin(), iter);                   
        dists[idx] = -1; // Remove the top result.                          
        counts[labels[idx]] += 1;                                           
    }                                                                       
    auto iter = std::max_element(counts.cbegin(), counts.cend(),            
            [](const counter::value_type& a, const counter::value_type& b){ return a.second < b.second; }); 

    return labels[cmp] == iter->first;;                                     
}

编辑:

这就是我的想法。但是它非常慢。我不确定我是否做错了什么

template<size_t N>
struct popl 
{
    typedef unsigned long word_type;
    std::bitset<N> _cmp;

    popl(const std::bitset<N>& cmp) : _cmp(cmp) {}

    __device__
    int operator()(const std::bitset<N>& x) const
    {
        int pop_total = 0;
        #pragma unroll
        for(size_t i = 0; i < N/64; i++)
            pop_total += __popcll(x._M_w[i] & _cmp._M_w[i]);

        return pop_total;
    }
}; 

int main(void) {
    const size_t N = 2048;

    thrust::host_vector<std::bitset<N> > h_vec;
    load_bits(h_vec);

    thrust::device_vector<std::bitset<N> > d_vec = h_vec;
    thrust::device_vector<int> r_vec(h_vec.size(), 0);
    for(int i = 0; i < h_vec.size(); i++){
        r_vec[i] = thrust::transform_reduce(d_vec.cbegin(), d_vec.cend(),  popl<N>(d_vec[i]), 0, thrust::maximum<int>());
    }

    return 0;
}
rsaldnfx

rsaldnfx1#

CUDA具有用于32位和64位类型的填充计数内部函数。(__popc()__popcll()
这些可以直接在CUDA内核中使用,或者通过thrust(在函子中)传递给thrust::transform_reduce
如果这是您想在GPU上执行的唯一功能,则可能很难获得净“胜利”,因为将数据传输到GPU/从GPU传输数据的“成本”。(100000个位长为65536的向量),但根据我的计算,输出数据集的大小似乎是10- 40 GB(每个结果100000 * 100000 * 1-4字节)。
无论是CUDA内核还是推力函数和数据布局都应该精心设计,目标是让代码运行仅受内存带宽的限制。数据传输的成本也可以通过复制和计算操作的重叠来减轻,主要是在输出数据集上。
乍一看,这个问题似乎有点类似于计算向量集之间的欧氏距离的问题,因此从CUDA的Angular 来看,this question/answer可能会感兴趣。

**编辑:**添加一些我用来研究这个问题的代码。我能够获得比单纯的单线程CPU实现显著的加速(包括数据复制时间在内大约25倍),但我不知道使用“分块多线程和SSE/AVX“的CPU版本会有多快,所以看到更多的实现或获得一些性能数据会很有趣。我也不认为我这里的CUDA代码是高度优化的,它只是“第一次削减”。

在这个例子中,为了验证概念,我关注了一个小问题大小,N =2048,10000位集。对于这个小问题大小,我可以在共享内存中放置足够的位集向量,对于“小”线程块大小,以利用共享内存。因此,对于更大的N,必须修改这个特定的方法。

$ cat t581.cu
#include <iostream>
#include <vector>
#include <bitset>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>

#define nTPB 128
#define OUT_CHUNK 250
#define N_bits 2048
#define N_vecs 10000
const size_t N = N_bits;

__global__ void comp_dist(unsigned *in, unsigned *out, unsigned numvecs, unsigned start_idx, unsigned end_idx){
  __shared__ unsigned sdata[(N/32)*nTPB];
  int idx = threadIdx.x+blockDim.x*blockIdx.x;
  if (idx < numvecs)
    for (int i = 0; i < (N/32); i++)
      sdata[(i*nTPB)+threadIdx.x] = in[(i*numvecs)+idx];
  __syncthreads();
  int vidx = start_idx;
  if (idx < numvecs)
    while (vidx < end_idx) {
      unsigned sum = 0;
      for (int i = 0; i < N/32; i++)
        sum += __popc(sdata[(i*nTPB)+ threadIdx.x] & in[(i*numvecs)+vidx]);
      out[((vidx-start_idx)*numvecs)+idx] = sum;
      vidx++;}
}

void cpu_test(std::vector<std::bitset<N> > &in, std::vector<std::vector<unsigned> > &out){

  for (int i=0; i < in.size(); i++)
    for (int j=0; j< in.size(); j++)
      out[i][j] = (in[i] & in[j]).count();
}

int check_data(unsigned *d1, unsigned start_idx, std::vector<std::vector<unsigned> > &d2){
  for (int i = start_idx; i < start_idx+OUT_CHUNK; i++)
    for (int j = 0; j<N_vecs; j++)
      if (d1[((i-start_idx)*N_vecs)+j] != d2[i][j]) {std::cout << "mismatch at " << i << "," << j << " was: " << d1[((i-start_idx)*N_vecs)+j] << " should be: " << d2[i][j] << std::endl;  return 1;}
  return 0;
}

unsigned long long get_time_usec(){
  timeval tv;
  gettimeofday(&tv, 0);
  return (unsigned long long)(((unsigned long long)tv.tv_sec*1000000ULL)+(unsigned long long)tv.tv_usec);
}

int main(){

  unsigned long long t1, t2;
  std::vector<std::vector<unsigned> > distances;
  std::vector<std::bitset<N> > bits;

  for (int i = 0; i < N_vecs; i++){
    std::vector<unsigned> dist_row(N_vecs, 0);
    distances.push_back(dist_row);
    std::bitset<N> data;
    for (int j =0; j < N; j++) data[j] = rand() & 1;
    bits.push_back(data);}
  t1 = get_time_usec();
  cpu_test(bits, distances);
  t1 = get_time_usec() - t1;
  unsigned *h_data = new unsigned[(N/32)*N_vecs];
  memset(h_data, 0, (N/32)*N_vecs*sizeof(unsigned));
  for (int i = 0; i < N_vecs; i++)
    for (int j = 0; j < N; j++)
        if (bits[i][j]) h_data[(i)+((j/32)*N_vecs)] |= 1U<<(31-(j&31));

  unsigned *d_in, *d_out1, *d_out2, *h_out1, *h_out2;
  cudaMalloc(&d_in, (N/32)*N_vecs*sizeof(unsigned));
  cudaMalloc(&d_out1, N_vecs*OUT_CHUNK*sizeof(unsigned));
  cudaMalloc(&d_out2, N_vecs*OUT_CHUNK*sizeof(unsigned));
  cudaStream_t stream1, stream2;
  cudaStreamCreate(&stream1);
  cudaStreamCreate(&stream2);
  h_out1 = new unsigned[N_vecs*OUT_CHUNK];
  h_out2 = new unsigned[N_vecs*OUT_CHUNK];
  t2 = get_time_usec();
  cudaMemcpy(d_in, h_data, (N/32)*N_vecs*sizeof(unsigned), cudaMemcpyHostToDevice);
  for (int i = 0; i < N_vecs; i += 2*OUT_CHUNK){
    comp_dist<<<(N_vecs + nTPB - 1)/nTPB, nTPB, 0, stream1>>>(d_in, d_out1, N_vecs, i, i+OUT_CHUNK);
    cudaStreamSynchronize(stream2);
    if (i > 0) if (check_data(h_out2, i-OUT_CHUNK, distances)) return 1;
    comp_dist<<<(N_vecs + nTPB - 1)/nTPB, nTPB, 0, stream2>>>(d_in, d_out2, N_vecs, i+OUT_CHUNK, i+2*OUT_CHUNK);
    cudaMemcpyAsync(h_out1, d_out1, N_vecs*OUT_CHUNK*sizeof(unsigned), cudaMemcpyDeviceToHost, stream1);
    cudaMemcpyAsync(h_out2, d_out2, N_vecs*OUT_CHUNK*sizeof(unsigned), cudaMemcpyDeviceToHost, stream2);
    cudaStreamSynchronize(stream1);
    if (check_data(h_out1, i, distances)) return 1;
    }
  cudaDeviceSynchronize();
  t2 = get_time_usec() - t2;
  std::cout << "cpu time: " << ((float)t1)/(float)1000 << "ms gpu time: " << ((float) t2)/(float)1000 << "ms" << std::endl;
  return 0;
}
$ nvcc -O3 -arch=sm_20 -o t581 t581.cu
$ ./t581
cpu time: 20324.1ms gpu time: 753.76ms
$

CUDA 6.5、Fedora20、Xeon X5560、Quadro5000(cc2.0)GPU。上面的测试用例包括CPU和GPU产生的距离数据之间的结果验证。我还将其分解为具有结果数据传输的分块算法(和验证)与计算操作重叠,以使其更容易扩展到输出数据量非常大的情况(例如,100000位集)。

**编辑2:**这里是一个“windows版本”的代码:

#include <iostream>
#include <vector>
#include <bitset>
#include <stdlib.h>
#include <time.h>

#define nTPB 128
#define OUT_CHUNK 250
#define N_bits 2048
#define N_vecs 10000
const size_t N = N_bits;

#define cudaCheckErrors(msg) \
    do { \
        cudaError_t __err = cudaGetLastError(); \
        if (__err != cudaSuccess) { \
            fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                msg, cudaGetErrorString(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)


__global__ void comp_dist(unsigned *in, unsigned *out, unsigned numvecs, unsigned start_idx, unsigned end_idx){
  __shared__ unsigned sdata[(N/32)*nTPB];
  int idx = threadIdx.x+blockDim.x*blockIdx.x;
  if (idx < numvecs)
    for (int i = 0; i < (N/32); i++)
      sdata[(i*nTPB)+threadIdx.x] = in[(i*numvecs)+idx];
  __syncthreads();
  int vidx = start_idx;
  if (idx < numvecs)
    while (vidx < end_idx) {
      unsigned sum = 0;
      for (int i = 0; i < N/32; i++)
        sum += __popc(sdata[(i*nTPB)+ threadIdx.x] & in[(i*numvecs)+vidx]);
      out[((vidx-start_idx)*numvecs)+idx] = sum;
      vidx++;}
}

void cpu_test(std::vector<std::bitset<N> > &in, std::vector<std::vector<unsigned> > &out){

  for (unsigned i=0; i < in.size(); i++)
    for (unsigned j=0; j< in.size(); j++)
      out[i][j] = (in[i] & in[j]).count();
}

int check_data(unsigned *d1, unsigned start_idx, std::vector<std::vector<unsigned> > &d2){
  for (unsigned i = start_idx; i < start_idx+OUT_CHUNK; i++)
    for (unsigned j = 0; j<N_vecs; j++)
      if (d1[((i-start_idx)*N_vecs)+j] != d2[i][j]) {std::cout << "mismatch at " << i << "," << j << " was: " << d1[((i-start_idx)*N_vecs)+j] << " should be: " << d2[i][j] << std::endl;  return 1;}
  return 0;
}

unsigned long long get_time_usec(){

  return (unsigned long long)((clock()/(float)CLOCKS_PER_SEC)*(1000000ULL));
}

int main(){

  unsigned long long t1, t2;
  std::vector<std::vector<unsigned> > distances;
  std::vector<std::bitset<N> > bits;

  for (int i = 0; i < N_vecs; i++){
    std::vector<unsigned> dist_row(N_vecs, 0);
    distances.push_back(dist_row);
    std::bitset<N> data;
    for (int j =0; j < N; j++) data[j] = rand() & 1;
    bits.push_back(data);}
  t1 = get_time_usec();
  cpu_test(bits, distances);
  t1 = get_time_usec() - t1;
  unsigned *h_data = new unsigned[(N/32)*N_vecs];
  memset(h_data, 0, (N/32)*N_vecs*sizeof(unsigned));
  for (int i = 0; i < N_vecs; i++)
    for (int j = 0; j < N; j++)
        if (bits[i][j]) h_data[(i)+((j/32)*N_vecs)] |= 1U<<(31-(j&31));

  unsigned *d_in, *d_out1, *d_out2, *h_out1, *h_out2;
  cudaMalloc(&d_in, (N/32)*N_vecs*sizeof(unsigned));
  cudaMalloc(&d_out1, N_vecs*OUT_CHUNK*sizeof(unsigned));
  cudaMalloc(&d_out2, N_vecs*OUT_CHUNK*sizeof(unsigned));
  cudaCheckErrors("cudaMalloc fail");
  cudaStream_t stream1, stream2;
  cudaStreamCreate(&stream1);
  cudaStreamCreate(&stream2);
   cudaCheckErrors("cudaStrem fail");
  h_out1 = new unsigned[N_vecs*OUT_CHUNK];
  h_out2 = new unsigned[N_vecs*OUT_CHUNK];
  t2 = get_time_usec();
  cudaMemcpy(d_in, h_data, (N/32)*N_vecs*sizeof(unsigned), cudaMemcpyHostToDevice);
   cudaCheckErrors("cudaMemcpy fail");
  for (int i = 0; i < N_vecs; i += 2*OUT_CHUNK){
    comp_dist<<<(N_vecs + nTPB - 1)/nTPB, nTPB, 0, stream1>>>(d_in, d_out1, N_vecs, i, i+OUT_CHUNK);
    cudaCheckErrors("cuda kernel loop 1 fail");
    cudaStreamSynchronize(stream2);
    if (i > 0) if (check_data(h_out2, i-OUT_CHUNK, distances)) return 1;
    comp_dist<<<(N_vecs + nTPB - 1)/nTPB, nTPB, 0, stream2>>>(d_in, d_out2, N_vecs, i+OUT_CHUNK, i+2*OUT_CHUNK);
    cudaCheckErrors("cuda kernel loop 2 fail");
    cudaMemcpyAsync(h_out1, d_out1, N_vecs*OUT_CHUNK*sizeof(unsigned), cudaMemcpyDeviceToHost, stream1);
    cudaMemcpyAsync(h_out2, d_out2, N_vecs*OUT_CHUNK*sizeof(unsigned), cudaMemcpyDeviceToHost, stream2);
    cudaCheckErrors("cuda kernel loop 3 fail");
    cudaStreamSynchronize(stream1);
    if (check_data(h_out1, i, distances)) return 1;
    }
  cudaDeviceSynchronize();
  cudaCheckErrors("cuda kernel loop 4 fail");
  t2 = get_time_usec() - t2;
  std::cout << "cpu time: " << ((float)t1)/(float)1000 << "ms gpu time: " << ((float) t2)/(float)1000 << "ms" << std::endl;
  return 0;
}

我在这段代码中添加了CUDA错误检查。请确保在Visual Studio中构建一个 release 项目,而不是调试。当我在配备Quadro 1000 M GPU的Windows 7笔记本电脑上运行这段代码时,CPU执行时间约为35秒,GPU执行时间约为1.5秒。

ssm49v7z

ssm49v7z2#

OpenCL 1.2有popcount,它似乎可以做你想要的事情。它可以在矢量上工作,所以一次最多可以做ulong16,也就是1024位。请注意,NVIDIA驱动程序只支持OpenCL 1.1,它不包括这个功能。
当然,您可以使用函数或表来快速计算它,因此OpenCL 1.1实现也是可能的,并且可能会在设备的内存带宽上运行。

相关问题