pycuda编程中线程、线程块与网格

持续创作，加速成长！这是我参与「掘金日新计划 · 10 月更文挑战」的第22天，点击查看活动详情

首先说一下线程，线程在这里的含义指的是在GPU内核(单个内核)上运行的一系列指令。在GPU上，多个线程的执行是以线程块(抽象单位)为单位进行。在前面的博客中有一行代码int i = threadId.x，这行代码的作用是获取线程的ID，后面的x表示维度，除了x我们还可以在y和z两个维度进行索引。在GPU编程这一块，线程块进一步可以组成网格，就如前面博客中使用的这一张图。

所以我们可以使用int i = ｂｌｏｃｋId.x，int i = ｂｌｏｃｋId.ｙ，int i = ｂｌｏｃｋId.ｚ三个维度来索引线程块。

所以，假设我们在一张二维的图像（灰度图）中索引线程时，此时我们将我们就必须按照下面的格式进行了，因为是二维，所以要计算两个方向的索引。再就是CUDA对于线程块的带下有一定的限制，一般来说一个线程块中最多含有１０２４个线程，这样的话对于二维图像，最多可以做出尺寸为３２ｘ３２＝１０２４大小的图像，这样的话作用太小了，像现在的图像动不动就是１９２０ｘ１０２４，完全不够用的，所以为了解决这个问题我们必须，在一个网格中的多个线程块上运行我们编写好的核函数，所以线程的索引必须覆盖到多个线程块中。

就如下面的代码，blockDim指的是线程块的大小，就像下面那张丑陋的图里面画的一样，整个这九个黑色的格子我们认为它是ｇｒｉｄ，其中一个黑色的格子就是ｂｌｏｃｋ，一个ｂｌｏｃｋ里面有九个线程，只用一个线程块（０，０）表示一个１９２０ｘ１２００的图像显然是不够的，所以我们就要用很多个ｂｌｏｃｋ一起来覆盖表示这个图像，因此，此时ｔｈｒｅａｄＩｄｘ就不再单纯的用一个块来表示，假设一个ｂｌｏｃｋ大小为ｂｌｏｃｋＤｉｍ，所以ｂｌｏｃｋＩｄｘ＊ｂｌｏｃｋＤｉｍ就是一块中ｔｈｒｅａｄ的数目，至于具体一共需要多少个，后面我们后面在实例化核函数时指定ｇｒｉｄ和ｂｌｏｃｋ大小即可。

#define _X (threadIdx.x + blockIdx.x * blockDim.x)
#define _Y (threadIdx.y + blockIdx.y * blockDim.y)

下面写一个核函数，进行举例

kernal = SourceModule(
    """
    #define _X (threadIdx.x + blockIdx.x * blockDim.x)
    #define _Y (threadIdx.y + blockIdx.y * blockDim.y)
    #define _WIDTH (blockDim.x * gridDim.x)
    #define _HEIGHT (blockDim.y * gridDim.y)
    #define _XM(x) ( (x + _WIDTH) % _WIDTH )
    #define _YM(y) ( (y + _HEIGHT) % _HEIGHT )
    #define _INDEX(x,y) (_XM(x) + _YM(y) * _WIDTH)
    
    __device__ int nbrs(int x, int y, int *in)
    {
        return ( in[ _INDEX(x-1 ,y+1) ] + in[ _INDEX(x-1 ,y) ] + in[ _INDEX(x-1 ,y-1) ] + in[ _INDEX(x ,y+1) ] + in[ _INDEX(x ,y-1) ] + in[ _INDEX(x+1 ,y+1) ] + in[ _INDEX(x+1 ,y) ] + in[ _INDEX(x+1 ,y-1) ]);
    }
    __global__ void conway_ker(int *lattice_out, int  *lattice)
    {
        int x = _X, y = _Y;
        int n = nbrs(x, y, lattice);
        if( lattice[_INDEX(x,y)] == 1 )
            switch(n)
            {
                case 2:
                case 3:lattice_out[ _INDEX(x,y) ] =1;break;
                default: lattice_out[ _INDEX(x,y) ] =0;
            }
        else if( lattice[ _INDEX(x,y) ] == 0 )
            switch(n)
            {
                case 3:lattice_out[ _INDEX(x,y) ] =1;break;
                default: lattice_out[ _INDEX(x,y) ] =0;
            }
    }
    """
)

在指定我们需要的大小时候我们可以进行下面的操作，这里的N是２５６，不一定非要设置为Ｎ＊Ｎ，Ｎ＊Ｍ也是可以的，都没有问题，上面程序的输出结果我就用Ｍ＊Ｎ进行举例。

conway_ker = kernal.get_function("conway_ker")

def update_gpu(frameNum, img, newLattice_gpu, lattice_gpu, N):
    conway_ker(newLattice_gpu, lattice_gpu, grid=(N//32, N//32, 1), block=(32, 32, 1))
    img.set_data(newLattice_gpu.get())
    lattice_gpu[:] = newLattice_gpu[:]
    return img


if __name__ == '__main__':
    N = 256
    lattice = np.int32(np.random.choice([1, 0], N * N, p=[0.25, 0.75]).reshape(N, N) )
    lattice_gpu = gpuarray.to_gpu(lattice)
    newLattice_gpu = gpuarray.empty_like(lattice_gpu)
    fig, ax = plt.subplots()
    img = ax.imshow(lattice_gpu.get(), interpolation='nearest')
    ani = animation.FuncAnimation(fig, update_gpu, fargs=(img, newLattice_gpu, lattice_gpu, N,), interval=0, frames=1000, save_count=1000)
    plt.show()

注意ｇｒｉｄ不要设置的太大，不然会报错，报错的结果可能如下，不要慌，这是因为设置的ｇｒｉｄ太大了，本身就没有那么多人，你非得安排那么多人，肯定不够呀。

最后把需要的头文件放一下，autoinit这个一定加上，不然可能会出错，只是可能，不一定。

import pycuda.autoinit
import pycuda.driver as drv
from pycuda import gpuarray
from pycuda.compiler import SourceModule
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import numpy as np