矩阵乘法运算(一)-使用OpenMP加速循环计算
说到矩阵，只要是理工科都会想到被线性代数课支配的恐惧。矩阵乘法运算，往大了说各种工业和科研数值计算离不开它，往小了说各种跑分软件也会用到它，矩阵乘法运算的耗时也是评判计算机浮点运算性能的重要指标。本文的目的是通过矩阵乘法运算验证各种实现方式的性能差异，并对比不同计算平台的性能差异，为高性能计算开发提供参考。
1. 矩阵乘法算法
矩阵乘法运算也成为矩阵点乘，可以用下图表示1，通常由A矩阵对应行和B矩阵对应列的元素相乘并累加，形成C矩阵对应行列的值。要求A矩阵的列数要与B矩阵的行数相等，C矩阵的尺寸则相当于B矩阵的行数×A矩阵的列数。
用C语言表示一个M×N矩阵点乘一个N×T矩阵，一般写成3层循环嵌套的形式。
...
for (int i = 0; i < M; i++)
{
for (int j = 0; j < T; j++)
{
C[i,j] = 0.0;
for (int k = 0; k < N; k++)
{
C[i,j] += A[i,k] * B[k,j];
}
}
}
...
通常认为一个N×N维矩阵的乘法次数为N的3次方，用O(n3)表示计算复杂度。也有一些算法，将其中某些乘法运算转化成加法，减少了乘法运算次数。比如Strassen算法2，理论上计算复杂度只有O(n2.807)，降低了大矩阵相乘的耗时。更新的算法比如Coppersmith-Winograd方法，计算复杂度更是号称只有O(n2.3727)。这种算法上的差异不在本文讨论之列。
2. 循环嵌套并行化
现代计算机处理器一般都为多核心，为充分利用处理器性能，使用openmp等并行库对计算程序进行并行优化不失为一个好的选项。使用单核心进行循环嵌套求解没有实际意义，这里不做演示，有兴趣的可以在下面的代码中单独去掉并行库后进行编译计算。
2.1 C实现
C实现很简单，通过三层循环嵌套实现矩阵相乘算法。这里通过openmp将前两层循环合并分配到不同线程中展开计算3，以充分利用多核心处理器的性能。以下是openmp.c的示例代码。
#include <omp.h>
void matrix_multiply_float(int n, float A[], float B[], float C[])
{
#pragma omp parallel for collapse(2) shared(A, B, C)
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
C[i * n + j] = 0;
for (int k = 0; k < n; k++)
{
C[i * n + j] += A[i * n + k] * B[k * n + j];
}
}
}
}
void matrix_multiply_double(int n, double A[], double B[], double C[])
{
#pragma omp parallel for collapse(2) shared(A, B, C)
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
C[i * n + j] = 0;
for (int k = 0; k < n; k++)
{
C[i * n + j] += A[i * n + k] * B[k * n + j];
}
}
}
}
在main.c中定义一些命令行开关，用户可以自定义矩阵尺寸、循环次数、计算精度等。这里定义了一个性能参数，GFLOPs，用来衡量单位时间内执行浮点运算的能力，1GFLOPs相当于1秒内执行109次浮点运算。
