c++ - 将两个不同长度的 vector 根据各自包含全局地址的索引 vector 拼接到具有推力的共同长度的新 vector

Question

这个问题困扰我好几年了。我从这个论坛学习了大量的 c++ 和 cuda。以前我用 fortran 串行代码编写了以下带有很多条件语句的代码，并使用了 gotos 因为我找不到一个聪明的方法来做到这一点。
问题就在这里。

给定 4 个 vector :

int indx(nshape);
float dnx(nshape); 
/* nshape > nord */
int indy(nord);
float dny(nord);

indx 和 indy 是包含全局坐标的索引 vector (分别为值 dnx、dny 的键)。在解析到这个所需的交错/拼接函数之前，它们的全局范围是未知的。所知道的是可能的局部范围的长度可以是 [0,nord*nord] 以及 vector indx 和 indy 中的最大值和最小值。

我想创建相同长度的新 vector dn1 和 dn2，其中包含 dnx 和 dny 的原始值，但被扩展为用零填充原始 vector dnx 和 dny，用于它们不包含的所有全局坐标另一个 vector 。它们将为需要对齐全局地址的外积形成 vector 。

我无法在网上找到任何关于在 c++ 中使用逻辑掩码(如在 fortran 中)进行并行化的引用资料。我的出发点是使用推力库 stable_sort 进行升序排序，使用 binary_search 比较数组、分区等。也许有一种清晰简洁的方法可以做到这一点。

下面的示例索引和值 vec 通常不会从 0 开始或与临时索引 vector 的本地寻址一致，也不会出现任何偶数模式 - 这些值只是为了帮助说明。)

indx[]={0,2,4,6,8,10,12};  indy[]={1, 2, 3, 4};
dnx[]={99,99,99,99,99,99,99};  dny[]={66,66,66,66};

ind[]={0,1,2,3,4,6,8,10,12}
dn1[]={99,0,99,0,99,99,99,99,99}
dn2[]={0,66,66,66,66,0,0,0,0}

之前我做了类似下面的事情，其中​​内核根据以下条件应用比较、填充和流，并继续返回通过这些条件行之一再次输入，直到最大局部索引超过最大 vector i 的长度， e i , j > nshape :

