理解 CUDA、Numba、Cupy 等的扩展示例

Posted 2023-03-12

【中文标题】理解 CUDA、Numba、Cupy 等的扩展示例

【英文标题】：Extended example to understand CUDA, Numba, Cupy, etc

【发布时间】：2020-12-16 16:39:16

【问题描述】：

大多数在线可用的 Numba、CuPy 等示例都是简单的数组添加，显示了从 cpu 单核/线程到 gpu 的加速。并且命令文档大多缺乏好的例子。这篇文章旨在提供一个更全面的示例。

提供初始代码here。它是经典元胞自动机的简单模型。最初，它甚至不使用 numpy，只是使用普通 python 和 Pyglet 模块进行可视化。

我的目标是将此代码扩展到特定问题（这将非常大），但首先我认为最好已经针对 GPU 使用进行优化。

game_of_life.py 是这样的：

import random as rnd
import pyglet
#import numpy as np
#from numba import vectorize, cuda, jit

class GameOfLife: 
 
    def __init__(self, window_width, window_height, cell_size, percent_fill):
        self.grid_width = int(window_width / cell_size) # cell_size 
        self.grid_height = int(window_height / cell_size) # 
        self.cell_size = cell_size
        self.percent_fill = percent_fill
        self.cells = []
        self.generate_cells()
  
    def generate_cells(self):
        for row in range(0, self.grid_height): 
            self.cells.append([])
            for col in range(0, self.grid_width):
                if rnd.random() < self.percent_fill:
                    self.cells[row].append(1)
                else:
                    self.cells[row].append(0)
                
    def run_rules(self): 
        temp = []
        for row in range(0, self.grid_height):
            temp.append([])
            for col in range(0, self.grid_width):
                cell_sum = sum([self.get_cell_value(row - 1, col),
                                self.get_cell_value(row - 1, col - 1),
                                self.get_cell_value(row,     col - 1),
                                self.get_cell_value(row + 1, col - 1),
                                self.get_cell_value(row + 1, col),
                                self.get_cell_value(row + 1, col + 1),
                                self.get_cell_value(row,     col + 1),
                                self.get_cell_value(row - 1, col + 1)])
                
                if self.cells[row][col] == 0 and cell_sum == 3:
                    temp[row].append(1)
                elif self.cells[row][col] == 1 and (cell_sum == 3 or cell_sum == 2):
                    temp[row].append(1)
                else:                 
                    temp[row].append(0)
        
        self.cells = temp

    def get_cell_value(self, row, col): 
        if row >= 0 and row < self.grid_height and col >= 0 and col < self.grid_width:
           return self.cells[row][col]
        return 0

    def draw(self): 
        for row in range(0, self.grid_height):
            for col in range(0, self.grid_width):
                if self.cells[row][col] == 1:
                    #(0, 0) (0, 20) (20, 0) (20, 20)
                    square_coords = (row * self.cell_size,                  col * self.cell_size,
                                     row * self.cell_size,                  col * self.cell_size + self.cell_size,
                                     row * self.cell_size + self.cell_size, col * self.cell_size,
                                     row * self.cell_size + self.cell_size, col * self.cell_size + self.cell_size)
                    pyglet.graphics.draw_indexed(4, pyglet.gl.GL_TRIANGLES,
                                         [0, 1, 2, 1, 2, 3],
                                         ('v2i', square_coords))

首先，我可以在generate_cells 这个self.cells = np.asarray(self.cells) 的末尾和run_rules 这个self.cells = np.asarray(temp) 的末尾使用numpy 添加，因为之前这样做不会带来加速，如here 所示。（实际上更改为 numpy 并没有带来明显的加速）

例如，关于gpu，我在每个函数之前添加了@jit，并且变得非常慢。 也尝试使用@vectorize(['float32(float32, float32)'], target='cuda')，但这提出了一个问题：如何在只有self 作为输入参数的函数中使用@vectorize？

