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 = 0; i < 10 * BLOCK_SIZE; i++)    {        sum += A[row * 10 * BLOCK_SIZE + i] * B[i * 20 * BLOCK_SIZE + col];    }    C[row * 20 * BLOCK_SIZE + col] = sum;}void parallel1(int *A, int *B, int *C){    int *CA, *CB, *CC;    cudaMalloc(&CA, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE);    cudaMalloc(&CB, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);     cudaMalloc(&CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);    cudaMemcpy(CA, A, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE,                                         cudaMemcpyHostToDevice);    cudaMemcpy(CB, B, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE,                                         cudaMemcpyHostToDevice);    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);     dim3 dimGrid(20, 10);    deviceParallel1<<<dimBlock, dimGrid>>>(CA, CB, CC);    cudaThreadSynchronize();    cudaMemcpy(C, CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE,                                         cudaMemcpyDeviceToHost);    cudaFree(CA);    cudaFree(CB);    cudaFree(CC);}

3.并行分块计算矩阵乘法，具体方法为根据线程块号和块内线程号，循环遍历得到所有的子矩阵，在主要循环计算步骤前后添加同步线程语句，最后将结果写回到C矩阵

__global__void deviceParallel2(int *A, int *B, int *C){    //获得线程块号    int blkRow = blockIdx.y;     int blkCol = blockIdx.x;    //获得块内的线程号     int row = threadIdx.y;     int col = threadIdx.x;    int var = 0;    //循环，遍历所有子矩阵    for (int i = 0; i < 10; i++)     {           const int *ASub = A + blkRow * BLOCK_SIZE * 10 + i * BLOCK_SIZE;         const int *BSub = B + i * BLOCK_SIZE * 20 + blkCol * BLOCK_SIZE;        __shared__ int Ads[BLOCK_SIZE][BLOCK_SIZE];         __shared__ int Bds[BLOCK_SIZE][BLOCK_SIZE];        Ads[row][col] = *(ASub + row * BLOCK_SIZE * 10 + col);         Bds[row][col] = *(BSub + row * BLOCK_SIZE * 20 + col);        __syncthreads();        for (int i = 0; i < BLOCK_SIZE; i++)         {            var += Ads[row][i] * Bds[i][col];         }        __syncthreads();    }    int *CSub = C + blkRow * BLOCK_SIZE * 20 + blkCol * BLOCK_SIZE;    *(CSub + row * BLOCK_SIZE * 20 + col) = var;}void parallel2(int *A, int *B, int *C){    int *CA, *CB, *CC;    cudaMalloc(&CA, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE);    cudaMalloc(&CB, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);     cudaMalloc(&CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);    cudaMemcpy(CA, A, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE,                                         cudaMemcpyHostToDevice);    cudaMemcpy(CB, B, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE,                                         cudaMemcpyHostToDevice);    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);     dim3 dimGrid(20, 10);    deviceParallel2<<<dimBlock, dimGrid>>>(CA, CB, CC);    cudaThreadSynchronize();    cudaMemcpy(C, CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE,                                         cudaMemcpyDeviceToHost);    cudaFree(CA);    cudaFree(CB);    cudaFree(CC);}

三.运行结果及加速比计算

计算加速比：(blocksize设置为32，若设置为16,则计算结果在时间上比较不明显，而64的话有可能会导致内存溢出)
第一种不分块并行算法：speedup = 656361/412700 = 1.59
第二种分块并行算法：speedup = 656361/5464 = 120.12

四.源代码

#include <stdio.h> #include <stdlib.h>#include <time.h>#include <sys/time.h>#include <iostream>#include <cuda_runtime.h>using namespace std;#define BLOCK_SIZE 32typedef void (*multiply)(int *A, int *B, int *C);double getTime(int *A, int *B, int *C, multiply mul){    timeval start, finish;    gettimeofday(&start, 0);    mul(A, B, C);    gettimeofday(&finish, 0);    double interval = 1e6 * (finish.tv_sec - start.tv_sec) + finish.tv_usec - start.tv_usec;    return interval;}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;        }    }}__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 = 0; i < 10 * BLOCK_SIZE; i++)    {        sum += A[row * 10 * BLOCK_SIZE + i] * B[i * 20 * BLOCK_SIZE + col];    }    C[row * 20 * BLOCK_SIZE + col] = sum;}void parallel1(int *A, int *B, int *C){    int *CA, *CB, *CC;    cudaMalloc(&CA, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE);    cudaMalloc(&CB, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);     cudaMalloc(&CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);    cudaMemcpy(CA, A, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE,                                         cudaMemcpyHostToDevice);    cudaMemcpy(CB, B, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE,                                         cudaMemcpyHostToDevice);    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);     dim3 dimGrid(20, 10);    deviceParallel1<<<dimBlock, dimGrid>>>(CA, CB, CC);    cudaThreadSynchronize();    cudaMemcpy(C, CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE,                                         cudaMemcpyDeviceToHost);    cudaFree(CA);    cudaFree(CB);    cudaFree(CC);}__global__void deviceParallel2(int *A, int *B, int *C){    //获得线程块号    int blkRow = blockIdx.y;     int blkCol = blockIdx.x;    //获得块内的线程号     int row = threadIdx.y;     int col = threadIdx.x;    int var = 0;    //循环，遍历所有子矩阵    for (int i = 0; i < 10; i++)     {           const int *ASub = A + blkRow * BLOCK_SIZE * 10 + i * BLOCK_SIZE;         const int *BSub = B + i * BLOCK_SIZE * 20 + blkCol * BLOCK_SIZE;        __shared__ int Ads[BLOCK_SIZE][BLOCK_SIZE];         __shared__ int Bds[BLOCK_SIZE][BLOCK_SIZE];        Ads[row][col] = *(ASub + row * BLOCK_SIZE * 10 + col);         Bds[row][col] = *(BSub + row * BLOCK_SIZE * 20 + col);        __syncthreads();        for (int i = 0; i < BLOCK_SIZE; i++)         {            var += Ads[row][i] * Bds[i][col];         }        __syncthreads();    }    int *CSub = C + blkRow * BLOCK_SIZE * 20 + blkCol * BLOCK_SIZE;    *(CSub + row * BLOCK_SIZE * 20 + col) = var;}void parallel2(int *A, int *B, int *C){    int *CA, *CB, *CC;    cudaMalloc(&CA, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE);    cudaMalloc(&CB, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);     cudaMalloc(&CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE);    cudaMemcpy(CA, A, sizeof(int) * 10 * BLOCK_SIZE * 10 * BLOCK_SIZE,                                         cudaMemcpyHostToDevice);    cudaMemcpy(CB, B, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE,                                         cudaMemcpyHostToDevice);    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);     dim3 dimGrid(20, 10);    deviceParallel2<<<dimBlock, dimGrid>>>(CA, CB, CC);    cudaThreadSynchronize();    cudaMemcpy(C, CC, sizeof(int) * 10 * BLOCK_SIZE * 20 * BLOCK_SIZE,                                         cudaMemcpyDeviceToHost);    cudaFree(CA);    cudaFree(CB);    cudaFree(CC);}void read(int *M, int row, int col){    srand((unsigned)time(NULL));    for (int i = 0; i < row; i++)    {        for (int j = 0; j < col; j++)        {            M[i * col + j] = rand() % 1000;        }    }}int main(int argc, char const *argv[]){    int *A = new int[10 * BLOCK_SIZE * 10 * BLOCK_SIZE];    if (A == NULL)    {        cerr << "can not malloc A" << endl; exit(1);    }    int *B = new int[10 * BLOCK_SIZE * 20 * BLOCK_SIZE];    if (B == NULL)    {        cerr << "can not malloc B" << endl; exit(1);    }    int *C1 = new int[10 * BLOCK_SIZE * 20 * BLOCK_SIZE];    if (C1 == NULL)    {        cerr << "can not malloc C" << endl; exit(1);    }    int *C2 = new int[10 * BLOCK_SIZE * 20 * BLOCK_SIZE];    if (C2 == NULL)    {        cerr << "can not malloc C" << endl; exit(1);    }    int *C3 = new int[10 * BLOCK_SIZE * 20 * BLOCK_SIZE];    if (C3 == NULL)    {        cerr << "can not malloc C" << endl; exit(1);    }    // 读取矩阵数据    read(A, 10 * BLOCK_SIZE, 10 * BLOCK_SIZE);    read(B, 10 * BLOCK_SIZE, 20 * BLOCK_SIZE);    cout << "Serial Time = " << getTime(A, B, C1, searial) << " ps." << endl;    cout << "Parallel1 Time = " << getTime(A, B, C2, parallel1) << " ps." << endl;    cout << "Parallel2 Time = " << getTime(A, B, C3, parallel2) << " ps." << endl;    delete[] A;    delete[] B;    delete[] C1;    delete[] C2;    delete[] C3;    return 0;}