程序性能优化指南——以矩阵乘法为例

© 2008–2018 by the MIT 6.172 Lecturers

衡量程序性能的指标

一般以FLOPS（Floating Point Operations Per Second，浮点运算每秒）作为衡量计算机性能的一种标准。

对于软件而言，FLOPS可以用以下公式表示：
$FLOPS= 浮点运算次数 / 运算时间（秒）$
如果一个程序需要执行 $1\times 10^9$ 次浮点运算，并且该程序在 2 秒内完成运算，则该程序的 FLOPS 为：

$FLOPS = \frac{ 1\times 10^9}{2} = 5\times 10^8 FLOPS = (500 \ MFLOPS)$

对于硬件而言，理论上的FLOPS的计算为：
$FLOPS（理论峰值）=每周期的浮点运算次数×时钟频率（GHz）×核心数$
假设一个处理器，每个周期可以执行 4 次浮点运算，时钟频率为 3.0 GHz，包含两个处理器，每个处理器有 8 个核心。
则理论峰值 FLOPS 为：
$FLOPS = 4 \times 3 \times 10^9 \times 8 \times 2 = 192 \times 10^9 \ FLOPS = 192 \ GFLOPS$

基于以上计算，要衡量程序性能可以用程序运算的实际FLOPS除以硬件峰值FLOPS, 该值越大，说明程序对硬件的利用率越高。

矩阵乘法的优化

一个典型的矩阵乘法，可以用以下过程表示：

图1

可以观察出，以上公式的运算复杂度为

O(n^3)

编程语言的影响

将图1中的矩阵运算写成最Python代码，如下：

'''
Example1-py
'''
import sys
import random
import time

n=512
A = [[random.random() for row in range(n)] for col in range(n)]
B = [[random.random() for row in range(n)] for col in range(n)]
C = [[0 for row in range(n)] for col in range(n)]

start = time.time()
for i in range(n):
    for j in range(n):
        for k in range(n):
            C[i][j] = A[i][k] * B[k][j]

end = time.time()
peak_FLOPS = 3.2 * 1e9 * 32 * 2 / 1e9
current_FLOPS = 2 * 512 * 512 * 512 * 64 /1e9 / (end-start)
utilization_rate = current_FLOPS / peak_FLOPS
print("Elapsed time is: {:.6f}s  ".format(end-start))
print("Current FLOPS is {} GFLOPS ".format(current_FLOPS))
print("Peak FLOPS is: {} GFLOPS".format(peak_FLOPS))
print("Utilization rate is: {:.6f}%".format(utilization_rate * 100) )

程序输出为：

Elapsed time is: 22.370289s  
Current FLOPS is 0.011999642055333309 GFLOPS 
Peak FLOPS is: 409.6 GFLOPS
Utilization rate is: 0.002930%