我也试过用numpy代替cupy，比如self.cells = cupy.asarray(self.cells)，但也变得很慢。

按照 gpu 使用扩展示例的初步想法，解决问题的正确方法是什么？放置修改/矢量化/并行化/numba/cupy等的正确位置在哪里？最重要的是，为什么？

附加信息：除了提供的代码，这里是 main.py 文件：

import pyglet
from game_of_life import GameOfLife 
 
class Window(pyglet.window.Window):
 
    def __init__(self):
        super().__init__(800,800)
        self.gameOfLife = GameOfLife(self.get_size()[0],
                                     self.get_size()[1],
                                     15,  # the lesser this value, more computation intensive will be
                                     0.5) 

        pyglet.clock.schedule_interval(self.update, 1.0/24.0) # 24 frames per second
 
    def on_draw(self):
        self.clear()
        self.gameOfLife.draw()
        
    def update(self, dt):
        self.gameOfLife.run_rules()
 
if __name__ == '__main__':
    window = Window()
    pyglet.app.run()

【问题讨论】：

我对使用 cuda.jit 装饰器的理解非常有限，但在我看来，这种内核性能不佳的主要原因是在 CPU 和 GPU 之间传输过多数据时。为避免这种情况，必须只传递必要的变量，尤其是在谈论大型数组时。我认为通过使用 self 作为每个函数（将是内核）的参数，您可能会传递不必要的数据。另外，请记住，每个线程都对数组的单个元素进行操作，因此使用 for 迭代数组将不会被并行化。希望这会有所帮助。

@boi ，感谢您指出这一点。我在 3 个月前开始使用 Python，这是我使用的第一种语言。我从来没有使用过，它对我来说是新的，即使我编码了 +10 年。 self、_init_ 等对我来说是新的。我会更仔细地查看以正确传递参数。关于for，不知道Python有没有类似parfor的东西，比如Matlab？

实际上是的，numba.prange 可能就是您要找的东西，尽管我认为不可能在 numba.cuda 中并行化循环。这是文档：numba.readthedocs.io/en/stable/user/…。我对这一切也很陌生:)。

@rod_CAE 对此有何更新？我还可以欣赏class 构造的新颖性。

@Sterling 不幸的是没有。该项目在没有并行化的情况下继续进行了一段时间，但由于我的其他项目变得比这个更重要，它目前被搁置了......

【参考方案1】：

我不太了解您的示例，但我只需要 GPU 计算。痛了几天，大概明白它的用法了，给大家演示一下，希望对大家有所帮助。 另外需要指出的是，在使用“...kernel(cuts,cuts”的时候，我会放两个。因为第一个在传入的时候指定了类型，所以会被核心用作遍历元素，不能被索引读取。所以我用第二个来计算空闲索引数据。

```
binsort_kernel = cp.ElementwiseKernel(
'int32 I,raw T cut,raw T ind,int32 row,int32 col,int32 q','raw T out,raw T bin,raw T num',    
'''
int i_x = i / col;                
int i_y = i % col;                
int b_f = i_x*col;                
int b_l = b_f+col;                
int n_x = i_x * q;                
int inx = i_x%row*col;            
////////////////////////////////////////////////////////////////////////////////////////
int r_x = 0; int adi = 0; int adb = 0;  
////////////////////////////////////////////////////////////////////////////////////////
if (i_y == 0)

for(size_t j=b_f; j<b_l; j++)
    if (cut[j]<q)                
        r_x = inx + j -b_f;       
        adb = n_x + cut[j];       
        adi = bin[adb] + num[adb];
        out[adi] = ind[r_x];      
        num[adb]+= 1;             
    

////////////////////////////////////////////////////////////////////////////////////////
''','binsort')

binsort_kernel(cuts,cuts,ind,row,col,q,iout,bins,bnum)

【讨论】：

这个答案解决了什么问题？ OP 的帖子中似乎没有使用 binsort，但看起来您正在使用 CuPy 并展示了如何将 ElementwiseKernel 与 C 代码一起使用？