本实验采用不同的方法来计算 8192 * 8192 的整型矩阵乘法运算。

C语言版

C语言是大家公认的高性能语言，那我们就从C语言开始吧。

// 用一位数组表示二维矩阵 mat1 = (int*) malloc(m_size * m_size * sizeof(int));

mat2 = (int*) malloc(m_size * m_size * sizeof(int));

result = (int*) malloc(m_size * m_size * sizeof(int));

// initialize for (int i = 0; i < m_size * m_size; i++) {

mat1[i] = rand()/1000000;

mat2[i] = rand()/1000000;

result[i] = 0;

}

for (int r = 0; r < m_size; r++) {

for (int c = 0; c < m_size; c++) {

for (int n = 0; n < m_size; n++) {

result[r*m_size + c] += mat1[r*m_size+n] * mat2[n*m_size+c];

}

}

}

这代码没什么可说的，值得一提的可能就是用一维数组来存二维矩阵，这样可以让矩阵在内存中的分布更连续。很容易估算出来这段代码需要进行万亿级别的乘法或加法运算，倘若使用1950年冯.诺依曼用来预测天气的计算机(每秒5000次计算)，大概要算个六七年吧。

还好现在是2019年，这段代码在我的 i7 cpu 上运行花了 7000 s，差不多两个小时了。(记得开编译器优化，不然可能今天都跑不完)

eigen 版

#include using namespace::Eigen;

typedef Eigen::Matrix IntMatrix;

IntMatrix mat1(m_size,m_size);

IntMatrix mat2(m_size,m_size);

IntMatrix result(m_size,m_size);

mat1.setRandom();

mat2.setRandom();

result = mat1 * mat2;

运行此段代码需要下载 eigen 库，并在编译时link上去。这段代码在我的机器(i7 单核)上运行花了 120s，从两个小时直接变成了两分钟，代码还简单了很多啊握草。(同样记得开编译器优化)

给 eigen 跪了啊，这比我快了快60倍啊，这告诉我们一个道理，轮子该用还是得用啊。( eigen 进行了指令级级别的优化，以后有机会在写篇文章深入讲下。)

CUDA

显然，矩阵乘法是一个天然完美的可并行问题。事实上，大量伟大的理论研究和工程实践都基于矩阵乘法可并行这一前提，比如google。

CPU并不擅长并行计算，但是GPU擅长啊。GPU可以一口气拉起一个网格(grid)的线程，像一个计算方针一样同时处理计算的不同部分。关于CUDA可以参考我的这篇文章陈敏华：cuda实战入门​zhuanlan.zhihu.com

在矩阵乘法这个问题中，结果矩阵中的每一个位置的数据都不依赖于其他位置数据的计算过程和结果。所以我们完全可以拉一堆进程起来，让他们分别计算某一些位置的结果。

你可以想象为你们学校要算两个100 * 100的矩阵的乘法，这个矩阵计算的结果是万级的，而运算量可是百万级的。于是你们校长找到50个班主任，给每个班主任分配两行结果计算任务，你们班呢，分到了前两行，也就是要算出结果矩阵中前两行的200个数字。班主任再把这两百个数字分给了50个同学，每个同学只要算4个数字，大约400次加法和400次乘法。于是一个巨复杂的计算问题，你们学校半天就完成了。

下面是我的代码实现

#include #include#include

#define BLOCK_NUM 32//块数量#define THREAD_NUM 256// 每个块中的线程数#define R_SIZE BLOCK_NUM * THREAD_NUM#define M_SIZE R_SIZE * R_SIZE

__global__ void mat_mul(int *mat1, int *mat2, int *result) {

const int bid = blockIdx.x;

const int tid = threadIdx.x;

// 每个线程计算一行 const int row = bid * THREAD_NUM + tid;

for (int c = 0; c < R_SIZE; c++) {

for (int n = 0; n < R_SIZE; n++) {

result[row*R_SIZE+c] += mat1[row*R_SIZE+n] * mat2[n*R_SIZE+c];

}

}

}

