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 32
typedef 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;
}