鲲鹏平台HPC优化：矩阵乘法从基础到专业数学库的演进-嵌云网-嵌入式AI开发资源站

鲲鹏平台HPC优化：矩阵乘法从基础到专业数学库的演进

青清朝阳

1. 鲲鹏平台HPC优化实战：从朴素实现到专业数学库的矩阵乘法演进

在HPC（高性能计算）领域，矩阵乘法是最基础也是最核心的运算之一。作为鲲鹏平台的开发者，我最近完成了一个从最基础实现到高度优化的完整优化过程，将1024×1024双精度矩阵乘法的性能从最初的0.26 GFLOPS提升到了142.35 GFLOPS，达到了鲲鹏920处理器理论峰值的85.5%。这个过程中积累的经验和教训，值得与各位同行分享。

1.1 为什么选择矩阵乘法作为优化案例

矩阵乘法之所以成为HPC优化的经典案例，主要基于三个重要特性：

首先，矩阵运算在科学计算和工程应用中无处不在。从深度学习训练中的卷积运算到流体力学仿真，矩阵乘法往往占据了60%-80%的计算时间。优化矩阵乘法可以直接提升这些应用的性能。

其次，矩阵乘法具有极大的优化空间。最简单的三层循环实现可能只能达到理论性能的0.1%，而经过系统优化后可以达到80%以上。这中间的数百倍差距，正好可以用来展示各种优化技术的效果。

最后，矩阵乘法的性能可以精确量化。通过GFLOPS（每秒十亿次浮点运算）这个指标，我们可以清楚地看到每个优化步骤带来的性能提升，便于分析和比较。

2. 矩阵乘法优化的五个阶段

2.1 版本1：朴素实现（基准版本）

朴素实现是理解算法本质的最佳起点。这个版本使用最直接的三层循环结构：

c复制void matrix_multiply(double *A, double *B, double *C, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            double sum = 0.0;
            for (int k = 0; k < n; k++) {
                sum += A[i * n + k] * B[k * n + j];
            }
            C[i * n + j] = sum;
        }
    }
}

在鲲鹏920上测试1024×1024矩阵，这个版本仅获得0.26 GFLOPS的性能。性能分析显示，主要瓶颈在于内存访问模式：

对矩阵B的访问是按列进行的（k循环在最内层），这导致严重的缓存未命中
每次内层循环都需要从内存加载新的数据，无法利用缓存局部性
编译器优化空间有限，因为内存访问模式限制了指令级并行

2.2 版本2：缓存优化（循环重排）

通过简单的循环顺序调整，我们实现了第一个重要优化：

c复制void matrix_multiply(double *A, double *B, double *C, int n) {
    for (int i = 0; i < n; i++) {
        for (int k = 0; k < n; k++) {
            double a_ik = A[i * n + k];
            for (int j = 0; j < n; j++) {
                C[i * n + j] += a_ik * B[k * n + j];
            }
        }
    }
}

关键改变是将k循环提到中间层，这使得：

对B矩阵的访问变为按行进行，大幅提高了空间局部性
a_ik可以在j循环中重复使用，减少了内存访问次数
C矩阵的写入模式更加连续

这个简单的调整带来了4.7倍的性能提升，达到1.22 GFLOPS。缓存命中率从15%提升到72%，验证了内存访问模式对性能的关键影响。

2.3 版本3：OpenMP并行化

鲲鹏920拥有32个物理核心，利用OpenMP可以充分发挥多核优势：

c复制void matrix_multiply(double *A, double *B, double *C, int n) {
    #pragma omp parallel for
    for (int i = 0; i < n; i++) {
        for (int k = 0; k < n; k++) {
            double a_ik = A[i * n + k];
            for (int j = 0; j < n; j++) {
                C[i * n + j] += a_ik * B[k * n + j];
            }
        }
    }
}

编译时需要添加-fopenmp选项，并设置合适的线程数：

bash复制gcc -O3 -fopenmp openmp_parallel.c -o openmp_parallel
export OMP_NUM_THREADS=32
./openmp_parallel

这个版本达到了29.38 GFLOPS，相比串行版本提升了24.1倍，并行效率为75.2%。性能损失主要来自：

线程创建和同步开销
不同行的计算量可能不均衡
多个线程同时写入相邻内存导致的伪共享（false sharing）

2.4 版本4：分块优化 + OpenMP

