CUDA编程——矩阵乘法的串行和两种并行实现

一.CUDA是什么

这里仅简单介绍一下主要概念，如下：
1.主机
将CPU及系统的内存（内存条）称为主机。
2.设备
将GPU及GPU本身的显示内存称为设备。
3.线程(Thread)
一般通过GPU的一个核进行处理。
4.线程块(Block)
1. 由多个线程组成。
2. 各block是并行执行的，block间无法通信，也没有执行顺序。
3. 注意线程块的数量限制为不超过65535（硬件限制）。
5.线程格(Grid)
由多个线程块组成。

6.线程束
在CUDA架构中，线程束是指一个包含32个线程的集合，这个线程集合被“编织在一起”并且“步调一致”的形式执行。在程序中的每一行，线程束中的每个线程都将在不同数据上执行相同的命令。
7.核函数（Kernel）
1. 在GPU上执行的函数通常称为核函数。
2. 一般通过标识符__global__修饰，调用通过<<<参数1,参数2>>>，用于说明内核函数中的线程数量，以及线程是如何组织的。
3. 以线程格（Grid）的形式组织，每个线程格由若干个线程块（block）组成，而每个线程块又由若干个线程（thread）组成。
4. 是以block为单位执行的。
5. 只能在主机端代码中调用。
6. 调用时必须声明内核函数的执行参数。
7. 在编程时，必须先为kernel函数中用到的数组或变量分配好足够的空间，再调用kernel函数，否则在GPU计算时会发生错误。
8.dim3结构类型
1. dim3是基于uint3定义的矢量类型，相当于由3个unsigned int型组成的结构体。uint3类型有三个数据成员unsigned int x; unsigned int y; unsigned int z;
2. 可使用于一维、二维或三维的索引来标识线程，构成一维、二维或三维线程块。
3. dim3结构类型变量用在核函数调用的<<<,>>>中。
4. 相关的几个内置变量
threadIdx，获取线程thread的ID索引；
blockIdx，线程块的ID索引；
blockDim，线程块的维度；
gridDim，线程格的维度。
5. 对于一维的block，线程的threadID=threadIdx.x。
6. 对于大小为（blockDim.x, blockDim.y）的 二维 block，线程的threadID=threadIdx.x+threadIdx.y*blockDim.x。
7. 对于大小为（blockDim.x, blockDim.y, blockDim.z）的 三维 block，线程的threadID=threadIdx.x+threadIdx.y*blockDim.x+threadIdx.z*blockDim.x*blockDim.y。
8. 对于计算线程索引偏移增量为已启动线程的总数。如stride = blockDim.x * gridDim.x; threadId += stride。
9.函数修饰符
1. __global__，表明被修饰的函数在设备上执行，但在主机上调用。
2. __device__，表明被修饰的函数在设备上执行，但只能在其他__device__函数或者__global__函数中调用。
10.同步方法__syncthreads()
确保线程块中的每个线程都执行完__syscthreads()前面的语句后，才会执行下一条语句。
注意：
1. 当线程块的数量为GPU中处理数量的2倍时，将达到最优性能。
2. 核函数执行的第一个计算就是计算输入数据的偏移。每个线程的起始偏移都是0到线程数量减1之间的某个值。然后，对偏移的增量为已启动线程的总数。

更具体的原理可参考这篇博文CUDA编程入门：向量加法和矩阵乘法

二.矩阵乘法算法设计思路

A的大小为[10*blocksize][10*blocksize]， B的大小为[10*blocksize][20*blocksize]，求C矩阵的结果。

1.串行矩阵乘法

void searial(int *A, int *B, int *C)
{
    for (int i = 0; i < 10 * BLOCK_SIZE; i++)
    {
        for (int j = 0; j < 20 * BLOCK_SIZE; j++)
        {
            int sum = 0;

            for (int k = 0; k < 10 * BLOCK_SIZE; k++)
            {
                sum += A[i * 10 * BLOCK_SIZE + k] * 
                             B[k * 20 * BLOCK_SIZE + j];
            }

            C[i * 20 * BLOCK_SIZE + j] = sum;
        }
    }
}

2.并行不分块计算矩阵乘法，具体方法为每个thread计算C矩阵中的一个元素，计算得到row值和col值后即可直接for循环计算

__global__
void deviceParallel1(int *A, int *B, int *C)
{

    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    int sum = 0;

    for (int i =