如何在 cuda 推力中实现嵌套循环

Posted 2023-02-22

【中文标题】如何在 cuda 推力中实现嵌套循环

【英文标题】：How to implement nested loops in cuda thrust

【发布时间】：2015-03-28 16:22:51

【问题描述】：

我目前必须按如下方式运行嵌套循环：

for(int i = 0; i < N; i++)
    for(int j = i+1; j <= N; j++)
        compute(...)//some calculation here
    

我尝试将第一个循环留在CPU 中，并在GPU 中执行第二个循环。结果是too many memory access。还有其他方法吗？例如thrust::reduce_by_key?

整个程序都在这里：

#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/generate.h>
#include <thrust/sort.h>
#include <thrust/binary_search.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/random.h>
#include <cmath>
#include <iostream>
#include <iomanip>

#define N 1000000

// define a 2d point pair
typedef thrust::tuple<float, float> Point;

// return a random Point in [0,1)^2
Point make_point(void)

  static thrust::default_random_engine rng(12345);
  static thrust::uniform_real_distribution<float> dist(0.0f, 1.0f);
  float x = dist(rng);
  float y = dist(rng);
  return Point(x,y);


struct sqrt_dis: public thrust::unary_function<Point, double>

  float x, y;
  double tmp;
  sqrt_dis(float _x, float _y): x(_x), y(_y)
  __host__ __device__
  float operator()(Point a)
 
    tmp =(thrust::get<0>(a)-x)*(thrust::get<0>(a)-x)+\
    (thrust::get<1>(a)-y)*(thrust::get<1>(a)-y);
    tmp = -1.0*(sqrt(tmp));
    return (1.0/tmp);
 

;


int main(void) 
  clock_t t1, t2;
  double result;

  t1 = clock();
  // allocate some random points in the unit square on the host
  thrust::host_vector<Point> h_points(N);
  thrust::generate(h_points.begin(), h_points.end(), make_point);

  // transfer to device
  thrust::device_vector<Point> points = h_points;

  thrust::plus<double> binary_op;
  float init = 0;

  for(int i = 0; i < N; i++)
    Point tmp_i = points[i];
    float x = thrust::get<0>(tmp_i);
    float y = thrust::get<1>(tmp_i);
    result += thrust::transform_reduce(points.begin()+i,\
                                       points.end(),sqrt_dis(x,y),\
                                       init,binary_op);
    std::cout<<"result"<<i<<": "<<result<<std::endl;
  
  t2 = clock()-t1;
  std::cout<<"result: ";
  std::cout.precision(10);
  std::cout<< result <<std::endl;
  std::cout<<"run time: "<<t2/CLOCKS_PER_SEC<<"s"<<std::endl;
  return 0;
 

【问题讨论】：

您需要result 中的总和还是需要每次循环的迭代总和？

结果的总和

【参考方案1】：

编辑：既然您已经发布了一个示例，那么您可以如何解决它：

您将n 2D 点存储在这样的线性数组中（此处为 n=4）

points = [p0 p1 p2 p3]

根据您的代码，我假设您要计算：

result = f(p0, p1) + f(p0, p2) + f(p0, p3) +
         f(p1, p2) + f(p1, p3) +
         f(p2, p3)

f() 是你需要执行的距离函数m times in total:

m = (n-1)*n/2

在本例中：m=6

你可以把这个问题看成一个三角矩阵：

[ p0 p1 p2 p3 ] 
[    p1 p2 p3 ]
[       p2 p3 ]
[          p3 ]

将此矩阵转换为包含m 元素的线性向量，同时省略对角线元素的结果是：

[p1 p2 p3 p2 p3 p3]

向量中元素的索引是k = [0,m-1]。 索引k 可以是remapped to columns and rows of the triangular matrix 到k -> (i,j)：

i = n - 2 - floor(sqrt(-8*k + 4*n*(n-1)-7)/2.0 - 0.5)
j = k + i + 1 - n*(n-1)/2 + (n-i)*((n-i)-1)/2

i 是行，j 是列。

在我们的示例中：

0 -> (0, 1)
1 -> (0, 2)
2 -> (0, 3)
3 -> (1, 2)
4 -> (1, 3)
5 -> (2, 3)

现在您可以将所有这些放在一起并执行修改后的距离函子 m 次，它应用上述映射来根据索引获取相应的对，然后总结所有内容。

我相应地修改了您的代码：

#include <thrust/device_vector.h>
#include <thrust/generate.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/transform_reduce.h>
#include <thrust/random.h>
#include <math.h>
#include <iostream>
#include <stdio.h>
#include <stdint.h>

#define PRINT_DEBUG

typedef float Float;

// define a 2d point pair
typedef thrust::tuple<Float, Float> Point;

// return a random Point in [0,1)^2
Point make_point(void)

    static thrust::default_random_engine rng(12345);
    static thrust::uniform_real_distribution<Float> dist(0.0, 1.0);
    Float x = dist(rng);
    Float y = dist(rng);
    return Point(x,y);



struct sqrt_dis_new

    typedef thrust::device_ptr<Point> DevPtr;

    DevPtr points;
    const uint64_t n;

    __host__
    sqrt_dis_new(uint64_t n, DevPtr p) : n(n), points(p)
    
    

    __device__ 
    Float operator()(uint64_t k) const
    
        // calculate indices in triangular matrix
        const uint64_t i = n - 2 - floor(sqrt((double)(-8*k + 4*n*(n-1)-7))/2.0 - 0.5);
        const uint64_t j = k + i + 1 - n*(n-1)/2 + (n-i)*((n-i)-1)/2;

#ifdef PRINT_DEBUG
        printf("%llu -> (%llu, %llu)\n", k,i,j);
#endif

        const Point& p1 = *(points.get()+j);
        const Point& p2 = *(points.get()+i);

        const Float xm = thrust::get<0>(p1)-thrust::get<0>(p2);
        const Float ym = thrust::get<1>(p1)-thrust::get<1>(p2);

        return 1.0/(-1.0 * sqrt(xm*xm + ym*ym));
    
;


int main()

    const uint64_t N = 4;

    // allocate some random points in the unit square on the host
    thrust::host_vector<Point> h_points(N);
    thrust::generate(h_points.begin(), h_points.end(), make_point);

    // transfer to device
    thrust::device_vector<Point> d_points = h_points;

    const uint64_t count = (N-1)*N/2;

    std::cout << count << std::endl;


    thrust::plus<Float> binary_op;
    const Float init = 0.0;

    Float result = thrust::transform_reduce(thrust::make_counting_iterator((uint64_t)0),
                                            thrust::make_counting_iterator(count),
                                            sqrt_dis_new(N, d_points.data()),
                                            init,
                                            binary_op);

    std::cout.precision(10);  
    std::cout<<"result: " << result << std::endl;

    return 0;
    

---

这取决于您未指定的 compute 函数。 通常，如果计算是独立的，则对于 i 和 j 的每个组合，您都可以展开循环并以 2D 方式启动内核。 查看 Thrust 示例并找出与您的问题类似的用例。

【讨论】：

我已经发布了我当前的版本。