logo标签列表

在GPU上使用popcnt

c++cudaopenclgputhrust
3
我需要计算一个大集合（>10000）的位向量（std::bitset），其中N从2^10到2^16不等，代码如下：
(a & b).count()
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 760m）上计算“(a & b).count()”函数？理想情况下，我只需传递两个大小为N的内存块。我正在考虑使用Thrust，但找不到“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;
}
似乎 char 不足以容纳 std::bitset<N> 的人口统计数据，其中 N 为1024到65536？您是否对按位与的结果做出了一些假设？即使它可以，我猜您正在生成大约10GB的数据（对于 bits(100000) 中的 distances）？ - Robert Crovella
@RobertCrovella 位集合非常稀疏。我可能不会在我的测试中使用2 ^ 16，但我想尝试2 ^ 14。2 ^ 12可能是最有可能的用例。此外，我不需要存储“距离”，那只是一个示例/简化。我只需要找到每个“bits[i]”中的“K”个最大距离，因此大部分数据都被丢弃了。我正在使用这种度量计算最近邻居。 - en4bz
2个回答
10
CUDA有种群计数内建函数，适用于32位和64位类型。(__popc() 和 __popcll())。
这些函数可以直接在CUDA内核中使用，或通过thrust（在函数对象中）传递给thrust::transform_reduce。
如果这是您想要在GPU上执行的唯一功能，由于数据传输的“成本”，可能很难获得净“胜利”。您的整体输入数据集似乎大小约为1GB（100000个长度为65536位的向量），但根据我的计算，输出数据集的大小似乎为10-40GB（100000 * 100000 * 1-4字节每个结果）。
CUDA核函数或thrust函数和数据布局应该精心设计，以便让代码仅受内存带宽限制。通过重叠复制和计算操作（主要针对输出数据集），数据传输的成本也可以得到缓解，可能会在很大程度上得到缓解。
乍一看，这个问题似乎与计算向量集之间的欧几里得距离的问题有些相似，因此从CUDA的角度来看，这个问题/答案可能会引起兴趣。 编辑：添加了一些我用来调查这个问题的代码。我能够获得一个显著的加速（包括数据复制时间在内的约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中构建一个“发行版”项目，而不是调试版。当我在配备Quadro1000M GPU的Windows 7笔记本电脑上运行此代码时，CPU执行时间约为35秒，GPU执行时间约为1.5秒。
我试图使用给定的度量来找到K个最近邻居。这意味着我不需要存储“距离”，只需存储K个最大值（其中K < 50）。另外，位集非常稀疏。我想将bits加载到GPU上，然后传输回一个std :: vector <char>以计算每个bits[i]的K个最大元素（及其索引）。我有3GB的VRAM和16GB的RAM。你认为将100000个std :: vector <char>传输回CPU进行处理会消除我在GPU上获得的任何速度提升吗？ - en4bz
1
我会尽可能将算法的大部分移动到GPU上。例如，查找K个最大元素。这将减少数据移动的成本。但是，作为一个开始，您可以只做您在此处展示的内容。在这种情况下，我会专注于复制和计算的重叠——由于输出数据集的大小，这几乎是必要的——并查看比较结果如何。您是否已经具体制定了当前CPU案例，并进行了性能测量？ - Robert Crovella
请看我的修改。我尝试了自己的实现，但速度极慢。由于某种原因，我在系统调用中花费了约1/4的时间（约7分钟），总共24分钟。通过CPU大约需要3分钟。在此期间，我会尝试你的方法。感谢你的帮助！ - en4bz
在我的代码中，当我将N_vecs增加到100000时，CPU时间约为40分钟，GPU时间约为45秒，因此速度提高了约50倍。知道你的CPU实现是什么样子的会很好，至少要知道与3分钟测量相对应的数据大小。这是针对N=2048吗？还是针对100000个位集？ - Robert Crovella
@RobertCorvella 我已经添加了CPU实现的核心部分。我已经切换到线程池以最大化吞吐量。我认为3分钟是使用bitset<4096>和K ~ 15。位集的数量为100000。我尝试了你的代码，但它的检查函数失败了。我认为这与你所假设的大小有关。std::bitset使用无符号长整型作为其内部表示。我看到你正在使用unsigned并且经常除以32。 - en4bz
我认为失败的原因与 bitset 的内部表示没有任何关系。我的所有访问方法都不会假设关于内部表示的任何内容。我每次只访问一个位。为了简明起见，我省略了 CUDA 错误检查。您可以使用 cuda-memcheck 运行我的代码，看看您的机器上是否报告任何错误。您是在 Windows 笔记本电脑上运行吗？并且您是按原样运行我的代码还是更改了什么？ - Robert Crovella
1
OpenCL 1.2有popcount函数，似乎可以实现您想要的功能。它可以在向量上工作，因此每次最多处理1024位，即ulong16。请注意，NVIDIA驱动程序仅支持不包括此函数的OpenCL 1.1版本。
当然，您也可以使用函数或表格来快速计算它，因此OpenCL 1.1的实现也是可行的，并且可能以设备的内存带宽运行。
1
对于NVIDIA案例，使用小位内联汇编和popc PTX指令也非常简单。请参见：https://github.com/kylelutz/compute/blob/master/include/boost/compute/functional/detail/nvidia_popcount.hpp - Kyle Lutz
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接