int main(int argc, char *argv[]) {

int *mat1, *mat2, *result;

int *g_mat1, *g_mat2, *g_mat_result;

// 用一位数组表示二维矩阵 mat1 = (int*) malloc(M_SIZE * sizeof(int));

mat2 = (int*) malloc(M_SIZE * sizeof(int));

result = (int*) malloc(M_SIZE * sizeof(int));

// initialize for (int i = 0; i < M_SIZE; i++) {

mat1[i] = rand()/1000000;

mat2[i] = rand()/1000000;

result[i] = 0;

}

cudaMalloc((void **)&g_mat1, sizeof(int) * M_SIZE);

cudaMalloc((void **)&g_mat2, sizeof(int) * M_SIZE);

cudaMalloc((void **)&g_mat_result, sizeof(int) * M_SIZE);

cudaMemcpy(g_mat1, mat1, sizeof(int) * M_SIZE, cudaMemcpyHostToDevice);

cudaMemcpy(g_mat2, mat2, sizeof(int) * M_SIZE, cudaMemcpyHostToDevice);

mat_mul<<>>(g_mat1, g_mat2, g_mat_result);

cudaMemcpy(result, g_mat_result, sizeof(int) * M_SIZE, cudaMemcpyDeviceToHost);

}

这段代码在我的机器(RTX 2080TI)上运行了多久呢？

22秒！！！

差不多比c写的第一版要快了300倍了。

Tensorflow

mat1 = tf.random_uniform([8192, 8192], minval=1, maxval=65536, dtype=tf.int32)

mat2 = tf.random_uniform([8192, 8192], minval=1, maxval=65536, dtype=tf.int32)

sess = tf.Session()

sess.run(mat1)

sess.run(mat2)

sess.run(tf.matmul(mat1, mat2))

sess.close()

这段代码在我的机器上也运行了 22秒，和上面我自己手写cuda程序性能相当。

CUBLAS

22秒已经到极限了吗？还早得很呢！我们可以使用CUBLAS库，不但可以屏蔽掉复杂的底层实现以及不同计算设备带来的参数设计的影响，还有机会把矩阵乘法的效率进一步提升。先贴出cublas的一份文档 。cublas doc​developer.nvidia.com

下面是我的代码实现

// 用一位数组表示二维矩阵 mat1 = (float*) malloc(m_size * sizeof(float));

mat2 = (float*) malloc(m_size * sizeof(float));

result = (float*) malloc(m_size * sizeof(float));

// initialize for (int i = 0; i < m_size; i++) {

mat1[i] = rand()/10000000;

mat2[i] = rand()/10000000;

result[i] = 0;

}

cudaMalloc((void **)&g_mat1, sizeof(*mat1) * m_size);

cudaMalloc((void **)&g_mat2, sizeof(*mat2) * m_size);

cudaMalloc((void **)&g_mat_result, sizeof(*result) * m_size);

// initialize CUBLAS context stat = cublasCreate(&handle);

stat = cublasSetMatrix(r_size, r_size, sizeof(*mat1), mat1, r_size, g_mat1, r_size);

stat = cublasSetMatrix(r_size, r_size, sizeof(*mat2), mat2, r_size, g_mat2, r_size);

stat = cublasSetMatrix(r_size, r_size, sizeof(*result), result, r_size, g_mat_result, r_size);

float al = 1.0f;

float bet = 0.0f;

stat = cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,

r_size, r_size, r_size, &al, g_mat1,

r_size, g_mat2, r_size, &bet, g_mat_result, r_size);

stat = cublasGetMatrix(r_size, r_size, sizeof(*result), g_mat_result, r_size, result, r_size);

其中的 cublasSgemm方法就是计算矩阵乘法的函数，具体请参考上面pdf的2.4小节。

这段代码计算8192 * 8192的矩阵乘法要多久呢？

答案是：1.35秒 ！！！我们在第一版c语言实现的基础上完成了5000倍的提升！！！

总结

本篇文章中，我们将 8192 * 8192 的矩阵运算从 7000秒一路提升到了1.35秒，提升了五千多倍。这是一个成就感与挫败感并存的过程，一方面不断提升计算效率很爽，一方面又感觉这种提升并非来自我的智慧而是他人完成的工作。

这篇文章虽然暂时告一段落，但是我相信，这里面依然有各种优化的空间。另外，如果矩阵更大呢？如果是稀疏矩阵呢？如果想利用多块GPU呢？

游戏才刚刚开始呢。

PS: 本文代码我都放到 github上啦。chenminhua/math​github.com