3
if(indx[i] < indy[j]{kernel_1; i++; if(i > nshape){return}; goto 3}
if(indx[i] == indy[j]){kernel_2;i++;j++; if(i || j > nshape) {return}; goto 3}
if(indx[i] > indy[j]{kernel_3, j++, if(j>nshape){return}; goto 3}

对于杂种伪代码感到抱歉。我真的很期待使用 c++、cuda、thrust 的任何想法或更好的解决方案。非常感谢。标记

Answer 1

我想到的一种方法是进行并行二进制搜索，从全局索引 vector 中获取每个值，然后查看它是否在键 vector 中具有匹配项。如果它在键 vector 中有匹配项，则将相应的值放入结果中。如果它在键 vector 中没有匹配项，则将 0 放入结果中。

所以对于每个位置:

ind[]={0,1,2,3,4,6,8,10,12}

查看是否有匹配的索引:

indy[]={1, 2, 3, 4};

我们将使用并行二进制搜索来执行此操作，并返回相应的索引(indy 中的匹配值)。

如果我们确实找到匹配项，那么在结果中的相关位置，我们将放置对应于 indy[matching_index] 的值，即。 dny[匹配索引]。否则在结果中置零。

在推力的情况下，我能够将其减少为两次推力调用。

第一个是 thrust::lower_bound 操作，它实际上是一个矢量化/并行二进制搜索。就像在 CUDA 案例中一样，我们使用二进制搜索来获取全局 vector (ind) 的每个元素，并查看关键 vector 中是否存在匹配项(例如 indx >), 返回键 vector 中匹配位置的索引(下界)。

第二个调用是对 thrust::for_each 的复杂使用。我们为 for_each 操作创建了一个特殊的仿函数 (extend_functor)，它使用指向键 vector 开始的指针进行初始化(例如 indx )，它的长度，以及指向值 vector 的指针(例如 dnx)。 extend_functor 然后采用全局 vector 、下界 vector 和结果 vector 值的 3 元组，并执行剩余的步骤。如果下限值在键 vector 的长度内，则检查下限是否在键 vector 和全局 vector 之间产生匹配。如果是，则将相应的值放入结果 vector ，否则将零放入结果 vector 。

以下代码使用 CUDA 和 thrust 实现了这一点。

#include <iostream>
#include <thrust/device_vector.h>
#include <thrust/binary_search.h>
#include <thrust/copy.h>

#define MAX_DSIZE 32768
#define nTPB 256


struct extend_functor
{
  int vec_len;
  int *vec1;
  float *vec2;
  extend_functor(int *_vec1, int _vec_len, float *_vec2) : vec1(_vec1), vec2(_vec2),  vec_len(_vec_len) {};
  template <typename Tuple>
  __device__ __host__
  void operator()(const Tuple &my_tuple) {
    float result = 0.0f;
    if (thrust::get<1>(my_tuple) < vec_len)
      if (thrust::get<0>(my_tuple) == vec1[thrust::get<1>(my_tuple)]) result = vec2[thrust::get<1>(my_tuple)];
    thrust::get<2>(my_tuple) = result;
  }

};

// binary search, looking for key in a

__device__ void bsearch_range(const int *a, const int key, const unsigned len_a, unsigned *idx){
  unsigned lower = 0;
  unsigned upper = len_a;
  unsigned midpt;
  while (lower < upper){
    midpt = (lower + upper)>>1;
    if (a[midpt] < key) lower = midpt +1;
    else upper = midpt;
    }
  *idx = lower;
  return;
  }


// k is the key vector
// g is the global index vector
// v is the value vector
// r is the result vector

__global__ void extend_kernel(const int *k, const unsigned len_k, const int *g,  const unsigned len_g, const float *v, float *r){
  unsigned idx = (blockDim.x * blockIdx.x) + threadIdx.x;
  if (idx < len_g) {
    unsigned my_idx;
    int g_key = g[idx];
    bsearch_range(k, g_key, len_k, &my_idx);
    int my_key = -1;
    if (my_idx < len_k)
      my_key = k[my_idx];
    float my_val;
    if (g_key == my_key) my_val = v[my_idx];
    else my_val = 0.0f;
    r[idx] = my_val;
    }
}


int main(){

  int len_x = 7;
  int len_y = 4;
  int len_g = 9;
  int indx[]={0,2,4,6,8,10,12};
  int indy[]={1, 2, 3, 4};
  float dnx[]={91.0f,92.0f,93.0f,94.0f,95.0f,96.0f,97.0f};
  float dny[]={61.0f,62.0f,63.0f,64.0f};
  int ind[]={0,1,2,3,4,6,8,10,12};

  int *h_k, *d_k, *h_g, *d_g;
  float *h_v, *d_v, *h_r, *d_r;

  h_k = (int *)malloc(MAX_DSIZE*sizeof(int));
  h_g = (int *)malloc(MAX_DSIZE*sizeof(int));
  h_v = (float *)malloc(MAX_DSIZE*sizeof(float));
  h_r = (float *)malloc(MAX_DSIZE*sizeof(float));

  cudaMalloc(&d_k, MAX_DSIZE*sizeof(int));
  cudaMalloc(&d_g, MAX_DSIZE*sizeof(int));
  cudaMalloc(&d_v, MAX_DSIZE*sizeof(float));
  cudaMalloc(&d_r, MAX_DSIZE*sizeof(float));

// test case x

  memcpy(h_k, indx, len_x*sizeof(int));
  memcpy(h_g, ind, len_g*sizeof(int));
  memcpy(h_v, dnx, len_x*sizeof(float));

  cudaMemcpy(d_k, h_k, len_x*sizeof(int), cudaMemcpyHostToDevice);
  cudaMemcpy(d_g, h_g, len_g*sizeof(int), cudaMemcpyHostToDevice);
  cudaMemcpy(d_v, h_v, len_x*sizeof(float), cudaMemcpyHostToDevice);
  extend_kernel<<<(len_g+nTPB-1)/nTPB, nTPB >>>(d_k, len_x, d_g, len_g, d_v, d_r);
  cudaMemcpy(h_r, d_r, len_g*sizeof(float), cudaMemcpyDeviceToHost);

  std::cout << "case x result: ";
  for (int i=0; i < len_g; i++)
    std::cout << h_r[i] << " ";
  std::cout << std::endl;


// test case y

  memcpy(h_k, indy, len_y*sizeof(int));
  memcpy(h_g, ind, len_g*sizeof(int));
  memcpy(h_v, dny, len_y*sizeof(float));

  cudaMemcpy(d_k, h_k, len_y*sizeof(int), cudaMemcpyHostToDevice);
  cudaMemcpy(d_g, h_g, len_g*sizeof(int), cudaMemcpyHostToDevice);
  cudaMemcpy(d_v, h_v, len_y*sizeof(float), cudaMemcpyHostToDevice);
  extend_kernel<<<(len_g+nTPB-1)/nTPB, nTPB >>>(d_k, len_y, d_g, len_g, d_v, d_r);
  cudaMemcpy(h_r, d_r, len_g*sizeof(float), cudaMemcpyDeviceToHost);

  std::cout << "case y result: ";
  for (int i=0; i < len_g; i++)
    std::cout << h_r[i] << " ";
  std::cout << std::endl;

// using thrust

  thrust::device_vector<int> tind(ind, ind+len_g);
  thrust::device_vector<int> tindx(indx, indx+len_x);
  thrust::device_vector<float> tdnx(dnx, dnx+len_x);
  thrust::device_vector<float> tresult(len_g);

  thrust::device_vector<int> tbound(len_g);
  thrust::lower_bound(tindx.begin(), tindx.end(), tind.begin(), tind.end(), tbound.begin());
  thrust::for_each(thrust::make_zip_iterator(thrust::make_tuple(tind.begin(), tbound.begin(), tresult.begin())), thrust::make_zip_iterator(thrust::make_tuple(tind.end(), tbound.end(), tresult.end())), extend_functor(thrust::raw_pointer_cast(tindx.data()), len_x, thrust::raw_pointer_cast(tdnx.data())));

  std::cout << "thrust case x result: ";
  thrust::copy(tresult.begin(), tresult.end(), std::ostream_iterator<float>(std::cout, " "));
  std::cout << std::endl;
  return 0;
}