当矩阵维度为512时，运行时间为22.370289s `, 随着矩阵维度的增大，计算时间将会呈几何级数的增加，但CPU占用率只有 $0.002930\%$ 左右。

如果将以上代码改为C语言版本呢？

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define N 512

double random_double() {
    return (double)rand() / RAND_MAX;
}

int main() {
    double A[N][N], B[N][N], C[N][N] = {0};
    clock_t start, end;
    double peak_FLOPS, current_FLOPS, utilization_rate;

   
    srand(time(NULL));


    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            A[i][j] = random_double();
            B[i][j] = random_double();
        }
    }


    start = clock();


    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            for (int k = 0; k < N; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

 
    end = clock();

   
    double elapsed_time = ((double)(end - start)) / CLOCKS_PER_SEC;
    long long total_operations = 2LL * N * N * N;
    peak_FLOPS = 3.2e9 * 32 * 2 / 1e9; 
    current_FLOPS = (double)total_operations * 64 / 1e9 / elapsed_time; 
    utilization_rate = current_FLOPS / peak_FLOPS;  

  
    printf("Elapsed time is: %.6fs\n", elapsed_time);
    printf("Current FLOPS is: %.6f GFLOPS\n", current_FLOPS);
    printf("Peak FLOPS is: %.6f GFLOPS\n", peak_FLOPS);
    printf("Utilization rate is: %.6f%%\n", utilization_rate * 100);

    return 0;
}

程序运行结果为：

Elapsed time is: 0.538117s
Current FLOPS is: 0.498842 GFLOPS
Peak FLOPS is: 409.600000 GFLOPS
Utilization rate is: 0.121788%

可以看到，相同的硬件平台，C语言版本的矩阵乘法速度是0.560995s, 利用率是0.121788%, 编译型语言的性能相比于解释型的Python版本提高了一个数量级。
造成这种差异的结果，相信大家都知道，在此不做过多赘述。

方法	运行时间（s）	GFLOPS	CPU利用率
Python	21.726495	0.7907335	0.3861%
C	0.560995	30.623926	14.953089%

改变循环顺序的影响

假设我们更改循环顺序，对矩阵乘法会有多大的影响呢？
如下，我们将i,j,k 的遍历顺序，更改为i,k,j，如下代码所示：

for (int i = 0; i < N; i++) {
        for (int k = 0; k < N; k++) {
            for (int j = 0; j < N; j++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

程序运行结果为：

Elapsed time is: 0.262266s
Current FLOPS is: 1.023524 GFLOPS
Peak FLOPS is: 409.600000 GFLOPS
Utilization rate is: 0.249884%

不断改变遍历的顺序，程序运行结果如下表

遍历顺序	运行时间（s）	GFLOPS	CPU利用率
i,j,k	0.560995	0.473690	0.115647%
i,k,j	0.251579s	1.023524	0.249884%
j,k,i	0.876779s	0.283216	0.069145%
k,i,j	0.287382s	1.023914	0.249979%
k,j,i	0.873372s	0.294188	0.071823%

可以看到，仅仅是改变循环的顺序，程序性能就提升了一倍以上。实际上，造成这一结果的原因就是缓存，即也就是数据访问中常常提到的缓存命中。

当处理器需要访问某个数据时，首先会在缓存中查找。如果数据已经在缓存中存在，这就称为缓存命中，相比于从RAM中读取数据，从缓存中读取数据可以极大减少数据的访问延迟，进而提升程序性能。

而对于多维度矩阵而言，数据的排列方式，会影响缓存中的暂时存储的数据分布，进而影响程序性能。即也就是，我们常说的行优先（Row-major order）和列优先（Column-major order）。

行优先和列优先

行优先即指在多维数组的存储方式中，内存中的数据按行顺序存储。

Matrix

行优先的数据内存分布

C/C++一般都是行优先分布，而Fortan一般采用列优先存储二维数组

`i,j,k`的数据访问模式

以行优先为例，当以i,j,k的方式遍历二维数组时，矩阵乘法每次计算将会读取矩阵A的一整行，和矩阵B的一个元素，此时，最内层循环的数据访问是不连续的，所以，程序缓存命中的效率大大降低，进而影响程序性能。

i,j,k访问

`i,k,j`的数据访问模式

当以i,k,j的方式遍历二维数组时，矩阵乘法每次计算将会读取矩阵A的一个元素，和矩阵B的，此时，最内层循环的数据访问是连续的，计算机会将可以将矩阵B 中连续的数据存放到缓存中，缓存命中的概率大大提升，降低了数据访问的延迟，进而提高了程序性能。

i,k,j访问

`j,k,i`的数据访问模式

当程序以j,k,i的方式进行访问时，可以看到此时，A,B,C三个矩阵都不是连续的内存访问，此时，程序性能最差。

j,k,i

查看程序的缓存命中率

在Linux系统上，可以通过valgrind 查看末级缓存未命中率运行如下命令:

末级缓存未命中率, 指在处理器的最后一级缓存（通常是 L3 缓存）中，发生缓存未命中的次数占总访问次数的比例)。

valgrind --tool=cachegrind ./SquareMatrix

不同循环顺序的缓存命中率如下：

遍历顺序	运行时间（s）	GFLOPS	CPU利用率	末级缓存未命中率
i,j,k	0.560995	0.473690	0.115647%	7.7%
i,k,j	0.251579s	1.023524	0.249884%	1.0%
j,k,i	0.876779s	0.283216	0.069145%	15.4%
k,i,j	0.287382s	1.023914	0.249979%	1.0%
k,j,i	0.873372s	0.294188	0.071823%	15.4%

总结

数据的访问方式影响程序性能，我们应该尽可能利用计算机的缓存命中原理，尽可能访问连续的内存数据，提高缓存命中率，减少数据访问的延迟，进而提升程序性能，就以二维矩阵乘法而言，应该尽可能保证最内层的循环访问连续的数据，提高缓存命中率，以提高程序性能。

循环并行（Parallel-Loop）

对于循环问题，实际上，可以通过并行处理加速实现，为了验证并行优化的性能，我们不妨增加矩阵的维度，将矩阵的维度增加至2048. 但是，使用并行计算必须注意两个问题：

在堆上初始化，当矩阵维度增大时，可能会超过栈的默认大小，因此使用堆内存可以解决运行过程中可能带来的Segmentation Fault，堆内存的初始化如下：

float **A = malloc(N * sizeof(float *));
float **B = malloc(N * sizeof(float *));
float **C = malloc(N * sizeof(float *));
for (int i = 0; i < N; i++) {
    A[i] = malloc(N * sizeof(float));
    B[i] = malloc(N * sizeof(float));
    C[i] = malloc(N * sizeof(float));
}

并行运算的时间计算，基于OpenMP的编程， clock_t start 和 clock_t end 测量的是 CPU 时间。在多线程情况下，使用 omp_get_wtime() 计算真实的会更加准确。

综上具体代码如下：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <omp.h>

#define N 2048

float random_double() {
    return (float)rand() / RAND_MAX;
}

int main() {
    //float A[N][N], B[N][N], C[N][N] = {0};
   // clock_t start, end;
    double peak_FLOPS, current_FLOPS, utilization_rate;

   
    srand(time(NULL));

    
float **A = malloc(N * sizeof(float *));
float **B = malloc(N * sizeof(float *));
float **C = malloc(N * sizeof(float *));
for (int i = 0; i < N; i++) {
    A[i] = malloc(N * sizeof(float));
    B[i] = malloc(N * sizeof(float));
    C[i] = malloc(N * sizeof(float));
}


    //clock_t start = clock();
    double start = omp_get_wtime();
    #pragma omp parallel for schedule(dynamic)
    for (int i = 0; i < N; i++) {
        
        for (int k = 0; k < N; k++) {
            
            for (int j = 0; j < N; j++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

 
    //clock_t end = clock();
    double end = omp_get_wtime();

    double elapsed_time = end - start;
    //double elapsed_time = ((double)(end - start)) / CLOCKS_PER_SEC;
    long long total_operations = 2LL * N * N * N; 
    //末级缓存未命中率
    peak_FLOPS = 16 * 3.2 * 8;  // 计算峰值 FLOPS
    current_FLOPS = (double)total_operations  / 1e9 / elapsed_time;  // 当前 FLOPS
    utilization_rate = current_FLOPS / peak_FLOPS;  // 利用率

  
    printf("Elapsed time is: %.6fs\n", elapsed_time);
    printf("Current FLOPS is: %.6f GFLOPS\n", current_FLOPS);
    printf("Peak FLOPS is: %.6f GFLOPS\n", peak_FLOPS);
    printf("Utilization rate is: %.6f%%\n", utilization_rate * 100);

    return 0;
}

执行以下编译指令

gcc -fopenmp -o SquareMatrix SquareMatrix.c

程序结果如下：

Elapsed time is: 4.058899s
Current FLOPS is: 4.232643 GFLOPS
Peak FLOPS is: 409.600000 GFLOPS
Utilization rate is: 1.033360%

对比未使用并行循环的结果

Elapsed time is: 27.237198s
Current FLOPS is: 0.630750 GFLOPS
Peak FLOPS is: 409.600000 GFLOPS
Utilization rate is: 0.153992%

对比并行前和并行后的程序，可以发现程序运行速度和峰值FLOPS大约提升了6倍。

对比并行循环的位置

如果我们调整并行循环的位置，程序性能会发生何种变化呢？

调整至第二层循环

for (int i = 0; i < N; i++) {
        #pragma omp parallel for schedule(dynamic)
        for (int k = 0; k < N; k++) {
            for (int j = 0; j < N; j++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

程序结果如下：

Elapsed time is: 12.778550s
Current FLOPS is: 1.344430 GFLOPS
Peak FLOPS is: 409.600000 GFLOPS
Utilization rate is: 0.328230%

调整至最内层循环

for (int i = 0; i < N; i++) {
  for (int k = 0; k < N; k++) {
            #pragma omp parallel for schedule(dynamic)
            for (int j = 0; j < N; j++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

Elapsed time is: 726.218460s
Current FLOPS is: 0.023657 GFLOPS
Peak FLOPS is: 409.600000 GFLOPS
Utilization rate is: 0.005776%

Amazing!，我们发现，最内层并行，运行时间达到了726s, 并行运算完全起到了副作用，这是因为，并行循环在最外层时，程序性能更高，具体原因在于：外层循环的并行化可以减少线程启动和调度的频率。因为内层循环被并行化时，每次进入外层循环时都会启动新的线程和分配任务，而外层循环被并行化后，线程的创建和分配只发生一次，减少了开销。