PyCUDA(3)-并行计算
本文将进一步介绍CUDA核函数在Python中的执行，以及线程索引在并行中的使用。
Array Add
让我们从一个实际的程序来入手。有两个随机数组A、B，他们拥有一致的长度，现在将他们按照元素对应位置相加得到数组C。在Python中可以这样：
import numpy
num = 4
A = numpy.random.rand(num)
B = numpy.random.rand(num)
C = A + B
print('A=', A)
print('B=', B)
print('C=', C)现在让我们在CUDA中以多种方式实现。
Array Add on GPU (easy版)
pyCUDA对于这种简单的功能甚至可以省略核函数，完全以python的方式实现。
import pycuda.autoinit
import pycuda.gpuarray as gpuarray
import numpy
num = 4
A = numpy.random.rand(num)
B = numpy.random.rand(num)
A_GPU = gpuarray.to_gpu(A.astype(numpy.float32))
B_GPU = gpuarray.to_gpu(B.astype(numpy.float32))
C_GPU = A_GPU + B_GPU
C = C_GPU.get()
print('A=', A)
print('B=', B)
print('C=', C)逐行讲解Python代码
在第2行我们导入了一个新的模块gpuarray，顾名思义GPU端的数组。
第7，8行通过模块内接口to_gpu上传了数组。astype是numpy模块用于转换数组格式的函数。
第9行，以如此简单自然的方式在GPU端将两个GPU端数组进行了相加。
第10行，通过get接口从GPU下载数据到CPU。
Array Add on GPU (normal版)
省略核函数在一些简单运算中是非常非常方便快捷的。但是对于一些复杂操作，我们仍然需要手动编写核函数，这也是本节内容的重点。如果你仔细阅读了之前的文章，那么这个略长的代码没有任何陌生的面孔。
import pycuda.autoinit
from pycuda.compiler import SourceModule
import pycuda.gpuarray as gpuarray
import numpy
mod = SourceModule(r"""
void __global__ add(const float *x, const float *y, float *z)
{
const int n = threadIdx.x;
z[n] = x[n] + y[n];
}
""")
add = mod.get_function("add")
num = 4
A = numpy.random.rand(num)
B = numpy.random.rand(num)
C = numpy.zeros(num)
A_GPU = gpuarray.to_gpu(A.astype(numpy.float32))
B_GPU = gpuarray.to_gpu(B.astype(numpy.float32))
C_GPU = gpuarray.to_gpu(B.astype(numpy.float32))
add(A_GPU, B_GPU, C_GPU, block=(num,1,1))
C = C_GPU.get()
print('A=', A)
print('B=', B)
print('C=', C)逐行讲解Python代码
为了使用核函数，我们增加了不少代码。现在回顾上一节的内容中两句话。
一个线程就是一个学生，一个block就是一间教室，一个grid就是一层教学楼，GPU就是这栋教学楼。
blockIdx(几年几班)与threadIdx（几排几列）一起共同确定了一个线程（学生）的位置。
首先我们观察核函数，在这里我们仅使用threadIdx.x这一个参数。
const int n = threadIdx.x;相应的我们在调用核函数时也仅给了参数block当中的第一个参数。且使其等于数组的大小。
add(A_GPU, B_GPU, C_GPU, block=(num,1,1))通过比喻来通俗化的描述这个数组相加的过程便是，我们将整个任务给了1年1班（grid默认等于1）。共需要4位同学来完成这个相加的任务，下面描述每个同学都做了什么。
0号同学(threadIdx.x=0)，将数组x，y位置为0的值相加后放入数组z的位置为0的地方。
1号同学(threadIdx.x=1)，将数组x，y位置为1的值相加后放入数组z的位置为1的地方。
2号同学(threadIdx.x=2)，将数组x，y位置为2的值相加后放入数组z的位置为2的地方。
3号同学(threadIdx.x=3)，将数组x，y位置为3的值相加后放入数组z的位置为3的地方。翻译成程序语言如下。
第n号线程将x[n]与y[n]相加后存入z[n]。下面我们考虑如果任务量比较大，一个班级又比较小，需要班级间并行。为了让不同班级的同学既不发生冲突（同时操作同一个位置的数据）又不漏掉某个位置没处理，我们需要安排好每个班级中的每个同学到哪个位置上去处理数据这一问题。希望下面的插图可以帮助你理解这个过程。
先给出修改后的代码。
import pycuda.autoinit
from pycuda.compiler import SourceModule
import pycuda.gpuarray as gpuarray
import numpy
mod = SourceModule(r"""
void __global__ add(const float *x, const float *y, float *z)
{
const int n = blockDim.x * blockIdx.x + threadIdx.x;
z[n] = x[n] + y[n];
}
""")
add = mod.get_function("add")
num = 6
A = numpy.random.rand(num)
B = numpy.random.rand(num)
C = numpy.zeros(num)
A_GPU = gpuarray.to_gpu(A.astype(numpy.float32))
B_GPU = gpuarray.to_gpu(B.astype(numpy.float32))
C_GPU = gpuarray.to_gpu(B.astype(numpy.float32))
add(A_GPU, B_GPU, C_GPU, grid=(2,), block=(4,1,1))
C = C_GPU.get()
print('A=', A)
print('B=', B)
print('C=', C)现将有改动的主要两句取出单独讲解（另外一处改动为将数组大小num改为6）。
const int n = blockDim.x * blockIdx.x + threadIdx.x;
add(A_GPU, B_GPU, C_GPU, grid=(2,), block=(4,1,1))我们在计算n的时候加入了blockDim.x * blockIdx.x。
blockDim.x 是一个所有同学看到都一样的值，既block在x方向上的数量（这里通俗翻译为总共有几个班级），等于参数grid的第一个值（这里等于2）。
blockIdx.x为当前线程（同学）在哪个block中（哪个班级）。
这里得益于C语言数组索引从0起始的便利性，使用下面的公式即可安排好每一个同学将要去处理数据的位置。
将要处理的位置 = 班级数量*当前班级号（从0开始）+ 当前同学号（从0开始）小结
CUDA核函数中最重要的概念线程索引，最大的意义便是确定每一个线程该去处理哪个一个位置上的数据。只要写核函数，这便是第一件要计算好的变量。而上面那个公式是最基础的计算方法，在更为复杂的问题中，会使用更复杂的计算公式。
Array Add on GPU (hard版)
如果你不想手动管理内存，那么下面的内容则可以略过。
首先让我们回顾一下如何在C语言中动态的分配两个随机数组并相加。
- 通过malloc动态分配内存。
- 初始化分配的数组。
- for循环遍历相加
- 释放动态分配的内存
在GPU中也是这样类似的过程，只是上面代码使用了gpuarray自动分配显存空间，如果你想手动管理显存，那么下面的代码就是为你准备的。
import pycuda.autoinit
import pycuda.driver as drv
from pycuda.compiler import SourceModule
import numpy
mod = SourceModule(r"""
void __global__ add(const float *x, const float *y, float *z)
{
const int n = blockDim.x * blockIdx.x + threadIdx.x;
z[n] = x[n] + y[n];
}
""")
add = mod.get_function("add")
num = 6
A = numpy.random.rand(num).astype(numpy.float32)
B = numpy.random.rand(num).astype(numpy.float32)
C = numpy.zeros(num).astype(numpy.float32)
A_GPU = drv.mem_alloc(A.nbytes)
B_GPU = drv.mem_alloc(B.nbytes)
C_GPU = drv.mem_alloc(C.nbytes)
drv.memcpy_htod(A_GPU, A)
drv.memcpy_htod(B_GPU, B)
add(A_GPU, B_GPU, C_GPU, grid=(2, 1), block=(4,1,1))
drv.memcpy_dtoh(C, C_GPU)
A_GPU.free()
B_GPU.free()
C_GPU.free()
print('A=', A)
print('B=', B)
print('C=', C)由于针对高级读者准备，所以不赘述说明，仅介绍几个要点。
numpy数组首先要做好格式转换，这样取numpy数组的nbytes才是准确的，或者可以手动计算空间长度。
memcpy_htod与CUDA C++语法类似，注意参数顺序虽然函数名为htod（主机到设备）但是参数顺序是（设备，主机）。memcpy_dtoh同理。
显存的释放是通过显存变量成员函数，而不是drv的函数，这点与C++不相似。