为进一步提高缓存利用率，我们引入了分块技术：

c复制#define BLOCK_SIZE 64

void matrix_multiply(double *A, double *B, double *C, int n) {
    #pragma omp parallel for collapse(2)
    for (int ii = 0; ii < n; ii += BLOCK_SIZE) {
        for (int jj = 0; jj < n; jj += BLOCK_SIZE) {
            for (int kk = 0; kk < n; kk += BLOCK_SIZE) {
                // 分块内计算
                for (int i = ii; i < ii + BLOCK_SIZE && i < n; i++) {
                    for (int k = kk; k < kk + BLOCK_SIZE && k < n; k++) {
                        double a_ik = A[i * n + k];
                        for (int j = jj; j < jj + BLOCK_SIZE && j < n; j++) {
                            C[i * n + j] += a_ik * B[k * n + j];
                        }
                    }
                }
            }
        }
    }
}

选择64×64的分块大小是为了适配鲲鹏920的64KB L1缓存。每个分块需要：

输入块A：64×64×8B = 32KB
输入块B：64×64×8B = 32KB
输出块C：64×64×8B = 32KB

这个版本性能提升到44.69 GFLOPS，L1缓存命中率达到94.3%。关键优化点包括：

更好的数据局部性，块内数据可以保留在缓存中
减少了NUMA架构下的远程内存访问
使用collapse(2)将两层循环并行化，改善负载均衡

2.5 版本5：鲲鹏数学库优化

最终，我们使用鲲鹏数学库(KML)中的BLAS实现：

c复制#include "kblas.h"

void matrix_multiply(double *A, double *B, double *C, int n) {
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                n, n, n, 1.0, A, n, B, n, 0.0, C, n);
}

这个高度优化的实现包含了：

手工调优的汇编内核
自动适配的最佳分块策略
NEON/SVE向量化指令
多级并行化（线程+指令级）

性能达到惊人的142.35 GFLOPS，是初始版本的548.9倍。这个案例充分展示了专业数学库的价值。

3. 性能分析与优化技术总结

3.1 各版本性能对比

优化版本	GFLOPS	相对提升	关键技术
朴素实现	0.26	1x	基础实现
缓存优化	1.22	4.7x	循环重排
OpenMP并行	29.38	24.1x	多线程
分块+OpenMP	44.69	1.5x	缓存分块
KML BLAS	142.35	3.2x	专业优化

3.2 优化技术协同效应

最关键的发现是各种优化技术之间存在协同效应：

缓存优化为多线程并行奠定了基础，减少了伪共享问题
分块技术进一步降低了NUMA架构下的远程内存访问
数学库集成了所有优化技术，并加入了SIMD向量化

这种技术叠加带来的不是简单的加法效应，而是乘法效应。这也解释了为什么专业数学库的性能如此出色。

4. 实战经验与常见问题

4.1 环境配置建议

在鲲鹏平台上获得最佳性能需要注意：

编译器选项：

bash复制gcc -O3 -march=armv8.2-a -fopenmp

线程绑定设置：

bash复制export OMP_PROC_BIND=close
export OMP_PLACES=cores

NUMA控制：

bash复制numactl --cpunodebind=0,1 --membind=0,1 ./program

4.2 常见问题解决

问题1：BLAS库找不到

解决方案：

bash复制sudo yum install openblas-devel
# 或者指定路径
gcc -I/usr/include/openblas -lopenblas

问题2：并行效率低

检查点：

使用perf stat分析缓存命中率
检查线程是否均匀分布在NUMA节点上
验证是否有其他进程占用CPU资源

问题3：分块大小选择

通过实验确定最佳分块大小：

c复制// 测试不同分块大小
for (int bs = 32; bs <= 128; bs *= 2) {
    // 测试性能
}

5. 进阶优化方向

对于追求极致性能的开发者，还可以考虑：

手动编写ARM汇编内核，充分利用SVE指令集
混合精度计算，在适当场景使用FP16
异步计算与数据预取重叠
针对特定矩阵形状（如稀疏矩阵）的特化优化

在鲲鹏平台上进行HPC优化是一段充满挑战和成就感的旅程。从最简单的实现开始，逐步应用各种优化技术，最终达到接近理论峰值的性能，这个过程不仅提升了代码效率，也深化了对计算机体系结构的理解。