登录/ 注册
 找回密码
 立即注册
关闭

CSDN热搜

程序园 精品问答 技术交流 资源下载
返回列表 发新帖
首页 › 业界区 › 安全 ›[C#] 使用 .NET 的跨平台SIMD硬件加速功能，将 GEMM（通 ...

[C#] 使用 .NET 的跨平台SIMD硬件加速功能，将 GEMM（通用矩阵乘法）

陆菊 4 天前
目录

  • 一、基本算法与测试办法

    • 1.1 矩阵乘法的定义

      • 1.1.1 矩阵形状与运算复杂度

    • 1.2 C++ 开发的矩阵乘法最基础实现
    • 1.3 C# 开发的矩阵乘法最基础实现（Basic）

      • 1.3.1 矩阵乘法的实现
      • 1.3.2 基准测试方法

    • 1.4 增加MathNet、MKL、OpenBLAS的基准测试

      • 1.4.1 引入库
      • 1.4.2 这些库的基准测试方法

        • 1.4.2.1 MathNet的基准测试方法（UseMathNet）
        • 1.4.2.2 OpenBLAS的基准测试方法（UseOpenBLAS）
        • 1.4.2.3 MKL的基准测试方法（UseMKL）

      • 1.4.3 目前的基准测试结果


  • 二、标量算法优化

    • 2.1 循环顺序由ijk改为ikj（TileRow）

      • 2.1.1 算法的实现

        • 2.1.1.1 矩阵清零

      • 2.1.2 基准测试方法
      • 2.1.3 基准测试结果

    • 2.2 性能为什么会有这么大提升呢？

      • 2.2.1 循环次序ijk（Basic）
      • 2.2.2 循环次序ikj（TileRow）

    • 2.3 使用Span及消除索引计算时的乘法开销（TileRowSpan）

      • 2.3.1 算法的实现
      • 2.3.2 基准测试方法
      • 2.3.3 基准测试结果

    • 2.4 使用托管指针进一步降低地址计算的开销（TileRowRef）

      • 2.4.1 算法的实现

        • 2.4.1.1 矩阵清零

      • 2.4.2 基准测试方法
      • 2.4.3 基准测试结果


  • 三、向量算法优化与并行加速

    • 3.1 使用SIMD硬件加速技术的向量算法（TileRowSimd）

      • 3.1.1 算法的实现
      • 3.1.2 基准测试方法
      • 3.1.3 基准测试结果

    • 3.2 基于TileRowSimd的并行加速（TileRowSimdParallel）

      • 3.2.1 基准测试方法
      • 3.2.2 基准测试结果

    • 3.3 使用FusedMultiplyAdd（TileRowSimdFmaBcl）

      • 3.1.1 算法的实现
      • 3.1.2 基准测试方法
      • 3.1.3 基准测试结果
      • 3.1.4 加上TileRowSimdFmaX86的测试

        • 3.1.4.1 TileRowSimdFmaX86的测试结果


    • 3.4 将乘加运算封装为函数且测试并行加速（TileRowSimdFma、TileRowSimdFmaParallel）

      • 3.4.1 乘加运算封装为函数
      • 3.4.2 矩阵乘法的代码
      • 3.4.3 基准测试方法
      • 3.4.4 基准测试结果


  • 四、矩阵分块优化

    • 4.1 4行、1倍向量宽度、ijk顺序的分块算法（BlockM4Nv1_ijk_M32、BlockM4Nv1_ijk_M32Parallel）

      • 4.1.1 公共方法的改进
      • 4.1.2 矩阵乘法内核算法
      • 4.1.3 矩阵乘法完整算法
      • 4.1.4 基准测试方法
      • 4.1.5 基准测试结果

    • 4.2 4行、1倍向量宽度、ikj顺序的分块算法（BlockM4Nv1_ikj_M32、BlockM4Nv1_ikj_M32Parallel、BlockM4Nv1_ikj_M4、BlockM4Nv1_ikj_M4Parallel）

      • 4.2.1 矩阵乘法内核算法
      • 4.2.2 矩阵乘法完整算法
      • 4.2.3 基准测试方法
      • 4.2.4 基准测试结果

    • 4.3 4行、3倍向量宽度、ijk顺序的分块算法（BlockM4Nv3_ikj_M32、BlockM4Nv3_ikj_M32Parallel、BlockM4Nv3_ikj_M4、BlockM4Nv3_ikj_M4Parallel）

      • 4.3.1 公共方法的改进
      • 4.3.2 矩阵乘法内核算法
      • 4.3.3 矩阵乘法完整算法
      • 4.3.4 基准测试方法
      • 4.3.5 基准测试结果

    • 4.4 4行、3倍512位向量宽度、ijk顺序的分块算法（BlockM4Nv3On512_ikj_M32、BlockM4Nv3On512_ikj_M32Parallel、BlockM4Nv3On512_ikj_M4、BlockM4Nv3On512_ikj_M4Parallel）

      • 4.4.1 公共方法的改进
      • 4.4.2 矩阵乘法内核算法
      • 4.4.3 矩阵乘法完整算法
      • 4.4.4 基准测试方法
      • 4.4.5 基准测试结果


  • 五、基准测试结果

    • 5.1 X86 架构 - .NET 9.0

      • 5.1.1 Single（单精度浮点型）
      • 5.1.2 Double（双精度浮点型）

    • 5.2 X86 架构 - .NET Framework

      • 5.2.1 Single（单精度浮点型）
      • 5.2.2 Double（双精度浮点型）

    • 5.3 Arm 架构 - .NET 9.0

      • 5.3.1 Single（单精度浮点型）
      • 5.3.2 Double（双精度浮点型）


  • 六、结语
  • 附录

[C#] 使用 .NET 的跨平台SIMD硬件加速技术，将 GEMM（通用矩阵乘法）程序性能提升1080倍，比肩 MKL、OpenBLAS
GEMM（General Matrix Multiply，通用矩阵乘法）是科学计算与深度学习等领域的核心算法。GEMM的性能越高，对科学计算与深度学习等领域的促进作用越大。对于CPU做GEMM运算，Intel MKL（Math Kernel Library。后面将统一简称为 MKL）、OpenBLAS 是第一梯队的行业翘楚。它们使用了并行、SIMD指令集、手动汇编优化等优化技术，使性能达到行业顶尖水平。
以前用 C# 开发的GEMM程序的性能，比MKL、OpenBLAS差得远，这是因为那时的 .NET 不支持SIMD硬件加速技术。从2014年开始， .NET 对SIMD硬件加速技术的支持越来越完善了。我潜心研究用该技术来改进 C# GEMM程序的性能，最近有了重大突破——对于1024尺寸矩阵的SGEMM（Single General Matrix Multiply，单精度浮点数通用矩阵乘法），我的算法比基础算法的性能提升1080倍，与 MKL、OpenBLAS的测试结果在同一梯队。
它是纯 C# 编写的托管代码程序，同一套源代码能在 X86、Arm 等架构上运行，性能均能达到第一梯队的水平。而且该算法能很容易地改造为DGEMM（Double General Matrix Multiply，双精度浮点数通用矩阵乘法）运算——因巧妙利用了 using指令，只用复制（SGEMM的）源代码，无需修改。

  • MKL 版本：2022.0.0.115
  • OpenBLAS 版本：2025.1.13
一、基本算法与测试办法

1.1 矩阵乘法的定义

通用矩阵乘法（General Matrix Multiply，GEMM）是计算两个矩阵 A 和 B 的乘积 C 的运算，要求 A 的列数等于 B 的行数。假设：

  • A 是 M * K 矩阵（M 行，K 列），
  • B 是 K * N 矩阵（K 行，N 列），
  • 则乘积 C = A * B 是 M * N 矩阵，其中每个元素 $C_{i,j}$ 为：$C_{i,j} = \sum_{k=1}^K A_{i,k} \cdot B_{k,j}$。即 A 的第 i 行与 B 的第 j 列的点积。
1.png

1.1.1 矩阵形状与运算复杂度


  • 时间复杂度：O(M * N * K)，共需 M * N * K 次乘法和加法运算。
  • 空间复杂度：输入矩阵需 O(M * K + K * N) 存储空间，输出矩阵需 O(M * N)。
若 M、N、K 都相同，即 A、B、C都是 N*N 的方阵时——

  • 时间复杂度：O(N^3)，共需 N^3 次乘法和加法运算。
  • 空间复杂度：输入矩阵需 O(2 * N^2) 存储空间，输出矩阵需 O(N^2)。
当 N 为 1024时，乘法总次数为：1024^3 = 1073741824 = 1 Giga。
浮点运算总次数是N^3 次乘法和加法运算，即 2 * 1024^3 = 2 * 1073741824 = 2 Giga。
此时（ N 为 1024时）可以比较方便的算出 每秒浮点运算次数 GFOPS (Giga Floating-point Operations Per Second)。假设程序耗时为 t 秒(s)，那么它每秒浮点运算次数为——
  1. 每秒浮点运算次数 = 浮点运算总次数 / 耗时 = 2 / t (GFOPS)
复制代码
1.2 C++ 开发的矩阵乘法最基础实现

根据上面的矩阵乘法定义，可以编程实现它。例如用 C++ 开发的矩阵乘法最基础实现如下。
  1. void gemmVanilla(const std::vector &matA, const std::vector &matB, std::vector &matC, size_t matSize) { // Clear matC. for (size_t i = 0; i < matC.size(); i++) { matC[i] = 0; } // Matrix matrix multiply. for (size_t i = 0; i < matSize; i++) { for (size_t j = 0; j < matSize; j++) { size_t cIdx = i * matSize + j; for (size_t k = 0; k < matSize; k++) { size_t aIdx = i * matSize + k; size_t bIdx = k * matSize + j; matC[cIdx] += matA[aIdx] * matB[bIdx]; } } } }
复制代码
这段代码出自 【深缓中字】为什么6层嵌套循环让这个算法快了120倍？。这个视频是学习矩阵算法优化的优秀教程，推荐读者观看。
这段代码使用 std::vector 类型来传递矩阵数据。虽然 std::vector 本身仅能存储一维数据，而这段代码隐式约定了其存储了 matSize * matSize 个数据，正好是二维的矩阵数据。由于 std::vector 是线性存储数据的，所以此时存储的矩阵，可看作行主序的矩阵。
这段代码的开头做了“清除矩阵C”（Clear matC）的工作。随后按照矩阵乘法定义，编写3层循环，实现了矩阵乘法的功能。
1.3 C# 开发的矩阵乘法最基础实现（Basic）

1.3.1 矩阵乘法的实现

将上述代码由 C++ 转为 C# 语言，便能得到 C# 版矩阵乘法。见下面的 StaticBasic 方法。
  1. // My type. using TMy = Single; partial class MultiplyMatrixStatic { ///  /// Basic on Array. ///  /// The number of rows in matrix A (矩阵A的行数). /// The number of columns in matrix B (矩阵B的列数). /// The number of columns in matrix A, or the number of rows in matrix B (矩阵A的列数, 或矩阵B的行数). /// Matrix A. /// Stride of A. /// Matrix B. /// Stride of B. /// Matrix C. /// Stride of C. public static void StaticBasic(int M, int N, int K, TMy[] A, int strideA, TMy[] B, int strideB, TMy[] C, int strideC) { // Matrix multiply. for (int i = 0; i < M; ++i) { for (int j = 0; j < N; ++j) { int cIdx = i * strideC + j; C[cIdx] = 0; for (int k = 0; k < K; ++k) { int aIdx = i * strideA + k; int bIdx = k * strideB + j; C[cIdx] += A[aIdx] * B[bIdx]; } } } } }
复制代码
上述代码位于 MultiplyMatrixStatic.Single.cs.
C# 版代码做了3个改造——

  • 使用using语句定义了别名（using TMy = Single）。这能便于将代码改造为DGEMM（Double General Matrix Multiply，双精度浮点数通用矩阵乘法），那时仅需将 TMy指向Double类型就行。
  • 简化了代码。移除了“清除矩阵C”（Clear matC）的工作，改为每次在 j 循环时将C矩阵的对应元素清零。
  • 修改了函数参数列表，参考 BLAS(Basic Linear Algebra Subprograms)的sgemm函数的参数列表。明确传递了 M、N、K 这3种尺寸参数，便于支持非方形的矩阵。且明确传递了各个矩阵的跨距（strideA），它类似 sgemm函数的 lda参数。
参数说明——

  • M：矩阵A和矩阵C的行数.
  • N：矩阵B和矩阵C的列数.
  • K：矩阵A和矩阵B的列数.
  • A：矩阵A. 即左矩阵.
  • strideA：矩阵A每一行的跨距.
  • B：矩阵B. 即右矩阵.
  • strideB ：矩阵B每一行的跨距.
  • C：矩阵C. 即结果矩阵.
  • strideC：矩阵C每一行的跨距.
由于二维数组的性能比较差，于是这里使用了一维数组（TMy[]）。随后像 C++ 代码那样，手动计算索引，从而实现了对二维矩阵数据的计算。
1.3.2 基准测试方法

本文统一使用 BenchmarkDotNet 做基准测试。
对于上面的StaticBasic方法，它的基准测试代码如下。
  1. [Benchmark(Baseline = true)] public void Basic() { StaticBasic(MatrixM, MatrixN, MatrixK, arrayA!, StrideA, arrayB!, StrideB, arrayC!, StrideC); if (CheckMode) { dstTMy = GetCheckSum(); baselineTMy = dstTMy; BenchmarkUtil.WriteItem("# Basic", string.Format("{0}", baselineTMy)); } }
复制代码
上述代码位于 MatrixNMultiplyBenchmark_Single.cs.
该方法内的第一行，调用了StaticBasic方法，从而实现了对它基准测试。MatrixM, MatrixN, MatrixK, arrayA 等参数都是事先初始化好的，这样能避免每次重新分配带来的开销。
后面的 CheckMode 分支，是在做矩阵运算结果检查。得益于现代CPU的分支预测能力，该分支不会影响基准测试。读者若对矩阵运算结果检查或arrayA初始化感兴趣的话，可自行查阅相关源代码。
注：本项目开启了 可空引用检查（Nullable），因arrayA等数组可能为空，若直接使用会有警告，但其实这些数组已经做好初始化了。于是可以在数组变量后面加“!”（感叹号），提醒编译器此时数组已经非空，以此消除警告。
1.4 增加MathNet、MKL、OpenBLAS的基准测试

本程序还引入了下列矩阵运算库，用于做基准测试成绩的对比。

  • MathNet 版本：6.0.0-beta2
  • OpenBLAS 版本：2025.1.13
  • MKL 版本：2022.0.0.115
1.4.1 引入库

在 nuget.org 上，可以找到这些库的 C# 包。

  • MathNet：使用“MathNet.Numerics”包。它由托管代码组成，支持在各种CPU架构、各种操作系统上运行。
  • MKL：使用“MKL.NET”等包。它仅提供了X86架构的原生动态库，支持在 Windows、MacOS等操作系统。
  • OpenBLAS：使用“OpenBlasSharp”等包。虽然OpenBLAS支持多种CPU架构，但该包仅提供了X86架构、Windows平台的原生动态库。
虽然使用 Visual Studio 的 Nuget包管理 功能可以添加这些包，但原生动态库仅支持部分平台。为了跨平台测试方便，于是我在 MatrixBenchmarkCs.csproj 里编写了条件分支来引入这些包。
  1. true false $(DefineConstants);USE_NATIVE_DLL     
复制代码
上面配置为——仅 .NET 9.0、Windows操作系统时，才使用原生动态库。若你想在其他平台上测试，可以修改 false 中的值。将它改为“true”。
版本号是在 Directory.Build.props 里定义的。
  1. 6.0.0-beta2 1.6.0 2022.0.0.115 2022.0.1.117 2022.0.0.105 0.3.4 2025.1.13
复制代码
1.4.2 这些库的基准测试方法

1.4.2.1 MathNet的基准测试方法（UseMathNet）

由于 MathNet 的矩阵运算函数不直接支持数组，而需要用它提供的矩阵类型。于是按它要求，定义了相关变量，并进行了初始化。
  1. protected MathNet.Numerics.LinearAlgebra.Matrix? matA; protected MathNet.Numerics.LinearAlgebra.Matrix? matB; protected MathNet.Numerics.LinearAlgebra.Matrix? matC; protected override void GlobalSetupCore() { base.GlobalSetupCore(); try { matA = MathNet.Numerics.LinearAlgebra.Matrix.Build.DenseOfRowMajor(MatrixM, MatrixK, arrayA); matB = MathNet.Numerics.LinearAlgebra.Matrix.Build.DenseOfRowMajor(MatrixK, MatrixN, arrayB); matC = MathNet.Numerics.LinearAlgebra.Matrix.Build.Dense(MatrixM, MatrixN); } catch (Exception ex) { BenchmarkUtil.WriteItem("# Setup-error", string.Format("{0}", ex.Message)); } }
复制代码
MathNet的基准测试的方法如下。
  1. [Benchmark] public void UseMathNet() { //MathNet.Numerics.Control.UseMultiThreading(); // Default is UseMultiThreading. matA!.Multiply(matB, matC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("UseMathNet"); } } //[Benchmark] public void UseMathNetSingleThread() { MathNet.Numerics.Control.UseSingleThread(); matA!.Multiply(matB, matC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("UseMathNet"); } }
复制代码
使用Matrix类型的Multiply方法进行矩阵乘法运算。
注意 MathNet 默认是执行多线程运算的。若想对比 MathNet的单线程时的性能，可以将UseMathNetSingleThread中 //[Benchmark] 的注释移除，并将UseMathNet中 //MathNet.Numerics.Control.UseMultiThreading 的注释移除。测试后会发现 UseMathNetSingleThread的成绩很差，且多线程的UseMathNet的性能也下降了。应该是 MathNet.Numerics.Control.UseMultiThreading 开销比较大，故本程序注释了UseMathNetSingleThread方法，仅测试多线程，不用切换。
由于 MathNet 的矩阵运算函数不直接支持数组，而是需要用它提供的矩阵类型，使用繁琐。CheckMode分支目前其实是不起作用的。因我们重点是做性能的基准测试，可忽略这一点。
1.4.2.2 OpenBLAS的基准测试方法（UseOpenBLAS）

OpenBLAS的基准测试的方法如下。
  1. [Benchmark] public unsafe void UseOpenBLAS() { fixed (TMy* pA = &arrayA![0], pB = &arrayB![0], pC = &arrayC![0]) { OpenBlasSharp.Blas.Sgemm(OpenBlasSharp.Order.RowMajor, OpenBlasSharp.Transpose.NoTrans, OpenBlasSharp.Transpose.NoTrans, MatrixM, MatrixN, MatrixK, 1, pA, StrideA, pB, StrideB, 0, pC, StrideC); } if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("UseOpenBLAS"); } }
复制代码
使用该库的Sgemm方法进行矩阵乘法运算。因它仅支持原生指针参数，于是需要使用 fixed 关键字来获取原生指针。
1.4.2.3 MKL的基准测试方法（UseMKL）

MKL的基准测试的方法如下。
  1. [Benchmark] public void UseMKL() { MKLNET.Blas.gemm(MKLNET.Layout.RowMajor, MKLNET.Trans.No, MKLNET.Trans.No, MatrixM, MatrixN, MatrixK, 1, arrayA!, StrideA, arrayB!, StrideB, 0, arrayC!, StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("UseMKL"); } }
复制代码
使用该库的gemm方法进行矩阵乘法运算。它支持数组参数，使用简单。
1.4.3 目前的基准测试结果

我的CPU是 “AMD Ryzen 7 7840H”，操作系统是 “Windows 11”。对于 1024 * 1024 矩阵的SGEMM，目前的基准测试结果如下。
  1. BenchmarkDotNet v0.15.2, Windows 11 (10.0.26100.4946/24H2/2024Update/HudsonValley) AMD Ryzen 7 7840H w/ Radeon 780M Graphics 3.80GHz, 1 CPU, 16 logical and 8 physical cores .NET SDK 9.0.301 [Host] : .NET 9.0.6 (9.0.625.26613), X64 AOT AVX-512F+CD+BW+DQ+VL+VBMI MediumRun : .NET 9.0.6 (9.0.625.26613), X64 RyuJIT AVX-512F+CD+BW+DQ+VL+VBMI
复制代码
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
UseMathNet102453,205,723.8 ns836,527.21 ns1,226,170.84 ns52,803,610.0 ns0.0110.007,575 B
UseOpenBLAS10242,932,070.9 ns117,359.68 ns160,643.26 ns2,864,316.2 ns0.0010.00NA
UseMKL10244,379,927.2 ns58,880.02 ns88,128.84 ns4,340,348.8 ns0.0010.00508 B
Basic 耗时4,712,905,103.3 ns，即 4.7 秒左右。而最快的UseOpenBLAS仅需 2,932,070.9 ns，即0.0029 秒左右，快了很多倍。
对测试结果进行仔细对比——

  • Basic: 耗时 4,712,905,103.3 ns, 故 GFOPS 为 2 / 4.7129051033 ≈ 0.42.
  • UseMathNet: 耗时 53,205,723.8 ns, 故 GFOPS 为 2 / 0.0532057238 ≈ 37.59, 性能是Basic的 4,712,905,103.3 / 53,205,723.8 ≈ 88.58 倍.
  • UseOpenBLAS: 耗时 2,932,070.9 ns, 故 GFOPS 为 2 / 0.0029320709 ≈ 682.11, 性能是Basic的 4,712,905,103.3 / 2,932,070.9 ≈ 1607.36 倍.
  • UseMKL: 耗时 4,379,927.2 ns, 故 GFOPS 为 2 / 0.0043799272 ≈ 456.63, 性能是Basic的 4,712,905,103.3 / 4,379,927.2 ≈ 1076.02 倍.
可以看到，基础算法比起MKL、OpenBLAS，有上千倍的性能差距。虽然差距大，但这也表示优化的空间很大。
二、标量算法优化

首先可以在标量算法的范畴内，多多考虑优化办法。
2.1 循环顺序由ijk改为ikj（TileRow）

对于矩阵乘法，有一种简单且实用的优化办法，它就是——循环顺序由ijk改为ikj。
这么简单的办法真的有效吗？实际测试一下就知道了。
2.1.1 算法的实现
  1. public static void StaticTileRow(int M, int N, int K, TMy[] A, int strideA, TMy[] B, int strideB, TMy[] C, int strideC) { // Clear matrix C. //C.AsSpan().Clear(); MatrixUtil.Fill((TMy)0, M, N, C, strideC); // Matrix multiply. for (int i = 0; i < M; ++i) { for (int k = 0; k < K; ++k) { int aIdx = i * strideA + k; for (int j = 0; j < N; ++j) { int bIdx = k * strideB + j; int cIdx = i * strideC + j; C[cIdx] += A[aIdx] * B[bIdx]; } } } }
复制代码
因现在循环顺序由ijk改为ikj，导致循环内不适合做 C矩阵元素的清零。于是需要在循环前对C矩阵整体做清零处理。
2.1.1.1 矩阵清零

对数组做清零处理，目前最方便的办法是将数组转为Span，再执行Clear方法。即 C.AsSpan().Clear()。
考虑到这个参数具有 strideC 参数，应根据该参数的含义，仅对范围内的数据执行清零。于是我写了一个工具方法 “MatrixUtil.Fill”来实现该功能，源代码如下。
  1. ///  /// Fill value (填充值). ///  /// The element type (元素的类型). /// The value (值). /// The number of rows in matrix (矩阵的行数). /// The number of columns in matrix (矩阵的列数). /// The matrix (矩阵). /// The stride of matrix. When it is 0, use cols (矩阵的跨距. 为 0 时 使用 cols). /// The start index of matrix (矩阵的开始索引). public static void Fill(T value, int rows, int cols, Span matrix, int stride = 0, int start = 0) { if (0 == stride) { stride = cols; } int idx0 = start; for (int i = 0; i < rows; i++) { Span span = matrix.Slice(idx0, cols); span.Fill(value); // Next. idx0 += stride; } }
复制代码
这个工具函数，是逐行对数组进行处理的。对于每一行的内容，转为Span做清零处理。当一行处理完毕时，会根据stride参数移动到下一行的起始索引。
2.1.2 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. [Benchmark] public void TileRow() { StaticTileRow(MatrixM, MatrixN, MatrixK, arrayA!, StrideA, arrayB!, StrideB, arrayC!, StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("TileRow"); } }
复制代码
2.1.3 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
TileRow1024780,626,575.0 ns5,970,550.33 ns8,562,785.59 ns778,917,250.0 ns0.1660.001,238 B
对测试结果进行对比——

  • Basic: 耗时 4,712,905,103.3 ns, 故 GFOPS 为 2 / 4.7129051033 ≈ 0.42.
  • TileRow: 耗时 780,626,575.0 ns, 故 GFOPS 为 2 / 0.7806265750 ≈ 2.56, 性能是Basic的 4,712,905,103.3 / 780,626,575.0 ≈ 6.04 倍.
简单来说，耗时从 4.7 秒，降低到 0.77 秒左右，速度提升了6.09倍。而改变仅是“循环顺序由ijk改为ikj”。
2.2 性能为什么会有这么大提升呢？

有些读者可能会疑惑——只是调换循环次序，甚至多了矩阵清零操作，性能为什么会有这么大提升呢？
现在来仔细分析一下。
2.2.1 循环次序ijk（Basic）

Basic算法的循环次序是ijk，即它的最内层循环是k的循环。此时分析一下3个矩阵的数据访问情况。

  • A（左矩阵）：优。最内层循环每处理一次，A的索引仅需+1，是连续的内存访问。
  • B（右矩阵）：劣。最内层循环每处理一次，B的索引仅需+strideB，即移动到下一行，不是连续的内存访问。
  • C（结果矩阵）：最优。最内层循环处理期间，C的索引保持不变。
因矩阵乘法的时间复杂度是 O(N^3)，使用立方体图片来表示矩阵乘法的运算过程，是比较直观的。例如下图表示了循环次序ijk时矩阵乘法的内循环。
2.png

A矩阵在左侧，旋转了90度；B矩阵在右侧；C矩阵在下侧。矩阵元素采取了“矩阵 列号,行号”的命名风格，例如“A1,2”表示“A矩阵的第1列第2行元素”。
绿色的长方体，表示最内层循环所影响的元素。目前展示了 C[0,0] = A[0,0]*B[0,0] + A[0,1]*B[1,0] + A[0,2]*B[2,0] + A[0,3]*B[3,0] + A[0,4]*B[4,0] 的计算。这个图中，A元素的投影与B元素的投影相交的位置，表示这2个元素值做乘法；而C元素的投影，正好覆盖了上述相交的位置，表示对这些元素值做累加。
同一个矩阵的相邻元素之间，虚线表示同一行的元素，实线表示同一列的元素。于是可以直观的看到——

  • 绿色的长方体在A矩阵上仅跨越了虚线，表示它是连续内存访问；
  • 而它在B矩阵上都是跨越实线，移动到下一行，这不是连续的内存访问。
2.2.2 循环次序ikj（TileRow）

TileRow算法的循环次序是ikj，即它的最内层循环是j的循环。此时分析一下3个矩阵的数据访问情况。

  • A（左矩阵）：最优。最内层循环处理期间，A的索引保持不变。
  • B（右矩阵）：优。最内层循环每处理一次，B的索引仅需+1，是连续的内存访问。
  • C（结果矩阵）：优。最内层循环每处理一次，C的索引仅需+1，是连续的内存访问。
我们再用立方体图片来观察循环次序ikj时矩阵乘法的内循环。
3.png

于是可以直观的看到——绿色的长方体在B矩阵、C矩阵上，均仅跨越了虚线，表示它们都是连续内存访问。
对比后可以发现，ikj比ijk次序在内存访问方面具有优势，这就是它性能更好的原因。
2.3 使用Span及消除索引计算时的乘法开销（TileRowSpan）

观察上一个算法，会发现在内循环里，除了矩阵乘法运算所需的乘法，还用到2次乘法——用乘法来计算索引（bIdx、cIdx）。也就是说上一个算法至少有 3 * N^3次乘法运算，而标准矩阵乘法是 N^3次乘法运算，它多了 2倍的 N^3次乘法运算。
消除索引计算时的乘法开销，能够提高性能。具体办法是根据循环的特点，在各级循环里做合适的累加操作，便能算出各索引的值。
另外，还可以用Span代替数组，使该方法不仅能支持数据等托管数据，还能支持非托管的数据。
2.3.1 算法的实现
  1. public static void StaticTileRowSpan(int M, int N, int K, Span A, int strideA, Span B, int strideB, Span C, int strideC) { // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, C, strideC); // Matrix multiply. int aIdx0 = 0; int cIdx0 = 0; for (int i = 0; i < M; ++i) { int aIdx = aIdx0; int bIdx0 = 0; for (int k = 0; k < K; ++k) { int bIdx = bIdx0; int cIdx = cIdx0; for (int j = 0; j < N; ++j) { C[cIdx] += A[aIdx] * B[bIdx]; ++bIdx; ++cIdx; } ++aIdx; bIdx0 += strideB; } aIdx0 += strideA; cIdx0 += strideC; } }
复制代码
可见，现在消除了索引计算时的乘法开销。仅在做矩阵计算时用了乘法。
2.3.2 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. [Benchmark] public void TileRowSpan() { StaticTileRowSpan(MatrixM, MatrixN, MatrixK, arrayA!, StrideA, arrayB!, StrideB, arrayC!, StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("TileRowSpan"); } }
复制代码
由于Span支持数组的隐式转换，所以基准测试方法的写法，与上一算法相同。
2.3.3 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
TileRowSpan1024613,236,262.1 ns6,326,567.15 ns9,273,400.86 ns610,494,600.0 ns0.1300.001,358 B
对测试结果进行对比——

  • TileRowSpan: 耗时 613,236,262.1 ns, 故 GFOPS 为 2 / 0.6132362621 ≈ 3.26, 性能是Basic的 4,712,905,103.3 / 613,236,262.1 ≈ 7.69 倍.
性能是Basic的 7.69 倍了。
2.4 使用托管指针进一步降低地址计算的开销（TileRowRef）

上面已经消除索引计算时的乘法开销了，但是使用索引访问数据的本身是有一定开销的。将索引访问，改造为指针访问，能进一步降低地址计算的开销。
C/C++ 的指针，在 C# 里叫做“原生指针”（Native pointer）或“非托管指针”（Unmanaged pointer），可以在unsafe关键字申明的“非安全代码”里使用，且需使用 fixed 语句锁定数据获取指针。
从 C# 7.3开始，可以用引用代替原生指针, 摆脱unsafe关键字与fixed语句，它也被称作 “托管指针”（Managed pointer）。其编程思路与原生指针是差不多的，只是个别地方需要换一种写法。具体办法可参考《C# 使用SIMD向量类型加速浮点数组求和运算（4）：用引用代替指针, 摆脱unsafe关键字，兼谈Unsafe类的使用》。
2.4.1 算法的实现
  1. public static void StaticTileRowRef(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC) { // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, ref C, strideC); // Matrix multiply. ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pC0 = ref C; for (int i = 0; i < M; ++i) { ref TMy pA = ref pA0; ref TMy pB0 = ref Unsafe.AsRef(in B); for (int k = 0; k < K; ++k) { TMy aValue = pA; ref TMy pB = ref pB0; ref TMy pC = ref pC0; for (int j = 0; j < N; ++j) { pC += aValue * pB; pB = ref Unsafe.Add(ref pB, 1); pC = ref Unsafe.Add(ref pC, 1); } pA = ref Unsafe.Add(ref pA, 1); pB0 = ref Unsafe.Add(ref pB0, strideB); } pA0 = ref Unsafe.Add(ref pA0, strideA); pC0 = ref Unsafe.Add(ref pC0, strideC); } }
复制代码
它的算法与上一算法（StaticTileRowSpan）完全相同，只是将索引计算，改为了托管指针的地址计算。
简单说明一下，“= ref”表示托管指针地址的赋值。例如 pA = ref Unsafe.Add(ref pA, 1)，相当于原生指针的 pA += 1。虽然语法啰嗦了一些，但功能与原生指针是一样的，且JIT编译生成的机器码是相同的，故性能也是相同的。
而“=”（等于符号）右边没有“ref”关键字时，这表示是普通赋值。会使用托管指针（ref）所指向的实际数据。例如 TMy aValue = pA，相当于原生指针的 TMy aValue = *pA。
Unsafe类的大部分方法仅支持“ref”，而不支持“ref readonly”。面对“ref readonly”，需使用Unsafe类的AsRef方法，对引用取消只读（类似C++的 const_cast）。例如 ref TMy pA0 = ref Unsafe.AsRef(in A)，相当于C++指针的 TMy* pA0 = const_cast(A)。
2.4.1.1 矩阵清零

因现在使用了托管指针，这就没法使用先前的Span版参数的“MatrixUtil.Fill”方法。于是新增了一个托管指针版参数的 “MatrixUtil.Fill”方法。
  1. public static void Fill(T value, nint rows, nint cols, ref T matrix, nint stride = 0) { if (0 == stride) { stride = cols; } int colsInt = (int)cols; #if NETSTANDARD2_1_OR_GREATER || NETCOREAPP2_1_OR_GREATER bool useSpan = (long)cols < int.MaxValue; #else bool useSpan = false; #endif ref T p0 = ref matrix; for (nint i = 0; i < rows; i++) { if (useSpan) { #if NETSTANDARD2_1_OR_GREATER || NETCOREAPP2_1_OR_GREATER Span span = MemoryMarshal.CreateSpan(ref p0, colsInt); span.Fill(value); #endif } else { ref T p = ref p0; for (nint j = 0; j < cols; j++) { p = value; p = ref Unsafe.Add(ref p, 1); } } // Next. p0 = ref Unsafe.Add(ref p0, stride); } }
复制代码
从 .NET Standard 2.1 或 .NET Core 2.1开始，能够将托管指针转为Span，随后便可以利用Span的Fill方法。而对于之前的版本，需自行使用托管指针来实现功能。
Span的索引是int类型的，故它最多仅能处理 int.MaxValue (2G) 范围内的数据。当 cols 大于该范围时，还得使用托管指针来处理。
2.4.2 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. [Benchmark] public void TileRowRef() { StaticTileRowRef(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("TileRowRef"); } }
复制代码
因StaticTileRowRef使用了托管指针，于是这里需传递托管指针。例如 ref arrayA![0] 表示传递“矩阵A首个元素的托管指针”。
2.4.3 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
TileRowRef1024332,063,973.2 ns4,375,391.05 ns6,275,055.62 ns329,953,500.0 ns0.0700.001,512 B
对测试结果进行对比——

  • TileRowRef: 耗时 332,063,973.2 ns, 故 GFOPS 为 2 / 0.3320639732 ≈ 6.02, 性能是Basic的 4,712,905,103.3 / 332,063,973.2 ≈ 14.19 倍.
性能是Basic的 14.19 倍了。可见，即使限制在标量算法的范畴内，优化也能取得不小成果的。
三、向量算法优化与并行加速

接下来介绍 SIMD向量算法、并行加速等优化手段。
3.1 使用SIMD硬件加速技术的向量算法（TileRowSimd）

回顾一下“[2.2.2 循环次序ikj（TileRow）](#2.2.2 循环次序ikj（TileRow）)”，可以注意到——每一次内循环中，A的索引保持不变，且B、C的索引仅需+1，是连续的内存访问。
这种情况，是能容易地改造为SIMD向量算法的。设向量宽度 W 为 Vector.Count 的值，那么SIMD向量算法的关键思路如下。

  • “A的索引保持不变”：使用 Vector vA = new Vector(pA)，就能根据当前A矩阵中的元素，广播到向量变量里（vA = pA）。
  • “B、C的索引仅需+1”：例如使用 ref Vector pB = ref Unsafe.As>(ref pB0)，就能直接将 TMy 的引用，重新解释为 Vector的引用。继而能够从 pB0 这个地址开始，依次将“W个连续值”加载到向量变量里（(pB) = pB0）。同理可这样加载pC。
随后执行 pC = Vector.Add(Vector.Multiply(vA, pB), pC)，便能对这“W个连续值”中的每个元素执行 pC += vA * pB 的操作。
标量算法每次仅能处理1个值，而SIMD向量算法每次能处理 W个值，故它的理论计算性能是标量算法的 W倍。
若读者对 .NET SIMD硬件加速技术不熟悉，可查阅 “C# 使用SIMD向量类型加速浮点数组求和运算（1）：使用Vector4、Vector” 等文章。
4.png

假设现在的向量类型是128位（如X86的SSE指令集、Arm的AdvSimd指令集），那么向量类型的字节尺寸为 128 / 8 = 16。因Single类型的尺寸是4个字节，那么向量类型里可以存储 16 / 4 = 4个元素。即 W（Vector.Count）的值为4，很适合处理N（B、C矩阵的列数）为4时的矩阵。例如上图，它就是4时的矩阵，绿色长方体在标量算法时需要写循环来实现，而用了SIMD向量算法后能一次性算出这4个值。
当N为W的整数倍时，可以多次执行上述的SIMD向量算法，使整行的数据处理完毕。若不是整数倍，则需专门处理尾部数据，下面的源代码里使用“尾部覆盖”的策略。
3.1.1 算法的实现
  1. public static void StaticTileRowSimd(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC) { if (N < Vector.Count || !Vector.IsHardwareAccelerated) { StaticTileRowRef(M, N, K, in A, strideA, in B, strideB, ref C, strideC); return; } int cntRem = N % Vector.Count; // Remainder count. int cntBlockRaw = N / Vector.Count; // Block count raw. int cntBlock = cntBlockRaw; if (0 == cntRem) { --cntBlock; // Use vCLast. } // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, ref C, strideC); // Matrix multiply. ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pC0 = ref C; for (int i = 0; i < M; ++i) { ref TMy pA = ref pA0; ref TMy pB0 = ref Unsafe.AsRef(in B); for (int k = 0; k < K; ++k) { TMy aValue = pA; Vector vA = new Vector(aValue); // Last. int pos = N - Vector.Count; ref Vector pBLast = ref Unsafe.As>(ref Unsafe.Add(ref pB0, pos)); ref Vector pCLast = ref Unsafe.As>(ref Unsafe.Add(ref pC0, pos)); Vector vCLast = Vector.Add(Vector.Multiply(vA, pBLast), pCLast); // SIMD for. if (cntBlock >= 0) { ref Vector pB = ref Unsafe.As>(ref pB0); ref Vector pC = ref Unsafe.As>(ref pC0); for (int j = 0; j < cntBlock; ++j) { pC = Vector.Add(Vector.Multiply(vA, pB), pC); // pC += vA * pB; pB = ref Unsafe.Add(ref pB, 1); pC = ref Unsafe.Add(ref pC, 1); } } pCLast = vCLast; // Overrride remainder items. // Next. pA = ref Unsafe.Add(ref pA, 1); pB0 = ref Unsafe.Add(ref pB0, strideB); } pA0 = ref Unsafe.Add(ref pA0, strideA); pC0 = ref Unsafe.Add(ref pC0, strideC); } }
复制代码
该方法的开头做了检查，若发现 N < Vector.Count（N小于向量宽度）或没有硬件加速时，便进行回退，调用标量算法 StaticTileRowRef。
注意 pB、pC 都是 Vector 类型的托管指针，所以在 Unsafe.Add(ref pB, 1)时，并不是对单个 TMy 元素移动指针（增加 sizeof(TMy)字节），而是对整块 Vector 移动指针（增加 sizeof(Vector)字节）。
3.1.2 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. [Benchmark] public void TileRowSimd() { StaticTileRowSimd(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("TileRowSimd"); } }
复制代码
3.1.3 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
TileRowSimd102471,736,249.5 ns2,768,453.89 ns3,880,986.03 ns69,775,500.0 ns0.0150.001,258 B
对测试结果进行对比——

  • TileRowSimd: 耗时 71,736,249.5 ns, 故 GFOPS 为 2 / 0.0717362495 ≈ 27.88, 性能是Basic的 4,712,905,103.3 / 71,736,249.5 ≈ 65.70 倍.
性能是Basic的 65.70 倍了。
3.2 基于TileRowSimd的并行加速（TileRowSimdParallel）

上面的代码是单线程的代码，仅能利用一个CPU核心，而现代CPU一般有多个核心。可使用并行计算，让每一个CPU核心都分摊计算工作，从而实现了加速。
3.2.1 基准测试方法

从 .NET Framework 4.0 开始，提供了 Parallel.For 方法，用它能方便的执行并行循环。
矩阵乘法是很容易并行化的，可以将最外层的i循环，分别交给不同的线程去处理，即可看作对C（结果）矩阵每一行去做并行处理。由于TileRowSimd算法的参数定义的很灵活，于是现在调整 M的值，以及 A、C矩阵的起始位置，便能实现让矩阵每一行去做并行处理。
所以我们不用修改TileRowSimd算法，仅需修改基准测试方法，便能增加它的并行加速的基准测试方法。
  1. [Benchmark] public void TileRowSimdParallel() { int M = MatrixM; bool allowParallel = (M >= 16) && (Environment.ProcessorCount > 1); if (allowParallel) { Parallel.For(0, M, i => { StaticTileRowSimd(1, MatrixN, MatrixK, ref arrayA![StrideA * i], StrideA, ref arrayB![0], StrideB, ref arrayC![StrideC * i], StrideC); }); } else { StaticTileRowSimd(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); } if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("TileRowSimdParallel"); } }
复制代码
在开始时做了一个检查，若M的值太小，或是仅有1个处理器时，便回退为单线程算法。
3.2.2 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
TileRowSimdParallel102411,576,096.3 ns447,403.26 ns669,652.19 ns11,727,672.3 ns0.0020.008,549 B
对测试结果进行对比——

  • TileRowSimdParallel: 耗时 11,576,096.3 ns, 故 GFOPS 为 2 / 0.0115760963 ≈ 172.77, 性能是Basic的 4,712,905,103.3 / 11,576,096.3 ≈ 407.12 倍.
性能是Basic的 407.12 倍了。
3.3 使用FusedMultiplyAdd（TileRowSimdFmaBcl）

从 .NET 9.0 开始，向量类型增加了 FusedMultiplyAdd（融合乘加）方法。它用于计算 (left * right) + addend。它的缩写是FMA。
现代X86及Arm处理器的向量指令集里，均支持FMA操作，可以在一条指令中完成FusedMultiplyAdd运算。它有这些优点——

  • 速度快。在1条指令中同时计算乘法与加法。若用它来处理数据的乘法与加法，理论上是分开计算的2倍。
  • 节省寄存器。分开计算时，需要再用一个寄存器来保存中间结果。而FMA因融合了乘法与加法，无需保存中间结果，故节省了1个寄存器。
  • 简化代码。分开计算时，需要分别调用向量类型的Multiply、Add方法。而改造后，仅需调用FusedMultiplyAdd这1个方法。
  • 精度高。通过单次舍入实现运算，相比分步计算具有更高精度。
矩阵运算中有很多乘法与加法，正好可以用上它。
3.1.1 算法的实现
  1. #if NET9_0_OR_GREATER public static void StaticTileRowSimdFmaBcl(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC) { if (N < Vector.Count || !Vector.IsHardwareAccelerated) { StaticTileRowRef(M, N, K, in A, strideA, in B, strideB, ref C, strideC); return; } int cntRem = N % Vector.Count; // Remainder count. int cntBlockRaw = N / Vector.Count; // Block count raw. int cntBlock = cntBlockRaw; if (0 == cntRem) { --cntBlock; // Use vCLast. } // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, ref C, strideC); // Matrix multiply. ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pC0 = ref C; for (int i = 0; i < M; ++i) { ref TMy pA = ref pA0; ref TMy pB0 = ref Unsafe.AsRef(in B); for (int k = 0; k < K; ++k) { TMy aValue = pA; Vector vA = new Vector(aValue); // Last. int pos = N - Vector.Count; ref Vector pBLast = ref Unsafe.As>(ref Unsafe.Add(ref pB0, pos)); ref Vector pCLast = ref Unsafe.As>(ref Unsafe.Add(ref pC0, pos)); Vector vCLast = Vector.FusedMultiplyAdd(vA, pBLast, pCLast); // SIMD for. if (cntBlock >= 0) { ref Vector pB = ref Unsafe.As>(ref pB0); ref Vector pC = ref Unsafe.As>(ref pC0); for (int j = 0; j < cntBlock; ++j) { pC = Vector.FusedMultiplyAdd(vA, pB, pC); // pC += vA * pB; pB = ref Unsafe.Add(ref pB, 1); pC = ref Unsafe.Add(ref pC, 1); } } pCLast = vCLast; // Overrride remainder items. // Next. pA = ref Unsafe.Add(ref pA, 1); pB0 = ref Unsafe.Add(ref pB0, strideB); } pA0 = ref Unsafe.Add(ref pA0, strideA); pC0 = ref Unsafe.Add(ref pC0, strideC); } } #endif // NET9_0_OR_GREATER
复制代码
代码改动很小，仅是将 Vector.Add(Vector.Multiply(vA, pB), pC) 改为了 Vector.FusedMultiplyAdd(vA, pB, pC)。
由于从 .NET 9.0 开始才支持 FusedMultiplyAdd，故需要使用条件编译判断一下 .NET 的版本。
3.1.2 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. #if NET9_0_OR_GREATER [Benchmark] public void TileRowSimdFmaBcl() { StaticTileRowSimdFmaBcl(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("TileRowSimdFmaBcl"); } } #endif // NET9_0_OR_GREATER
复制代码
3.1.3 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
TileRowSimd102471,736,249.5 ns2,768,453.89 ns3,880,986.03 ns69,775,500.0 ns0.0150.001,258 B
TileRowSimdFmaBcl102471,320,997.3 ns1,540,394.38 ns2,056,382.42 ns70,551,614.3 ns0.0150.001,262 B
对测试结果进行对比——

  • TileRowSimd: 耗时 71,736,249.5 ns, 故 GFOPS 为 2 / 0.0717362495 ≈ 27.88, 性能是Basic的 4,712,905,103.3 / 71,736,249.5 ≈ 65.70 倍.
  • TileRowSimdFmaBcl: 耗时 71,320,997.3 ns, 故 GFOPS 为 2 / 0.0713209973 ≈ 28.04, 性能是Basic的 4,712,905,103.3 / 71,320,997.3 ≈ 66.08 倍.
发现一个不对劲的地方——“理论上是分开计算的2倍”，但是为什么 TileRowSimdFmaBcl 与 TileRowSimd 的性能差不多，且个别时候还更慢了？
3.1.4 加上TileRowSimdFmaX86的测试

为了确认是不是 .NET 的编译质量问题，我又编写了一个基准测试函数 TileRowSimdFmaX86。它直接使用 X86的 Fma.MultiplyAdd 指令，使其不再受到编译质量的影响。源代码如下。
  1. #if NETCOREAPP3_0_OR_GREATER public static void StaticTileRowSimdFmaX86(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC) { if (N < Vector.Count || !Vector.IsHardwareAccelerated || !Fma.IsSupported) { StaticTileRowRef(M, N, K, in A, strideA, in B, strideB, ref C, strideC); return; } // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, ref C, strideC); // Matrix multiply. int cntRem = N % Vector.Count; // Remainder count. int cntBlockRaw = N / Vector.Count; // Block count raw. int cntBlock = cntBlockRaw; if (0 == cntRem) { --cntBlock; // Use vCLast. } ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pC0 = ref C; for (int i = 0; i < M; ++i) { ref TMy pA = ref pA0; ref TMy pB0 = ref Unsafe.AsRef(in B); for (int k = 0; k < K; ++k) { TMy aValue = pA; Vector vA = new Vector(aValue); // Last. int pos = N - Vector.Count; ref Vector pBLast = ref Unsafe.As>(ref Unsafe.Add(ref pB0, pos)); ref Vector pCLast = ref Unsafe.As>(ref Unsafe.Add(ref pC0, pos)); Vector vCLast = Vector.Add(Vector.Multiply(vA, pBLast), pCLast); // SIMD for. if (cntBlock >= 0) { ref Vector pB = ref Unsafe.As>(ref pB0); ref Vector pC = ref Unsafe.As>(ref pC0); for (int j = 0; j < cntBlock; ++j) { // pC += vA * pB; if (Vector.Count == Vector256.Count) { pC = Fma.MultiplyAdd(vA.AsVector256(), pB.AsVector256(), pC.AsVector256()).AsVector(); } else if (Vector.Count == Vector256.Count) { pC = Fma.MultiplyAdd(vA.AsVector128(), pB.AsVector128(), pC.AsVector128()).AsVector(); } else { pC = Vector.Add(Vector.Multiply(vA, pB), pC); } pB = ref Unsafe.Add(ref pB, 1); pC = ref Unsafe.Add(ref pC, 1); } } pCLast = vCLast; // Overrride remainder items. // Next. pA = ref Unsafe.Add(ref pA, 1); pB0 = ref Unsafe.Add(ref pB0, strideB); } pA0 = ref Unsafe.Add(ref pA0, strideA); pC0 = ref Unsafe.Add(ref pC0, strideC); } } #endif // NETCOREAPP3_0_OR_GREATER
复制代码
3.1.4.1 TileRowSimdFmaX86的测试结果

测试之后发现 TileRowSimdFmaBcl 的性能 与 TileRowSimdFmaX86 很接近。且均与 TileRowSimd 差不多。
这就表示先前并不是编译质量问题。而是因为现在的CPU流水线调度效率已经非常高了，能够掩盖“乘法与加法依次计算”时的延迟，使其与FMA的性能差不多。
测试结果的稍快与稍慢，仅是测试误差。
因耗时差不多，为了减少总测试时间，于是默认屏蔽了TileRowSimdFmaX86的基准测试。有兴趣的读者可以自行尝试。
3.4 将乘加运算封装为函数且测试并行加速（TileRowSimdFma、TileRowSimdFmaParallel）

虽然 FusedMultiplyAdd 很好用，能简化矩阵乘法的代码。可它需要 .NET 9.0 或更高版本，早期版本的 .NET 用不了它。
于是可以将乘加运算封装为函数， .NET 9.0使用FusedMultiplyAdd，早期的 .NET 便回退为乘法与加法。这样便能在矩阵运算法时放心的使用该函数了。
且对它进行并行加速的基准测试，检查它会不会拖慢性能。
3.4.1 乘加运算封装为函数

由于回退为分开计算乘法与及加法，不再是融合（Fused）计算。故这个函数命名为“MultiplyAdd”比较好。
建立 VectorHelper 类，增加“MultiplyAdd”方法。
  1. [MethodImpl(MethodImplOptions.AggressiveInlining)] public static Vector MultiplyAdd(Vector left, Vector right, Vector addend) { #if NET9_0_OR_GREATER return Vector.FusedMultiplyAdd(left, right, addend); #else return Vector.Add(addend, Vector.Multiply(left, right)); #endif // NET9_0_OR_GREATER }
复制代码
该函数设置了MethodImpl特性，并传递了 MethodImplOptions.AggressiveInlining 标识，使该函数能尽可能地内联。这能避免调用函数的开销与内存传值的开销。
3.4.2 矩阵乘法的代码

使用 VectorHelper.MultiplyAdd，能像FusedMultiplyAdd一样简化矩阵乘法的代码。
  1. public static void StaticTileRowSimdFma(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC) { if (N < Vector.Count || !Vector.IsHardwareAccelerated) { StaticTileRowRef(M, N, K, in A, strideA, in B, strideB, ref C, strideC); return; } int cntRem = N % Vector.Count; // Remainder count. int cntBlockRaw = N / Vector.Count; // Block count raw. int cntBlock = cntBlockRaw; if (0 == cntRem) { --cntBlock; // Use vCLast. } // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, ref C, strideC); // Matrix multiply. ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pC0 = ref C; for (int i = 0; i < M; ++i) { ref TMy pA = ref pA0; ref TMy pB0 = ref Unsafe.AsRef(in B); for (int k = 0; k < K; ++k) { TMy aValue = pA; Vector vA = new Vector(aValue); // Last. int pos = N - Vector.Count; ref Vector pBLast = ref Unsafe.As>(ref Unsafe.Add(ref pB0, pos)); ref Vector pCLast = ref Unsafe.As>(ref Unsafe.Add(ref pC0, pos)); //Vector vCLast = Vector.FusedMultiplyAdd(vA, pBLast, pCLast); Vector vCLast = VectorHelper.MultiplyAdd(vA, pBLast, pCLast); // SIMD for. if (cntBlock >= 0) { ref Vector pB = ref Unsafe.As>(ref pB0); ref Vector pC = ref Unsafe.As>(ref pC0); for (int j = 0; j < cntBlock; ++j) { //pC = Vector.FusedMultiplyAdd(vA, pB, pC); // pC += vA * pB; pC = VectorHelper.MultiplyAdd(vA, pB, pC); // pC += vA * pB; pB = ref Unsafe.Add(ref pB, 1); pC = ref Unsafe.Add(ref pC, 1); } } pCLast = vCLast; // Overrride remainder items. // Next. pA = ref Unsafe.Add(ref pA, 1); pB0 = ref Unsafe.Add(ref pB0, strideB); } pA0 = ref Unsafe.Add(ref pA0, strideA); pC0 = ref Unsafe.Add(ref pC0, strideC); } }
复制代码
3.4.3 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. [Benchmark] public void TileRowSimdFma() { StaticTileRowSimdFma(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("TileRowSimdFma"); } } public void TileRowSimdFmaParallel() { int M = MatrixM; bool allowParallel = (M >= 16) && (Environment.ProcessorCount > 1); if (allowParallel) { Parallel.For(0, M, i => { StaticTileRowSimdFma(1, MatrixN, MatrixK, ref arrayA![StrideA * i], StrideA, ref arrayB![0], StrideB, ref arrayC![StrideC * i], StrideC); }); } else { StaticTileRowSimdFma(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); } if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("TileRowSimdParallel"); } }
复制代码
上面同时给出了单线程与并行版的基准测试代码。
3.4.4 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
TileRowSimd102471,736,249.5 ns2,768,453.89 ns3,880,986.03 ns69,775,500.0 ns0.0150.001,258 B
TileRowSimdFmaBcl102471,320,997.3 ns1,540,394.38 ns2,056,382.42 ns70,551,614.3 ns0.0150.001,262 B
TileRowSimdFma102472,957,369.2 ns3,088,652.97 ns4,329,860.46 ns71,072,300.0 ns0.0150.001,262 B
TileRowSimdParallel102411,576,096.3 ns447,403.26 ns669,652.19 ns11,727,672.3 ns0.0020.008,549 B
TileRowSimdFmaParallel102412,234,503.5 ns337,579.48 ns505,273.12 ns12,306,358.6 ns0.0030.008,565 B
对测试结果进行对比——

  • TileRowSimd: 耗时 71,736,249.5 ns, 故 GFOPS 为 2 / 0.0717362495 ≈ 27.88, 性能是Basic的 4,712,905,103.3 / 71,736,249.5 ≈ 65.70 倍.
  • TileRowSimdFmaBcl: 耗时 71,320,997.3 ns, 故 GFOPS 为 2 / 0.0713209973 ≈ 28.04, 性能是Basic的 4,712,905,103.3 / 71,320,997.3 ≈ 66.08 倍.
  • TileRowSimdFma: 耗时 72,957,369.2 ns, 故 GFOPS 为 2 / 0.0729573692 ≈ 27.41, 性能是Basic的 4,712,905,103.3 / 72,957,369.2 ≈ 64.60 倍.
  • TileRowSimdParallel: 耗时 11,576,096.3 ns, 故 GFOPS 为 2 / 0.0115760963 ≈ 172.77, 性能是Basic的 4,712,905,103.3 / 11,576,096.3 ≈ 407.12 倍.
  • TileRowSimdFmaParallel: 耗时 12,234,503.5 ns, 故 GFOPS 为 2 / 0.0122345035 ≈ 163.47, 性能是Basic的 4,712,905,103.3 / 12,234,503.5 ≈ 385.21 倍.
2种Parallel算法的耗时差不多，且3种单线程算法的耗时也差不多，都在误差范围内。看来 VectorHelper.MultiplyAdd 与手动使用 FusedMultiplyAdd 的性能差不多。考虑到它具有简化代码等优点，故值得使用。
四、矩阵分块优化

要想进一步提高矩阵乘法的性能，需要对矩形进行分块（Block）。对矩形进行分块后，能利用缓存局部性减少内存访问次数，从而提高的了性能。
首先分块策略有许多种，且同一种分块策略对于不同的矩阵尺寸，效果是不一样的。即使是同一种分块策略、同一种矩阵尺寸，也会因为CPU微架构、高速缓存尺寸等情况的不同，而表现不同，甚至还有性能下降的。
需要经过大量测试，才能找到各个处理器的最佳分块策略。
4.1 4行、1倍向量宽度、ijk顺序的分块算法（BlockM4Nv1_ijk_M32、BlockM4Nv1_ijk_M32Parallel）

首先回顾一下 [3.1 使用SIMD硬件加速技术的向量算法（TileRowSimd）](#3.1 使用SIMD硬件加速技术的向量算法（TileRowSimd）)
5.png

设向量宽度 W 为 Vector.Count 的值，此时向量算法内循环是处理了 1行、W个列的数据。
若向量算法不再仅处理1行数据，而是处理4行数据，便能进一步提高性能。将它称作“4行、1倍向量宽度”算法。
6.png

上面的立方图图片展示了 M、N都是4时的矩阵，且采用的是128位向量，使W的值为4。于是，一笔“4行、1倍向量宽度”数据正好完成了水平方向1整层数据处理，随后图中的垂直方向（K的方向）循环对每一层结果的累加操作，便能算出完整的矩阵乘法结果。这种算法被称作“外积累加”算法，其中水平方向1整层的运算是向量的“外积”运算。该算法也被称作“kij”或“K-M-N”算法，因最外层是K循环，随后水平方向1整层的运算可看作“M循环嵌套N循环”。
若直接对大矩阵按“kij”的顺序做矩阵乘法的运算，需要频繁读写C矩阵的元素，导致性能很差。但对于“4行、1倍向量宽度”的算法来说，结果矩阵C正好可以用4个向量寄存器来存储。于是在k循环过程中，完全不需要对C矩阵进行读写，全是高效向量寄存器运算操作。直到循环处理完成后，才将向量寄存器中的数据保存到C矩阵，这便大大减少了内存访问的开销。
“4行、1倍向量宽度”算法的内循环 在行、列这2个方向上都是处理多个数据，导致经典的3层循环将难以处理完整矩阵乘法。此时可将内循环看作“4行、1倍向量宽度”的子矩阵，然后利用矩阵分块的策略来处理数据。
由于此时矩阵乘法的完整算法比较复杂，故拆分为3个部分来处理比较好——公共方法、内核算法、完整算法。
4.1.1 公共方法的改进

从 .NET 8.0 开始，Vector静态类增加了 LoadUnsafe、StoreUnsafe 方法，能便于向量类型的加载。
本程序为了对比测试各个 .NET 版本，于是需要对这些方法进行兼容处理。于是在 VectorHelper 类中增加了这些函数
  1. #if NET8_0_OR_GREATER #else #define VECTOR_WHERE_STRUCT// Since .NET8, Vector type not have `where T : struct`. #endif // NET8_0_OR_GREATER [MethodImpl(MethodImplOptions.AggressiveInlining)] public static Vector LoadUnsafe(ref readonly T source) #if VECTOR_WHERE_STRUCT where T : struct #endif // VECTOR_WHERE_STRUCT { #if NET8_0_OR_GREATER return Vector.LoadUnsafe(in source); #else return Unsafe.As>(ref Unsafe.AsRef(in source)); #endif // NET8_0_OR_GREATER } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void StoreUnsafe( #if !NET8_0_OR_GREATER this #endif // NET8_0_OR_GREATER Vector source, ref T destination) #if VECTOR_WHERE_STRUCT where T : struct #endif // VECTOR_WHERE_STRUCT { #if NET8_0_OR_GREATER Vector.StoreUnsafe(source, ref destination); #else Unsafe.As>(ref destination) = source; #endif // NET8_0_OR_GREATER }
复制代码
从 .NET 8.0 开始，向量类型取消了 struct 的约束。于是本类也遵从了这个规则。
4.1.2 矩阵乘法内核算法

将矩阵进行分块后，对于每个分块，专门用一个方法来处比较好。这种方法，不再时处理整个矩阵，而仅是单个分块，故一般称它为内核算法。
内核算法不用管理矩阵的清零操作，且一般仅需支持固定尺寸的数据。
为了便于测试各种矩阵分块策略，我这边会传递固定尺寸的整数倍的数据。于是m、n、k这3个参数也需做循环处理，只是不用考虑非整数倍数据。
  1. private static void GemmKernelM4Nv1(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC) { const int blockM = 4; int vectorWidth = Vector.Count; int blockN = vectorWidth * 1; if (0 != (M % blockM)) { throw new ArgumentException(string.Format("The M({0}) of block is not an integer multiple of {1}!", M, blockM)); } if (0 != (N % blockN)) { throw new ArgumentException(string.Format("The N({0}) of block is not an integer multiple of {1}!", N, blockN)); } // Matrix multiply. ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pB0 = ref Unsafe.AsRef(in B); ref TMy pC0 = ref C; for (int i = 0; i < M; i += blockM) { ref TMy pCLine = ref pC0; for (int j = 0; j < N; j += blockN) { ref TMy pC = ref pCLine; Vector c0; Vector c1; Vector c2; Vector c3; c0 = VectorHelper.LoadUnsafe(ref pC); pC = ref Unsafe.Add(ref pC, strideC); c1 = VectorHelper.LoadUnsafe(ref pC); pC = ref Unsafe.Add(ref pC, strideC); c2 = VectorHelper.LoadUnsafe(ref pC); pC = ref Unsafe.Add(ref pC, strideC); c3 = VectorHelper.LoadUnsafe(ref pC); // Add. ref TMy pALine = ref pA0; ref TMy pB = ref Unsafe.Add(ref pB0, j); for (int k = 0; k < K; ++k) { // pC += pA * pB; Vector b = VectorHelper.LoadUnsafe(ref pB); ref TMy pA = ref pALine; c0 = VectorHelper.MultiplyAdd(new Vector(pA), b, c0); pA = ref Unsafe.Add(ref pA, strideA); c1 = VectorHelper.MultiplyAdd(new Vector(pA), b, c1); pA = ref Unsafe.Add(ref pA, strideA); c2 = VectorHelper.MultiplyAdd(new Vector(pA), b, c2); pA = ref Unsafe.Add(ref pA, strideA); c3 = VectorHelper.MultiplyAdd(new Vector(pA), b, c3); // Next. pALine = ref Unsafe.Add(ref pALine, 1); pB = ref Unsafe.Add(ref pB, strideB); } // Store. pC = ref pCLine; VectorHelper.StoreUnsafe(c0, ref pC); pC = ref Unsafe.Add(ref pC, strideC); VectorHelper.StoreUnsafe(c1, ref pC); pC = ref Unsafe.Add(ref pC, strideC); VectorHelper.StoreUnsafe(c2, ref pC); pC = ref Unsafe.Add(ref pC, strideC); VectorHelper.StoreUnsafe(c3, ref pC); // Next. pCLine = ref Unsafe.Add(ref pCLine, blockN); } pA0 = ref Unsafe.Add(ref pA0, strideA * blockM); pC0 = ref Unsafe.Add(ref pC0, strideC * blockM); } }
复制代码
4.1.3 矩阵乘法完整算法

为了便于测试各种分块尺寸，完整算法增加了 blockM、blockN、blockK 参数。且为了便于并行计算，增加了 allowParallel 参数。
  1. internal unsafe static void StaticBlockM4Nv1_ijk(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, int blockM, int blockN, int blockK, bool allowParallel = false) { if (0 != (M % blockM)) { throw new ArgumentException(string.Format("The M({0}) is not an integer multiple of {1}!", M, blockM)); } if (0 != (N % blockN)) { throw new ArgumentException(string.Format("The N({0}) is not an integer multiple of {1}!", N, blockN)); } if (0 != (K % blockK)) { throw new ArgumentException(string.Format("The K({0}) is not an integer multiple of {1}!", K, blockK)); } // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, ref C, strideC); // Body. int strideB_blockK = strideB * blockK; ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pB0 = ref Unsafe.AsRef(in B); if (allowParallel) { fixed (TMy* pA = &pA0, pB = &pB0, pC = &C) { nint addressA = (nint)pA; nint addressB = (nint)pB; nint addressC = (nint)pC; int cnt = M / blockM; Parallel.For(0, cnt, idx => { int i = idx * blockM; ref TMy refA = ref Unsafe.AsRef((void*)addressA); ref TMy refB = ref Unsafe.AsRef((void*)addressB); ref TMy refC = ref Unsafe.AsRef((void*)addressC); RunByI(M, N, K, ref refA, strideA, ref refB, strideB, ref refC, strideC, i); }); } } else { for (int i = 0; i < M; i += blockM) { RunByI(M, N, K, ref pA0, strideA, ref pB0, strideB, ref C, strideC, i); } } void RunByI(int M, int N, int K, ref TMy A, int strideA, ref TMy B, int strideB, ref TMy C, int strideC, int i) { ref TMy pA0 = ref Unsafe.Add(ref A, strideA * i); ref TMy pCLine = ref Unsafe.Add(ref C, strideC * i); for (int j = 0; j < N; j += blockN) { ref TMy pALine = ref pA0; ref TMy pB = ref Unsafe.Add(ref B, j); for (int k = 0; k < K; k += blockK) { GemmKernelM4Nv1(blockM, blockN, blockK, in pALine, strideA, in pB, strideB, ref pCLine, strideC); // Next. pALine = ref Unsafe.Add(ref pALine, blockK); pB = ref Unsafe.Add(ref pB, strideB_blockK); } // Next. pCLine = ref Unsafe.Add(ref pCLine, blockN); } } } internal static void StaticBlockM4Nv1_ijk_M32(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, bool allowParallel = false) { StaticBlockM4Nv1_ijk(M, N, K, in A, strideA, in B, strideB, ref C, strideC, 32, 32, 32, allowParallel); }
复制代码
由于并行算法比较复杂，可以先看单线程的算法，即 allowParallel 为 false的分支。它先使用变量i对M进行循环处理，步长为blockM，逐个调用了 RunByI 方法。RunByI 是一个本地方法，它用了2层循环分别在 N、K方向处理数据，将每一个分块交给GemmKernelM4Nv1去处理。
随后再来看并行版算法，即 allowParallel 为 true的分支。它的原理并不复杂，就是使用 Parallel.For 并发的处理数据。但是为了避免 C# 安全检测“托管指针或原生指针都不能跨越异步方法的边界”报错，加了一些繁琐的转换，解决了指针传递问题。

  • 使用 fixed语句，锁定数据获取原生指针（如 pA）。
  • 使用强制类型转换，将原生指针转为nint（如 addressA）。
  • 在匿名方法内，使用强制类型转换，将nint转为原生指针（如 (void*)addressA）。
  • 最后使用 Unsafe.AsRef，将原生指针转为托管指针 （如 refA）。从而能够调用 RunByI 等方法了。
这些步骤用 C# 写起来是比较繁琐的，但这些不同类型的变量其实是同一个值，运行起来的实际开销很低。
并增加了 StaticBlockM4Nv1_ijk_M32 方法，专门对 32*32*32 的分块尺寸进行基准测试。
4.1.4 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. [Benchmark] public void BlockM4Nv1_ijk_M32() { StaticBlockM4Nv1_ijk_M32(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv1_ijk_M32"); } } [Benchmark] public void BlockM4Nv1_ijk_M32Parallel() { StaticBlockM4Nv1_ijk_M32(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC, true); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv1_ijk_M32Parallel"); } }
复制代码
这次写了2个基准测试方法，分别是单线程与并行计算的。由于 StaticBlockM4Nv1_ijk_M32 提供了allowParallel参数，于是基准测试方法写起来比较简单，仅需处理好 allowParallel 参数的传递。
4.1.5 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
BlockM4Nv1_ijk_M321024155,479,981.2 ns925,593.09 ns1,327,458.06 ns155,197,087.5 ns0.0330.00NA
BlockM4Nv1_ijk_M32Parallel102418,389,778.0 ns391,156.50 ns573,352.18 ns18,200,718.8 ns0.0040.00NA
对测试结果进行对比——

  • BlockM4Nv1_ijk_M32: 耗时 155,479,981.2 ns, 故 GFOPS 为 2 / 0.1554799812 ≈ 12.86, 性能是Basic的 4,712,905,103.3 / 155,479,981.2 ≈ 30.31 倍.
  • BlockM4Nv1_ijk_M32Parallel: 耗时 18,389,778.0 ns, 故 GFOPS 为 2 / 0.0183897780 ≈ 108.76, 性能是Basic的 4,712,905,103.3 / 18,389,778.0 ≈ 256.28 倍.
虽然这些算法也是比较快的，但它未能突破 未分块TileRowSimdParallel算法407.12倍的记录。这个一方面是分块策略的影响，另一方面与CPU微架构、高速缓存大小等因素有关。在一些Intel的CPU上测试过，BlockM4Nv1_ijk_M32Parallel与TileRowSimdParallel的性能差不多的，且矩阵越大时BlockM4Nv1_ijk_M32Parallel的领先程度更大。
4.2 4行、1倍向量宽度、ikj顺序的分块算法（BlockM4Nv1_ikj_M32、BlockM4Nv1_ikj_M32Parallel、BlockM4Nv1_ikj_M4、BlockM4Nv1_ikj_M4Parallel）

上面算法在进行分块时，用的是 ijk顺序。而ikj顺序的对连续内存访问更加友好，于是可以测试一下 。
4.2.1 矩阵乘法内核算法

矩阵乘法内核算法不用修改，仍然沿用。
4.2.2 矩阵乘法完整算法

将 ijk 顺序，改为 ikj 顺序。
  1. public const int CommonBlockK = 32; internal unsafe static void StaticBlockM4Nv1_ikj(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, int blockM, int blockN, int blockK, bool allowParallel = false) { if (0 != (M % blockM)) { throw new ArgumentException(string.Format("The M({0}) is not an integer multiple of {1}!", M, blockM)); } if (0 != (N % blockN)) { throw new ArgumentException(string.Format("The N({0}) is not an integer multiple of {1}!", N, blockN)); } if (0 != (K % blockK)) { throw new ArgumentException(string.Format("The K({0}) is not an integer multiple of {1}!", K, blockK)); } // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, ref C, strideC); // Body. int strideB_blockK = strideB * blockK; ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pB0 = ref Unsafe.AsRef(in B); if (allowParallel) { fixed (TMy* pA = &pA0, pB = &pB0, pC = &C) { nint addressA = (nint)pA; nint addressB = (nint)pB; nint addressC = (nint)pC; int cnt = M / blockM; Parallel.For(0, cnt, idx => { int i = idx * blockM; ref TMy refA = ref Unsafe.AsRef((void*)addressA); ref TMy refB = ref Unsafe.AsRef((void*)addressB); ref TMy refC = ref Unsafe.AsRef((void*)addressC); RunByI(M, N, K, ref refA, strideA, ref refB, strideB, ref refC, strideC, i); }); } } else { for (int i = 0; i < M; i += blockM) { RunByI(M, N, K, ref pA0, strideA, ref pB0, strideB, ref C, strideC, i); } } void RunByI(int M, int N, int K, ref TMy A, int strideA, ref TMy B, int strideB, ref TMy C, int strideC, int i) { ref TMy pALine = ref Unsafe.Add(ref A, strideA * i); ref TMy pCLine = ref Unsafe.Add(ref C, strideC * i); ref TMy pBLine = ref B; for (int k = 0; k < K; k += blockK) { ref TMy pB = ref pBLine; ref TMy pC = ref pCLine; for (int j = 0; j < N; j += blockN) { GemmKernelM4Nv1(blockM, blockN, blockK, in pALine, strideA, in pB, strideB, ref pC, strideC); // Next. pB = ref Unsafe.Add(ref pB, blockN); pC = ref Unsafe.Add(ref pC, blockN); } // Next. pALine = ref Unsafe.Add(ref pALine, blockK); pBLine = ref Unsafe.Add(ref pBLine, strideB_blockK); } } } internal static void StaticBlockM4Nv1_ikj_M32(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, bool allowParallel = false) { //int blockN = 32; int blockN = Vector.Count * 16; if (blockN > N) { blockN = N; } StaticBlockM4Nv1_ikj(M, N, K, in A, strideA, in B, strideB, ref C, strideC, 32, blockN, CommonBlockK, allowParallel); } internal static void StaticBlockM4Nv1_ikj_M4(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, bool allowParallel = false) { //int blockN = 32; int blockN = Vector.Count * 16; if (blockN > N) { blockN = N; } StaticBlockM4Nv1_ikj(M, N, K, in A, strideA, in B, strideB, ref C, strideC, 4, blockN, CommonBlockK, allowParallel); }
复制代码
这次提供了2种分块策略——

  • M4：M为4，N为16倍向量宽度。
  • M32：M为32，N为16倍向量宽度。
4.2.3 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. [Benchmark] public void BlockM4Nv1_ikj_M32() { StaticBlockM4Nv1_ikj_M32(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv1_ikj_M32"); } } [Benchmark] public void BlockM4Nv1_ikj_M32Parallel() { StaticBlockM4Nv1_ikj_M32(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC, true); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv1_ikj_M32Parallel"); } } [Benchmark] public void BlockM4Nv1_ikj_M4() { StaticBlockM4Nv1_ikj_M4(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv1_ikj_M4"); } } [Benchmark] public void BlockM4Nv1_ikj_M4Parallel() { StaticBlockM4Nv1_ikj_M4(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC, true); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv1_ikj_M4Parallel"); } }
复制代码
4.2.4 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
BlockM4Nv1_ikj_M321024158,203,231.7 ns1,098,377.31 ns1,643,999.58 ns157,511,900.0 ns0.0340.002,303 B
BlockM4Nv1_ikj_M32Parallel102415,417,779.4 ns528,972.18 ns775,360.62 ns15,101,640.6 ns0.0030.009,453 B
BlockM4Nv1_ikj_M41024159,175,464.8 ns2,913,579.61 ns4,084,432.04 ns161,698,000.0 ns0.0340.002,188 B
BlockM4Nv1_ikj_M4Parallel102417,944,453.9 ns584,658.22 ns856,984.51 ns17,830,587.5 ns0.0040.009,311 B
对测试结果进行对比——

  • BlockM4Nv1_ikj_M32: 耗时 158,203,231.7 ns, 故 GFOPS 为 2 / 0.1582032317 ≈ 12.64, 性能是Basic的 4,712,905,103.3 / 158,203,231.7 ≈ 29.79 倍.
  • BlockM4Nv1_ikj_M32Parallel: 耗时 15,417,779.4 ns, 故 GFOPS 为 2 / 0.0154177794 ≈ 129.72, 性能是Basic的 4,712,905,103.3 / 15,417,779.4 ≈ 305.68 倍.
  • BlockM4Nv1_ikj_M4: 耗时 159,175,464.8 ns, 故 GFOPS 为 2 / 0.1591754648 ≈ 12.56, 性能是Basic的 4,712,905,103.3 / 159,175,464.8 ≈ 29.61 倍.
  • BlockM4Nv1_ikj_M4Parallel: 耗时 17,944,453.9 ns, 故 GFOPS 为 2 / 0.0179444539 ≈ 111.46, 性能是Basic的 4,712,905,103.3 / 17,944,453.9 ≈ 262.64 倍.
比起 ijk顺序，ikj顺序要稍微快一些，尤其是在并行运算时。
另外还发现，大多数情况下 M32比M4快一些，证明了32的分块尺寸是合适的。有兴趣的读者可以自己试试其他的尺寸。
4.3 4行、3倍向量宽度、ijk顺序的分块算法（BlockM4Nv3_ikj_M32、BlockM4Nv3_ikj_M32Parallel、BlockM4Nv3_ikj_M4、BlockM4Nv3_ikj_M4Parallel）

sgemm_hsw 是一个用汇编语言编写的高性能矩阵乘法程序。它的关键改进是采用了3倍向量宽度的列数。
为什么3倍向量宽度是更好的选择呢？这与向量寄存器的数量有关。我们来梳理一下内循环需使用多少个向量寄存器——

  • 首先看C矩阵：因现在每次处理4行，于是在3倍向量宽度时，共需 4*3=12 个寄存器。
  • 随后看B矩阵：C矩阵处理下一行的数据时，B矩阵还在同一行，故仅需跟列数匹配的向量寄存器就行了。在3倍向量宽度时，共需 3个寄存器。
  • 随后看A矩阵：每次仅需读取一个浮点值，广播到1个向量寄存器里。故仅需1个向量寄存器。
由于FMA指令能够节省寄存器，没有保存中间结果的寄存器开销。故总共需要 12 + 3 + 1 = 16 的向量寄存器。
如今的X86处理器已经普及了AVX指令，它有16个向量寄存器。且Arm的AdvSimd指令集里，现在是最少16个向量寄存器。故上面的算法，正好能将这16个向量寄存器充分利用。
接下来我们用 C# 语言来实现这个算法，并进行基准测试。
4.3.1 公共方法的改进

为了简化"一次性读写3个向量的操作"，VectorHelper 增加了这些公共方法。
  1. [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void LoadX3Unsafe(ref readonly T source, out Vector data0, out Vector data1, out Vector data2) #if VECTOR_WHERE_STRUCT where T : struct #endif // VECTOR_WHERE_STRUCT { ref Vector p = ref Unsafe.As>(ref Unsafe.AsRef(in source)); data0 = p; data1 = Unsafe.Add(ref p, 1); data2 = Unsafe.Add(ref p, 2); } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void StoreX3Unsafe(ref T destination, Vector data0, Vector data1, Vector data2) #if VECTOR_WHERE_STRUCT where T : struct #endif // VECTOR_WHERE_STRUCT { ref Vector p = ref Unsafe.As>(ref destination); p = data0; Unsafe.Add(ref p, 1) = data1; Unsafe.Add(ref p, 2) = data2; }
复制代码
4.3.2 矩阵乘法内核算法

随后改进内核算法，让它一次性处理3倍向量宽度。
  1. private static void GemmKernelM4Nv3(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC) { const int blockM = 4; int vectorWidth = Vector.Count; int blockN = vectorWidth * 3; if (0 != (M % blockM)) { throw new ArgumentException(string.Format("The M({0}) of block is not an integer multiple of {1}!", M, blockM)); } if (0 != (N % blockN)) { throw new ArgumentException(string.Format("The N({0}) of block is not an integer multiple of {1}!", N, blockN)); } // Matrix multiply. ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pB0 = ref Unsafe.AsRef(in B); ref TMy pC0 = ref C; for (int i = 0; i < M; i += blockM) { ref TMy pCLine = ref pC0; for (int j = 0; j < N; j += blockN) { ref TMy pC = ref pCLine; Vector c0_0, c0_1, c0_2; Vector c1_0, c1_1, c1_2; Vector c2_0, c2_1, c2_2; Vector c3_0, c3_1, c3_2; VectorHelper.LoadX3Unsafe(ref pC, out c0_0, out c0_1, out c0_2); pC = ref Unsafe.Add(ref pC, strideC); VectorHelper.LoadX3Unsafe(ref pC, out c1_0, out c1_1, out c1_2); pC = ref Unsafe.Add(ref pC, strideC); VectorHelper.LoadX3Unsafe(ref pC, out c2_0, out c2_1, out c2_2); pC = ref Unsafe.Add(ref pC, strideC); VectorHelper.LoadX3Unsafe(ref pC, out c3_0, out c3_1, out c3_2); // Add. ref TMy pALine = ref pA0; ref TMy pB = ref Unsafe.Add(ref pB0, j); for (int k = 0; k < K; ++k) { // pC += pA * pB; Vector b0, b1, b2; Vector va; VectorHelper.LoadX3Unsafe(ref pB, out b0, out b1, out b2); ref TMy pA = ref pALine; va = new Vector(pA); c0_0 = VectorHelper.MultiplyAdd(va, b0, c0_0); c0_1 = VectorHelper.MultiplyAdd(va, b1, c0_1); c0_2 = VectorHelper.MultiplyAdd(va, b2, c0_2); pA = ref Unsafe.Add(ref pA, strideA); va = new Vector(pA); c1_0 = VectorHelper.MultiplyAdd(va, b0, c1_0); c1_1 = VectorHelper.MultiplyAdd(va, b1, c1_1); c1_2 = VectorHelper.MultiplyAdd(va, b2, c1_2); pA = ref Unsafe.Add(ref pA, strideA); va = new Vector(pA); c2_0 = VectorHelper.MultiplyAdd(va, b0, c2_0); c2_1 = VectorHelper.MultiplyAdd(va, b1, c2_1); c2_2 = VectorHelper.MultiplyAdd(va, b2, c2_2); pA = ref Unsafe.Add(ref pA, strideA); va = new Vector(pA); c3_0 = VectorHelper.MultiplyAdd(va, b0, c3_0); c3_1 = VectorHelper.MultiplyAdd(va, b1, c3_1); c3_2 = VectorHelper.MultiplyAdd(va, b2, c3_2); // Next. pALine = ref Unsafe.Add(ref pALine, 1); pB = ref Unsafe.Add(ref pB, strideB); } // Store. pC = ref pCLine; VectorHelper.StoreX3Unsafe(ref pC, c0_0, c0_1, c0_2); pC = ref Unsafe.Add(ref pC, strideC); VectorHelper.StoreX3Unsafe(ref pC, c1_0, c1_1, c1_2); pC = ref Unsafe.Add(ref pC, strideC); VectorHelper.StoreX3Unsafe(ref pC, c2_0, c2_1, c2_2); pC = ref Unsafe.Add(ref pC, strideC); VectorHelper.StoreX3Unsafe(ref pC, c3_0, c3_1, c3_2); // Next. pCLine = ref Unsafe.Add(ref pCLine, blockN); } pA0 = ref Unsafe.Add(ref pA0, strideA * blockM); pC0 = ref Unsafe.Add(ref pC0, strideC * blockM); } }
复制代码
4.3.3 矩阵乘法完整算法

因内核算法现在是 3倍向量宽度，导致它无法覆盖尺寸是“2的整数幂”的矩阵。例如 1024 无法被3整除。
而尺寸是“2的整数幂”的矩阵是经常使用的，故需想办法支持它。
可以采用这个办法——优先用GemmKernelM4Nv3处理数据，剩余数据用GemmKernelM4Nv1处理。于是RunByI的里面有2个j循环，且实现计算好了各自的循环次数。
  1. internal unsafe static void StaticBlockM4Nv3_ikj(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, int blockM, int blockN, int blockK, bool allowParallel = false) { int vectorWidth = Vector.Count; int vectorWidth3 = vectorWidth * 3; if (0 != (M % blockM)) { throw new ArgumentException(string.Format("The M({0}) is not an integer multiple of {1}!", M, blockM)); } if (blockN < vectorWidth3) { throw new ArgumentException(string.Format("The blockN({0}) is less {1}!", blockN, vectorWidth3)); } if (0 != (K % blockK)) { throw new ArgumentException(string.Format("The K({0}) is not an integer multiple of {1}!", K, blockK)); } int blockN3 = (blockN / vectorWidth3) * vectorWidth3; int n; int nEnd3; int nEnd1; n = N / blockN3; nEnd3 = n * blockN3; n = (N - nEnd3) / vectorWidth; nEnd1 = nEnd3 + n * vectorWidth; if (nEnd1 != N) { throw new ArgumentException(string.Format("The N({0}) is not an integer multiple of {1}!", N, vectorWidth)); } // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, ref C, strideC); // Body. int strideB_blockK = strideB * blockK; ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pB0 = ref Unsafe.AsRef(in B); if (allowParallel) { fixed (TMy* pA = &pA0, pB = &pB0, pC = &C) { nint addressA = (nint)pA; nint addressB = (nint)pB; nint addressC = (nint)pC; int cnt = M / blockM; Parallel.For(0, cnt, idx => { int i = idx * blockM; ref TMy refA = ref Unsafe.AsRef((void*)addressA); ref TMy refB = ref Unsafe.AsRef((void*)addressB); ref TMy refC = ref Unsafe.AsRef((void*)addressC); RunByI(M, N, K, ref refA, strideA, ref refB, strideB, ref refC, strideC, i); }); } } else { for (int i = 0; i < M; i += blockM) { RunByI(M, N, K, ref pA0, strideA, ref pB0, strideB, ref C, strideC, i); } } void RunByI(int M, int N, int K, ref TMy A, int strideA, ref TMy B, int strideB, ref TMy C, int strideC, int i) { ref TMy pALine = ref Unsafe.Add(ref A, strideA * i); ref TMy pCLine = ref Unsafe.Add(ref C, strideC * i); ref TMy pBLine = ref B; for (int k = 0; k < K; k += blockK) { ref TMy pB = ref pBLine; ref TMy pC = ref pCLine; int j = 0; for (; j < nEnd3; j += blockN3) { GemmKernelM4Nv3(blockM, blockN3, blockK, in pALine, strideA, in pB, strideB, ref pC, strideC); // Next. pB = ref Unsafe.Add(ref pB, blockN3); pC = ref Unsafe.Add(ref pC, blockN3); } for (; j < nEnd1; j += vectorWidth) { GemmKernelM4Nv1(blockM, vectorWidth, blockK, in pALine, strideA, in pB, strideB, ref pC, strideC); // Next. pB = ref Unsafe.Add(ref pB, vectorWidth); pC = ref Unsafe.Add(ref pC, vectorWidth); } // Next. pALine = ref Unsafe.Add(ref pALine, blockK); pBLine = ref Unsafe.Add(ref pBLine, strideB_blockK); } } } internal static void StaticBlockM4Nv3_ikj_M32(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, bool allowParallel = false) { //int blockN = 32; int blockN = Vector.Count * 3 * 4; if (blockN > N) { blockN = N; } StaticBlockM4Nv3_ikj(M, N, K, in A, strideA, in B, strideB, ref C, strideC, 32, blockN, CommonBlockK, allowParallel); } internal static void StaticBlockM4Nv3_ikj_M4(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, bool allowParallel = false) { //int blockN = 32; int blockN = Vector.Count * 3 * 4; if (blockN > N) { blockN = N; } StaticBlockM4Nv3_ikj(M, N, K, in A, strideA, in B, strideB, ref C, strideC, 4, blockN, CommonBlockK, allowParallel); }
复制代码
此时也提供了 M4、M32这2种分块策略的方法。
4.3.4 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. [Benchmark] public void BlockM4Nv3_ikj_M4() { StaticBlockM4Nv3_ikj_M4(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv3_ikj_M4"); } } [Benchmark] public void BlockM4Nv3_ikj_M4Parallel() { StaticBlockM4Nv3_ikj_M4(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC, true); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv3_ikj_M4Parallel"); } } [Benchmark] public void BlockM4Nv3_ikj_M32() { StaticBlockM4Nv3_ikj_M32(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv3_ikj_M32"); } } [Benchmark] public void BlockM4Nv3_ikj_M32Parallel() { StaticBlockM4Nv3_ikj_M32(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC, true); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv3_ikj_M32Parallel"); } }
复制代码
4.3.5 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
BlockM4Nv3_ikj_M4102452,756,227.9 ns1,133,938.16 ns1,626,260.36 ns53,519,345.5 ns0.0110.003,766 B
BlockM4Nv3_ikj_M4Parallel10246,234,993.0 ns450,007.94 ns645,388.00 ns5,976,738.7 ns0.0010.0010,419 B
BlockM4Nv3_ikj_M32102451,007,023.1 ns1,010,156.07 ns1,448,735.77 ns51,688,581.8 ns0.0110.003,819 B
BlockM4Nv3_ikj_M32Parallel10245,851,851.6 ns266,588.14 ns390,761.47 ns5,628,710.9 ns0.0010.0010,563 B
对测试结果进行对比——

  • BlockM4Nv3_ikj_M4: 耗时 52,756,227.9 ns, 故 GFOPS 为 2 / 0.0527562279 ≈ 37.91, 性能是Basic的 4,712,905,103.3 / 52,756,227.9 ≈ 89.33 倍.
  • BlockM4Nv3_ikj_M4Parallel: 耗时 6,234,993.0 ns, 故 GFOPS 为 2 / 0.0062349930 ≈ 320.77, 性能是Basic的 4,712,905,103.3 / 6,234,993.0 ≈ 755.88 倍.
  • BlockM4Nv3_ikj_M32: 耗时 51,007,023.1 ns, 故 GFOPS 为 2 / 0.0510070231 ≈ 39.21, 性能是Basic的 4,712,905,103.3 / 51,007,023.1 ≈ 92.40 倍.
  • BlockM4Nv3_ikj_M32Parallel: 耗时 5,851,851.6 ns, 故 GFOPS 为 2 / 0.0058518516 ≈ 341.77, 性能是Basic的 4,712,905,103.3 / 5,851,851.6 ≈ 805.37 倍.
性能最好的是 BlockM4Nv3_ikj_M32Parallel，它性能是Basic的 805.37 倍。且它突破了 TileRowSimdFmaParallel 407.12倍的记录。
4.4 4行、3倍512位向量宽度、ijk顺序的分块算法（BlockM4Nv3On512_ikj_M32、BlockM4Nv3On512_ikj_M32Parallel、BlockM4Nv3On512_ikj_M4、BlockM4Nv3On512_ikj_M4Parallel）

目前 MKL、OpenBLAS 都使用了 AVX-512指令集优化。
在 .NET 8.0 开始，C# 也能是用 AVX-512指令集了。大多数情况时，使用 Vector512 静态类的方法，会更加方便——对于大多数代码，仅需将 Vector 替换为 Vector512。
4.4.1 公共方法的改进

Vector512 也需要 LoadX3Unsafe等公共方法。于是新增了Vector512Helper类，增加了这些公共方法。
  1. [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void LoadX3Unsafe(ref readonly T source, out Vector512 data0, out Vector512 data1, out Vector512 data2) #if VECTOR_WHERE_STRUCT where T : struct #endif // VECTOR_WHERE_STRUCT { ref Vector512 p = ref Unsafe.As>(ref Unsafe.AsRef(in source)); data0 = p; data1 = Unsafe.Add(ref p, 1); data2 = Unsafe.Add(ref p, 2); } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static void StoreX3Unsafe(ref T destination, Vector512 data0, Vector512 data1, Vector512 data2) #if VECTOR_WHERE_STRUCT where T : struct #endif // VECTOR_WHERE_STRUCT { ref Vector512 p = ref Unsafe.As>(ref destination); p = data0; Unsafe.Add(ref p, 1) = data1; Unsafe.Add(ref p, 2) = data2; } [MethodImpl(MethodImplOptions.AggressiveInlining)] public static Vector512 MultiplyAdd(Vector512 left, Vector512 right, Vector512 addend) { #if NET9_0_OR_GREATER return Vector512.FusedMultiplyAdd(left, right, addend); #else return Vector512.Add(addend, Vector512.Multiply(left, right)); #endif // NET9_0_OR_GREATER }
复制代码
4.4.2 矩阵乘法内核算法

仅需将 Vector 改为 Vector512，便能实现 GemmKernelM4Nv3On512。
  1. private static void GemmKernelM4Nv3On512(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC) { const int blockM = 4; int vectorWidth = Vector512.Count; int blockN = vectorWidth * 3; if (BenchmarkUtil.IsLastRun) { Volatile.Write(ref blockN, 0); //Debugger.Break(); } blockN = vectorWidth * 3; if (0 != (M % blockM)) { throw new ArgumentException(string.Format("The M({0}) of block is not an integer multiple of {1}!", M, blockM)); } if (0 != (N % blockN)) { throw new ArgumentException(string.Format("The N({0}) of block is not an integer multiple of {1}!", N, blockN)); } // Matrix multiply. ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pB0 = ref Unsafe.AsRef(in B); ref TMy pC0 = ref C; for (int i = 0; i < M; i += blockM) { ref TMy pCLine = ref pC0; for (int j = 0; j < N; j += blockN) { ref TMy pC = ref pCLine; Vector512 c0_0, c0_1, c0_2; Vector512 c1_0, c1_1, c1_2; Vector512 c2_0, c2_1, c2_2; Vector512 c3_0, c3_1, c3_2; Vector512Helper.LoadX3Unsafe(ref pC, out c0_0, out c0_1, out c0_2); pC = ref Unsafe.Add(ref pC, strideC); Vector512Helper.LoadX3Unsafe(ref pC, out c1_0, out c1_1, out c1_2); pC = ref Unsafe.Add(ref pC, strideC); Vector512Helper.LoadX3Unsafe(ref pC, out c2_0, out c2_1, out c2_2); pC = ref Unsafe.Add(ref pC, strideC); Vector512Helper.LoadX3Unsafe(ref pC, out c3_0, out c3_1, out c3_2); // Add. ref TMy pALine = ref pA0; ref TMy pB = ref Unsafe.Add(ref pB0, j); for (int k = 0; k < K; ++k) { // pC += pA * pB; Vector512 b0, b1, b2; Vector512 va; Vector512Helper.LoadX3Unsafe(ref pB, out b0, out b1, out b2); ref TMy pA = ref pALine; va = Vector512.Create(pA); c0_0 = Vector512Helper.MultiplyAdd(va, b0, c0_0); c0_1 = Vector512Helper.MultiplyAdd(va, b1, c0_1); c0_2 = Vector512Helper.MultiplyAdd(va, b2, c0_2); pA = ref Unsafe.Add(ref pA, strideA); va = Vector512.Create(pA); c1_0 = Vector512Helper.MultiplyAdd(va, b0, c1_0); c1_1 = Vector512Helper.MultiplyAdd(va, b1, c1_1); c1_2 = Vector512Helper.MultiplyAdd(va, b2, c1_2); pA = ref Unsafe.Add(ref pA, strideA); va = Vector512.Create(pA); c2_0 = Vector512Helper.MultiplyAdd(va, b0, c2_0); c2_1 = Vector512Helper.MultiplyAdd(va, b1, c2_1); c2_2 = Vector512Helper.MultiplyAdd(va, b2, c2_2); pA = ref Unsafe.Add(ref pA, strideA); va = Vector512.Create(pA); c3_0 = Vector512Helper.MultiplyAdd(va, b0, c3_0); c3_1 = Vector512Helper.MultiplyAdd(va, b1, c3_1); c3_2 = Vector512Helper.MultiplyAdd(va, b2, c3_2); // Next. pALine = ref Unsafe.Add(ref pALine, 1); pB = ref Unsafe.Add(ref pB, strideB); } // Store. pC = ref pCLine; Vector512Helper.StoreX3Unsafe(ref pC, c0_0, c0_1, c0_2); pC = ref Unsafe.Add(ref pC, strideC); Vector512Helper.StoreX3Unsafe(ref pC, c1_0, c1_1, c1_2); pC = ref Unsafe.Add(ref pC, strideC); Vector512Helper.StoreX3Unsafe(ref pC, c2_0, c2_1, c2_2); pC = ref Unsafe.Add(ref pC, strideC); Vector512Helper.StoreX3Unsafe(ref pC, c3_0, c3_1, c3_2); // Next. pCLine = ref Unsafe.Add(ref pCLine, blockN); } pA0 = ref Unsafe.Add(ref pA0, strideA * blockM); pC0 = ref Unsafe.Add(ref pC0, strideC * blockM); } }
复制代码
4.4.3 矩阵乘法完整算法

可以将大片区域的数据交给 GemmKernelM4Nv3On512 处理，而尾部剩余交给 GemmKernelM4Nv1 处理，这便能支持各种“2的幂次”的矩阵。
  1. internal unsafe static void StaticBlockM4Nv3On512_ikj(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, int blockM, int blockN, int blockK, bool allowParallel = false) { if (!Vector512.IsHardwareAccelerated || !Vector512.IsSupported) { throw new NotSupportedException("Vector512.IsHardwareAccelerated is false!"); } int vectorWidth = Vector.Count; int vectorWidth3 = Vector512.Count * 3; if (0 != (M % blockM)) { throw new ArgumentException(string.Format("The M({0}) is not an integer multiple of {1}!", M, blockM)); } if (blockN < vectorWidth3) { throw new ArgumentException(string.Format("The blockN({0}) is less {1}!", blockN, vectorWidth3)); } if (0 != (K % blockK)) { throw new ArgumentException(string.Format("The K({0}) is not an integer multiple of {1}!", K, blockK)); } int blockN3 = (blockN / vectorWidth3) * vectorWidth3; int n; int nEnd3; int nEnd1; n = N / blockN3; nEnd3 = n * blockN3; n = (N - nEnd3) / vectorWidth; nEnd1 = nEnd3 + n * vectorWidth; if (nEnd1 != N) { throw new ArgumentException(string.Format("The N({0}) is not an integer multiple of {1}!", N, vectorWidth)); } // Clear matrix C. MatrixUtil.Fill((TMy)0, M, N, ref C, strideC); // Body. int strideB_blockK = strideB * blockK; ref TMy pA0 = ref Unsafe.AsRef(in A); ref TMy pB0 = ref Unsafe.AsRef(in B); if (allowParallel) { fixed (TMy* pA = &pA0, pB = &pB0, pC = &C) { nint addressA = (nint)pA; nint addressB = (nint)pB; nint addressC = (nint)pC; int cnt = M / blockM; Parallel.For(0, cnt, idx => { int i = idx * blockM; ref TMy refA = ref Unsafe.AsRef((void*)addressA); ref TMy refB = ref Unsafe.AsRef((void*)addressB); ref TMy refC = ref Unsafe.AsRef((void*)addressC); RunByI(M, N, K, ref refA, strideA, ref refB, strideB, ref refC, strideC, i); }); } } else { for (int i = 0; i < M; i += blockM) { RunByI(M, N, K, ref pA0, strideA, ref pB0, strideB, ref C, strideC, i); } } void RunByI(int M, int N, int K, ref TMy A, int strideA, ref TMy B, int strideB, ref TMy C, int strideC, int i) { ref TMy pALine = ref Unsafe.Add(ref A, strideA * i); ref TMy pCLine = ref Unsafe.Add(ref C, strideC * i); ref TMy pBLine = ref B; for (int k = 0; k < K; k += blockK) { ref TMy pB = ref pBLine; ref TMy pC = ref pCLine; int j = 0; for (; j < nEnd3; j += blockN3) { GemmKernelM4Nv3On512(blockM, blockN3, blockK, in pALine, strideA, in pB, strideB, ref pC, strideC); // Next. pB = ref Unsafe.Add(ref pB, blockN3); pC = ref Unsafe.Add(ref pC, blockN3); } for (; j < nEnd1; j += vectorWidth) { GemmKernelM4Nv1(blockM, vectorWidth, blockK, in pALine, strideA, in pB, strideB, ref pC, strideC); // Next. pB = ref Unsafe.Add(ref pB, vectorWidth); pC = ref Unsafe.Add(ref pC, vectorWidth); } // Next. pALine = ref Unsafe.Add(ref pALine, blockK); pBLine = ref Unsafe.Add(ref pBLine, strideB_blockK); } } } internal static void StaticBlockM4Nv3On512_ikj_M32(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, bool allowParallel = false) { //int blockN = 32; int blockN = Vector512.Count * 3 * 4; if (blockN > N) { blockN = N; } StaticBlockM4Nv3On512_ikj(M, N, K, in A, strideA, in B, strideB, ref C, strideC, 32, blockN, CommonBlockK, allowParallel); } internal static void StaticBlockM4Nv3On512_ikj_M4(int M, int N, int K, ref readonly TMy A, int strideA, ref readonly TMy B, int strideB, ref TMy C, int strideC, bool allowParallel = false) { //int blockN = 32; int blockN = Vector512.Count * 3 * 4; if (blockN > N) { blockN = N; } StaticBlockM4Nv3On512_ikj(M, N, K, in A, strideA, in B, strideB, ref C, strideC, 4, blockN, CommonBlockK, allowParallel); }
复制代码
4.4.4 基准测试方法

对于上面的方法，它的基准测试代码如下。
  1. #if NET8_0_OR_GREATER [Benchmark] public void BlockM4Nv3On512_ikj_M4() { StaticBlockM4Nv3On512_ikj_M4(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv3On512_ikj_M4"); } } [Benchmark] public void BlockM4Nv3On512_ikj_M4Parallel() { StaticBlockM4Nv3On512_ikj_M4(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC, true); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv3On512_ikj_M4Parallel"); } } [Benchmark] public void BlockM4Nv3On512_ikj_M32() { if (BenchmarkUtil.IsLastRun) { Volatile.Write(ref dstTMy, 0); //Debugger.Break(); } StaticBlockM4Nv3On512_ikj_M32(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv3On512_ikj_M32"); } } [Benchmark] public void BlockM4Nv3On512_ikj_M32Parallel() { StaticBlockM4Nv3On512_ikj_M32(MatrixM, MatrixN, MatrixK, ref arrayA![0], StrideA, ref arrayB![0], StrideB, ref arrayC![0], StrideC, true); if (CheckMode) { dstTMy = GetCheckSum(); CheckResult("BlockM4Nv3On512_ikj_M32Parallel"); } } #endif // NET8_0_OR_GREATER
复制代码
4.4.5 基准测试结果

基准测试结果如下。
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
BlockM4Nv3On512_ikj_M4102444,685,858.1 ns4,745,584.17 ns6,652,646.79 ns41,372,192.9 ns0.0090.003,835 B
BlockM4Nv3On512_ikj_M4Parallel10244,842,750.9 ns306,368.06 ns429,485.26 ns4,669,524.2 ns0.0010.0010,499 B
BlockM4Nv3On512_ikj_M32102432,179,510.0 ns206,091.41 ns275,126.14 ns32,157,606.2 ns0.0070.003,666 B
BlockM4Nv3On512_ikj_M32Parallel10244,363,112.2 ns73,021.25 ns97,481.28 ns4,339,956.2 ns0.0010.0010,848 B
UseMathNet102453,205,723.8 ns836,527.21 ns1,226,170.84 ns52,803,610.0 ns0.0110.007,575 B
UseOpenBLAS10242,932,070.9 ns117,359.68 ns160,643.26 ns2,864,316.2 ns0.0010.00NA
UseMKL10244,379,927.2 ns58,880.02 ns88,128.84 ns4,340,348.8 ns0.0010.00508 B
对测试结果进行对比——

  • BlockM4Nv3On512_ikj_M4: 耗时 44,685,858.1 ns, 故 GFOPS 为 2 / 0.0446858581 ≈ 44.76, 性能是Basic的 4,712,905,103.3 / 44,685,858.1 ≈ 105.47 倍.
  • BlockM4Nv3On512_ikj_M4Parallel: 耗时 4,842,750.9 ns, 故 GFOPS 为 2 / 0.0048427509 ≈ 412.99, 性能是Basic的 4,712,905,103.3 / 4,842,750.9 ≈ 973.19 倍.
  • BlockM4Nv3On512_ikj_M32: 耗时 32,179,510.0 ns, 故 GFOPS 为 2 / 0.0321795100 ≈ 62.15, 性能是Basic的 4,712,905,103.3 / 32,179,510.0 ≈ 146.46 倍.
  • BlockM4Nv3On512_ikj_M32Parallel: 耗时 4,363,112.2 ns, 故 GFOPS 为 2 / 0.0043631122 ≈ 458.39, 性能是Basic的 4,712,905,103.3 / 4,363,112.2 ≈ 1080.17 倍.
  • UseMathNet: 耗时 53,205,723.8 ns, 故 GFOPS 为 2 / 0.0532057238 ≈ 37.59, 性能是Basic的 4,712,905,103.3 / 53,205,723.8 ≈ 88.58 倍.
  • UseOpenBLAS: 耗时 2,932,070.9 ns, 故 GFOPS 为 2 / 0.0029320709 ≈ 682.11, 性能是Basic的 4,712,905,103.3 / 2,932,070.9 ≈ 1607.36 倍.
  • UseMKL: 耗时 4,379,927.2 ns, 故 GFOPS 为 2 / 0.0043799272 ≈ 456.63, 性能是Basic的 4,712,905,103.3 / 4,379,927.2 ≈ 1076.02 倍.
此时性能最好的并行算法是BlockM4Nv3On512_ikj_M32Parallel，它的性能是Basic的 1080.17 倍。而 MKL是 1076.02 倍，即该算法一般比MKL还稍微快一点。该算法比起 OpenBLAS ，差距在1倍以内，可看作同一梯队。
此时性能最好的单线程算法是BlockM4Nv3On512_ikj_M32，它的性能是Basic的 146.46 倍，超过了具有并行计算加速的UseMathNet。这说明算法很重要，否则并行计算加速的程序性能，有可能会被单线程程序超过。
五、基准测试结果

现代 .NET的跨平台功能非常成熟，能跨平台使用SIMD硬件加速。本章进行全面的基准测试，分别测试 X86、Arm处理器的性能。且不仅会测试Single类型，还会测试 Double类型的性能。
5.1 X86 架构 - .NET 9.0

当前CPU是 “AMD Ryzen 7 7840H”，操作系统是 “Windows 11”。
5.1.1 Single（单精度浮点型）

X86架构上、.NET 9.0 程序、Single类型的基准测试结果如下。
  1. BenchmarkDotNet v0.15.2, Windows 11 (10.0.26100.4946/24H2/2024Update/HudsonValley) AMD Ryzen 7 7840H w/ Radeon 780M Graphics 3.80GHz, 1 CPU, 16 logical and 8 physical cores .NET SDK 9.0.301 [Host] : .NET 9.0.6 (9.0.625.26613), X64 AOT AVX-512F+CD+BW+DQ+VL+VBMI MediumRun : .NET 9.0.6 (9.0.625.26613), X64 RyuJIT AVX-512F+CD+BW+DQ+VL+VBMI
复制代码
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,712,905,103.3 ns35,792,824.68 ns53,573,019.05 ns4,722,569,850.0 ns1.0000.02954 B
TileRow1024780,626,575.0 ns5,970,550.33 ns8,562,785.59 ns778,917,250.0 ns0.1660.001,238 B
TileRowSpan1024613,236,262.1 ns6,326,567.15 ns9,273,400.86 ns610,494,600.0 ns0.1300.001,358 B
TileRowRef1024332,063,973.2 ns4,375,391.05 ns6,275,055.62 ns329,953,500.0 ns0.0700.001,512 B
TileRowSimd102471,736,249.5 ns2,768,453.89 ns3,880,986.03 ns69,775,500.0 ns0.0150.001,258 B
TileRowSimdFmaBcl102471,320,997.3 ns1,540,394.38 ns2,056,382.42 ns70,551,614.3 ns0.0150.001,262 B
TileRowSimdFma102472,957,369.2 ns3,088,652.97 ns4,329,860.46 ns71,072,300.0 ns0.0150.001,262 B
TileRowSimdParallel102411,576,096.3 ns447,403.26 ns669,652.19 ns11,727,672.3 ns0.0020.008,549 B
TileRowSimdFmaParallel102412,234,503.5 ns337,579.48 ns505,273.12 ns12,306,358.6 ns0.0030.008,565 B
BlockM4Nv1_ijk_M321024155,479,981.2 ns925,593.09 ns1,327,458.06 ns155,197,087.5 ns0.0330.00NA
BlockM4Nv1_ijk_M32Parallel102418,389,778.0 ns391,156.50 ns573,352.18 ns18,200,718.8 ns0.0040.00NA
BlockM4Nv1_ikj_M321024158,203,231.7 ns1,098,377.31 ns1,643,999.58 ns157,511,900.0 ns0.0340.002,303 B
BlockM4Nv1_ikj_M32Parallel102415,417,779.4 ns528,972.18 ns775,360.62 ns15,101,640.6 ns0.0030.009,453 B
BlockM4Nv1_ikj_M41024159,175,464.8 ns2,913,579.61 ns4,084,432.04 ns161,698,000.0 ns0.0340.002,188 B
BlockM4Nv1_ikj_M4Parallel102417,944,453.9 ns584,658.22 ns856,984.51 ns17,830,587.5 ns0.0040.009,311 B
BlockM4Nv3_ikj_M4102452,756,227.9 ns1,133,938.16 ns1,626,260.36 ns53,519,345.5 ns0.0110.003,766 B
BlockM4Nv3_ikj_M4Parallel10246,234,993.0 ns450,007.94 ns645,388.00 ns5,976,738.7 ns0.0010.0010,419 B
BlockM4Nv3_ikj_M32102451,007,023.1 ns1,010,156.07 ns1,448,735.77 ns51,688,581.8 ns0.0110.003,819 B
BlockM4Nv3_ikj_M32Parallel10245,851,851.6 ns266,588.14 ns390,761.47 ns5,628,710.9 ns0.0010.0010,563 B
BlockM4Nv3On512_ikj_M4102444,685,858.1 ns4,745,584.17 ns6,652,646.79 ns41,372,192.9 ns0.0090.003,835 B
BlockM4Nv3On512_ikj_M4Parallel10244,842,750.9 ns306,368.06 ns429,485.26 ns4,669,524.2 ns0.0010.0010,499 B
BlockM4Nv3On512_ikj_M32102432,179,510.0 ns206,091.41 ns275,126.14 ns32,157,606.2 ns0.0070.003,666 B
BlockM4Nv3On512_ikj_M32Parallel10244,363,112.2 ns73,021.25 ns97,481.28 ns4,339,956.2 ns0.0010.0010,848 B
UseMathNet102453,205,723.8 ns836,527.21 ns1,226,170.84 ns52,803,610.0 ns0.0110.007,575 B
UseOpenBLAS10242,932,070.9 ns117,359.68 ns160,643.26 ns2,864,316.2 ns0.0010.00NA
UseMKL10244,379,927.2 ns58,880.02 ns88,128.84 ns4,340,348.8 ns0.0010.00508 B
测试结果对比如下。
Method 耗时 GFOPS 性能倍数
Basic4,712,905,103.3 ns0.421.00
TileRow780,626,575.0 ns2.566.04
TileRowSpan613,236,262.1 ns3.267.69
TileRowRef332,063,973.2 ns6.0214.19
TileRowSimd71,736,249.5 ns27.8865.70
TileRowSimdFmaBcl71,320,997.3 ns28.0466.08
TileRowSimdFma72,957,369.2 ns27.4164.60
TileRowSimdParallel11,576,096.3 ns172.77407.12
TileRowSimdFmaParallel12,234,503.5 ns163.47385.21
BlockM4Nv1_ijk_M32155,479,981.2 ns12.8630.31
BlockM4Nv1_ijk_M32Parallel18,389,778.0 ns108.76256.28
BlockM4Nv1_ikj_M32158,203,231.7 ns12.6429.79
BlockM4Nv1_ikj_M32Parallel15,417,779.4 ns129.72305.68
BlockM4Nv1_ikj_M4159,175,464.8 ns12.5629.61
BlockM4Nv1_ikj_M4Parallel17,944,453.9 ns111.46262.64
BlockM4Nv3_ikj_M452,756,227.9 ns37.9189.33
BlockM4Nv3_ikj_M4Parallel6,234,993.0 ns320.77755.88
BlockM4Nv3_ikj_M3251,007,023.1 ns39.2192.40
BlockM4Nv3_ikj_M32Parallel5,851,851.6 ns341.77805.37
BlockM4Nv3On512_ikj_M444,685,858.1 ns44.76105.47
BlockM4Nv3On512_ikj_M4Parallel4,842,750.9 ns412.99973.19
BlockM4Nv3On512_ikj_M3232,179,510.0 ns62.15146.46
BlockM4Nv3On512_ikj_M32Parallel4,363,112.2 ns458.391080.17
UseMathNet53,205,723.8 ns37.5988.58
UseOpenBLAS2,932,070.9 ns682.111607.36
UseMKL4,379,927.2 ns456.631076.02
我们表现最好的算法是 BlockM4Nv3On512_ikj_M32Parallel，它的性能是基础算法的 1080.17 倍（约为 1080），此时的每秒浮点计算次数是 458.39 GFOPS。
它的性能比 MKL 稍微快了一些。4,379,927.2 / 4,363,112.2 ≈ 1.003854，即快了 0.3854%。
且它与 OpenBLAS的差距在一倍以内，属于同一梯队。（2,932,070.9 / 4,363,112.2 ≈ 0.672013，即它能达到 OpenBLAS 67.2012% 的性能 ）
5.1.2 Double（双精度浮点型）

通过巧妙地使用 using指令，本程序还能测试DGEMM（Double General Matrix Multiply，双精度浮点数通用矩阵乘法）的性能。只是为了减少基准测试的时间，本程序默认关闭了它的基准测试。
若想测试 DGEMM，可以修改“MatrixNMultiplyBenchmark_Double.cs”，将第一行的 //#undef BENCHMARKS_OFF 去掉注释（//），开启DGEMM的测试。并修改“MatrixNMultiplyBenchmark_Single.cs”，将第一行的 #undef BENCHMARKS_OFF 的前面加注释（//），关闭SGEMM的测试。
X86架构上、.NET 9.0 程序、Double类型的基准测试结果如下。
  1. BenchmarkDotNet v0.15.2, Windows 11 (10.0.26100.4946/24H2/2024Update/HudsonValley) AMD Ryzen 7 7840H w/ Radeon 780M Graphics 3.80GHz, 1 CPU, 16 logical and 8 physical cores .NET SDK 9.0.301 [Host] : .NET 9.0.6 (9.0.625.26613), X64 AOT AVX-512F+CD+BW+DQ+VL+VBMI MediumRun : .NET 9.0.6 (9.0.625.26613), X64 RyuJIT AVX-512F+CD+BW+DQ+VL+VBMI
复制代码
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,667,707,243.3 ns101,856,765.08 ns152,454,422.51 ns4,722,324,950.0 ns1.0010.05955 B
TileRow1024785,539,882.8 ns6,951,632.14 ns10,189,613.10 ns781,816,000.0 ns0.1680.011,238 B
TileRowSpan1024612,766,200.0 ns5,098,750.48 ns7,312,476.18 ns611,489,700.0 ns0.1310.001,358 B
TileRowRef1024331,644,776.0 ns7,293,618.64 ns9,736,772.22 ns329,782,100.0 ns0.0710.001,632 B
TileRowSimd1024175,592,405.9 ns27,847,032.37 ns39,037,653.46 ns163,639,500.0 ns0.0380.01946 B
TileRowSimdFmaBcl1024165,678,173.1 ns5,407,093.39 ns7,401,291.09 ns163,372,412.5 ns0.0360.00946 B
TileRowSimdFma1024164,456,937.7 ns5,192,151.95 ns6,931,374.30 ns162,801,575.0 ns0.0350.00946 B
TileRowSimdParallel102427,533,111.2 ns351,339.32 ns525,868.19 ns27,299,954.7 ns0.0060.008,146 B
TileRowSimdFmaParallel102427,163,961.1 ns474,347.60 ns709,981.21 ns26,996,106.2 ns0.0060.008,279 B
BlockM4Nv1_ijk_M321024314,843,932.1 ns13,838,266.66 ns19,846,430.18 ns329,845,200.0 ns0.0680.00NA
BlockM4Nv1_ijk_M32Parallel102434,141,308.7 ns578,416.74 ns865,747.01 ns34,121,580.0 ns0.0070.00NA
BlockM4Nv1_ikj_M321024308,316,596.2 ns1,245,009.33 ns1,704,182.97 ns308,077,750.0 ns0.0660.00NA
BlockM4Nv1_ikj_M32Parallel102429,421,452.7 ns563,599.63 ns843,569.46 ns29,088,087.5 ns0.0060.009,392 B
BlockM4Nv1_ikj_M41024332,366,150.0 ns1,675,601.19 ns2,293,582.02 ns332,275,700.0 ns0.0710.00NA
BlockM4Nv1_ikj_M4Parallel102435,847,960.5 ns352,075.27 ns493,560.39 ns35,929,240.0 ns0.0080.009,268 B
BlockM4Nv3_ikj_M4102497,701,496.6 ns4,138,544.95 ns5,801,662.51 ns98,073,216.7 ns0.0210.003,830 B
BlockM4Nv3_ikj_M4Parallel102411,446,186.9 ns446,374.98 ns640,177.71 ns11,182,593.0 ns0.0020.0010,587 B
BlockM4Nv3_ikj_M32102491,633,206.1 ns1,885,472.12 ns2,704,088.00 ns91,787,078.6 ns0.0200.003,857 B
BlockM4Nv3_ikj_M32Parallel102410,177,112.8 ns351,523.74 ns492,787.22 ns9,954,175.0 ns0.0020.0010,737 B
BlockM4Nv3On512_ikj_M41024103,019,217.0 ns1,234,728.32 ns1,648,326.98 ns102,917,933.3 ns0.0220.003,831 B
BlockM4Nv3On512_ikj_M4Parallel10249,567,307.4 ns314,506.73 ns430,500.40 ns9,411,418.8 ns0.0020.0010,654 B
BlockM4Nv3On512_ikj_M32102466,571,730.6 ns722,460.47 ns1,058,973.28 ns66,147,812.5 ns0.0140.003,689 B
BlockM4Nv3On512_ikj_M32Parallel10249,291,286.4 ns267,082.18 ns365,585.13 ns9,083,106.2 ns0.0020.0010,906 B
UseMathNet102457,825,146.5 ns823,123.70 ns1,153,904.57 ns57,571,922.2 ns0.0120.007,554 B
UseOpenBLAS10245,834,910.0 ns419,034.34 ns600,966.58 ns5,468,444.9 ns0.0010.00NA
UseMKL102413,417,242.3 ns585,101.67 ns857,634.51 ns13,104,798.4 ns0.0030.00508 B
测试结果对比如下。
Method 耗时 GFOPS 性能倍数
Basic4,667,707,243.3 ns0.431.00
TileRow785,539,882.8 ns2.555.94
TileRowSpan612,766,200.0 ns3.267.62
TileRowRef331,644,776.0 ns6.0314.07
TileRowSimd175,592,405.9 ns11.3926.58
TileRowSimdFmaBcl165,678,173.1 ns12.0728.17
TileRowSimdFma164,456,937.7 ns12.1628.38
TileRowSimdParallel27,533,111.2 ns72.64169.53
TileRowSimdFmaParallel27,163,961.1 ns73.63171.83
BlockM4Nv1_ijk_M32314,843,932.1 ns6.3514.83
BlockM4Nv1_ijk_M32Parallel34,141,308.7 ns58.58136.72
BlockM4Nv1_ikj_M32308,316,596.2 ns6.4915.14
BlockM4Nv1_ikj_M32Parallel29,421,452.7 ns67.98158.65
BlockM4Nv1_ikj_M4332,366,150.0 ns6.0214.04
BlockM4Nv1_ikj_M4Parallel35,847,960.5 ns55.79130.21
BlockM4Nv3_ikj_M497,701,496.6 ns20.4747.78
BlockM4Nv3_ikj_M4Parallel11,446,186.9 ns174.73407.80
BlockM4Nv3_ikj_M3291,633,206.1 ns21.8350.94
BlockM4Nv3_ikj_M32Parallel10,177,112.8 ns196.52458.65
BlockM4Nv3On512_ikj_M4103,019,217.0 ns19.4145.31
BlockM4Nv3On512_ikj_M4Parallel9,567,307.4 ns209.05487.88
BlockM4Nv3On512_ikj_M3266,571,730.6 ns30.0470.12
BlockM4Nv3On512_ikj_M32Parallel9,291,286.4 ns215.26502.37
UseMathNet57,825,146.5 ns34.5980.72
UseOpenBLAS5,834,910.0 ns342.76799.96
UseMKL13,417,242.3 ns149.06347.89
我们表现最好的算法是 BlockM4Nv3On512_ikj_M32Parallel，它的性能是基础算法的 502.37 倍，此时的每秒浮点计算次数是 215.26 GFOPS。
它的性能比 MKL 强了不少。13,417,242.3 / 9,291,286.4 ≈ 1.444067，即快了 44.4067%。
且它与 OpenBLAS的差距在一倍以内，属于同一梯队。（5,834,910.0 / 9,291,286.4 ≈ 0.627998，即它能达到 OpenBLAS 62.7998% 的性能 ）
BlockM4Nv3On512_ikj_M32Parallel算法对于 Single 类型是458.39 GFOPS，Double类型是 215.26 GFOPS。Single的性能是Double的2倍左右（458.39 / 215.26 ≈ 2.1295），这与理论值相同（512位向量能同时计算16个Single，或8个Double）。
5.2 X86 架构 - .NET Framework

从 .NET Framework 4.5开始，开始提供了基本的SIMD硬件加速功能。本程序作了一些兼容处理，使得大部分算法，不用修改源代码，便能在 .NET Framework 上运行。现在来做一下基准测试。
由于 MKL、OpenBLAS 都是使用自己的动态库，与 .NET 版本无关。故在测试.NET Framework时，没必要测试 MKL、OpenBLAS。
当前CPU是 “AMD Ryzen 7 7840H”，操作系统是 “Windows 11”。
5.2.1 Single（单精度浮点型）

X86架构上、.NET Framework 程序、Single类型的基准测试结果如下。
  1. BenchmarkDotNet v0.15.2, Windows 11 (10.0.26100.4946/24H2/2024Update/HudsonValley) AMD Ryzen 7 7840H w/ Radeon 780M Graphics 3.80GHz, 1 CPU, 16 logical and 8 physical cores [Host] : .NET Framework 4.8.1 (4.8.9310.0), X64 RyuJIT VectorSize=256 MediumRun : .NET Framework 4.8.1 (4.8.9310.0), X64 RyuJIT VectorSize=256
复制代码
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10244,454,702.840 us96,843.2650 us144,950.4510 us4,487,360.900 us1.0010.055,687 B
TileRow1024966,327.050 us11,170.0099 us16,019.6958 us962,294.050 us0.2170.011,145 B
TileRowSpan10241,661,727.444 us16,190.1043 us22,696.2670 us1,657,810.200 us0.3730.011,793 B
TileRowRef1024313,248.250 us4,727.8263 us6,929.9871 us310,946.100 us0.0700.00512 B
TileRowSimd102483,757.697 us5,095.7616 us7,308.1896 us81,437.767 us0.0190.00736 B
TileRowSimdFma102490,579.311 us1,557.4982 us2,183.3951 us90,123.900 us0.0200.00746 B
TileRowSimdParallel102411,748.538 us192.8809 us288.6951 us11,706.321 us0.0030.007,799 B
TileRowSimdFmaParallel102412,838.369 us227.1818 us340.0350 us12,874.963 us0.0030.007,809 B
BlockM4Nv1_ijk_M32102479,924.653 us523.0694 us733.2703 us79,925.500 us0.0180.00NA
BlockM4Nv1_ijk_M32Parallel102410,929.294 us317.3551 us475.0022 us10,924.080 us0.0020.00NA
BlockM4Nv1_ikj_M32102479,275.543 us393.1614 us563.8604 us79,307.186 us0.0180.00NA
BlockM4Nv1_ikj_M32Parallel10249,784.844 us304.0108 us436.0032 us9,644.648 us0.0020.00NA
BlockM4Nv1_ikj_M4102479,669.252 us665.6703 us975.7309 us79,380.671 us0.0180.00NA
BlockM4Nv1_ikj_M4Parallel102411,542.263 us193.4875 us283.6115 us11,452.883 us0.0030.00NA
BlockM4Nv3_ikj_M4102441,454.059 us544.8342 us745.7753 us41,140.408 us0.0090.00NA
BlockM4Nv3_ikj_M4Parallel10247,906.391 us427.4508 us639.7883 us7,604.690 us0.0020.00NA
BlockM4Nv3_ikj_M32102441,720.557 us659.2914 us924.2345 us41,472.992 us0.0090.00NA
BlockM4Nv3_ikj_M32Parallel10246,674.869 us290.9010 us407.8025 us6,521.323 us0.0020.00NA
UseMathNet102454,745.115 us783.4180 us1,172.5833 us54,466.025 us0.0120.00279 B
测试结果对比如下。
Method 耗时 GFOPS 性能倍数
Basic4,454,702.840 us0.451.00
TileRow966,327.050 us2.074.61
TileRowSpan1,661,727.444 us1.202.68
TileRowRef313,248.250 us6.3814.22
TileRowSimd83,757.697 us23.8853.19
TileRowSimdFma90,579.311 us22.0849.18
TileRowSimdParallel11,748.538 us170.23379.17
TileRowSimdFmaParallel12,838.369 us155.78346.98
BlockM4Nv1_ijk_M3279,924.653 us25.0255.74
BlockM4Nv1_ijk_M32Parallel10,929.294 us182.99407.59
BlockM4Nv1_ikj_M3279,275.543 us25.2356.19
BlockM4Nv1_ikj_M32Parallel9,784.844 us204.40455.27
BlockM4Nv1_ikj_M479,669.252 us25.1055.91
BlockM4Nv1_ikj_M4Parallel11,542.263 us173.28385.95
BlockM4Nv3_ikj_M441,454.059 us48.25107.46
BlockM4Nv3_ikj_M4Parallel7,906.391 us252.96563.43
BlockM4Nv3_ikj_M3241,720.557 us47.94106.77
BlockM4Nv3_ikj_M32Parallel6,674.869 us299.63667.38
UseMathNet54,745.115 us36.5381.37
我们表现最好的算法是 BlockM4Nv3_ikj_M32Parallel，它的性能是基础算法的 667.38 倍，此时的每秒浮点计算次数是 299.63 GFOPS。而 .NET 9.0 是 458.39 GFOPS。.NET Framework 下能达到 65%左右的性能（299.63 / 458.39 ≈ 0.653657），表现已经很不错了。
还可以注意到，TileRowSpan的性能比较慢。这是因为 .NET Framework 下的Span为了保障兼容性，只能采取性能慢的保守做法。改为托管指针后，能避免这一瓶颈。
5.2.2 Double（双精度浮点型）

X86架构上、.NET Framework 程序、Double类型的基准测试结果如下。
  1. BenchmarkDotNet v0.15.2, Windows 11 (10.0.26100.4946/24H2/2024Update/HudsonValley) AMD Ryzen 7 7840H w/ Radeon 780M Graphics 3.80GHz, 1 CPU, 16 logical and 8 physical cores [Host] : .NET Framework 4.8.1 (4.8.9310.0), X64 RyuJIT VectorSize=256 MediumRun : .NET Framework 4.8.1 (4.8.9310.0), X64 RyuJIT VectorSize=256
复制代码
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10245,219,216.603 us555,505.4296 us814,252.7230 us4,593,251.700 us1.0240.225,687 B
TileRow1024976,825.167 us8,050.6241 us11,285.8515 us973,636.700 us0.1920.031,145 B
TileRowSpan10241,710,756.415 us15,104.7271 us20,675.5227 us1,712,805.450 us0.3360.051,793 B
TileRowRef1024329,519.486 us4,073.9472 us5,842.7338 us326,634.675 us0.0650.01512 B
TileRowSimd1024226,359.507 us8,571.1670 us12,292.5126 us221,183.317 us0.0440.01736 B
TileRowSimdFma1024233,129.988 us9,004.9541 us12,623.6891 us230,386.867 us0.0460.01746 B
TileRowSimdParallel102428,560.033 us251.5022 us376.4366 us28,418.494 us0.0060.007,799 B
TileRowSimdFmaParallel102429,554.618 us340.4174 us509.5208 us29,640.494 us0.0060.007,809 B
BlockM4Nv1_ijk_M321024162,907.984 us570.6867 us781.1624 us162,883.600 us0.0320.00NA
BlockM4Nv1_ijk_M32Parallel102420,194.872 us437.9025 us628.0267 us19,967.569 us0.0040.00NA
BlockM4Nv1_ikj_M321024162,750.319 us1,913.2804 us2,618.9201 us162,468.392 us0.0320.00NA
BlockM4Nv1_ikj_M32Parallel102419,956.415 us356.5837 us533.7177 us19,780.269 us0.0040.00NA
BlockM4Nv1_ikj_M41024162,956.796 us783.5110 us1,018.7856 us163,084.538 us0.0320.00NA
BlockM4Nv1_ikj_M4Parallel102422,891.121 us381.4015 us570.8639 us22,858.144 us0.0040.00NA
BlockM4Nv3_ikj_M4102483,848.352 us3,621.1132 us5,419.9123 us80,793.929 us0.0160.00NA
BlockM4Nv3_ikj_M4Parallel102415,405.993 us592.3353 us886.5797 us15,284.463 us0.0030.00NA
BlockM4Nv3_ikj_M32102480,101.363 us2,405.1717 us3,449.4257 us79,586.886 us0.0160.00NA
BlockM4Nv3_ikj_M32Parallel102413,063.492 us528.6187 us758.1291 us12,659.439 us0.0030.00NA
UseMathNet102462,791.007 us962.8950 us1,411.3992 us62,336.978 us0.0120.00279 B
测试结果对比如下。
Method 耗时 GFOPS 性能倍数
Basic5,219,216.603 us0.381.00
TileRow976,825.167 us2.055.34
TileRowSpan1,710,756.415 us1.173.05
TileRowRef329,519.486 us6.0715.84
TileRowSimd226,359.507 us8.8423.06
TileRowSimdFma233,129.988 us8.5822.39
TileRowSimdParallel28,560.033 us70.03182.75
TileRowSimdFmaParallel29,554.618 us67.67176.60
BlockM4Nv1_ijk_M32162,907.984 us12.2832.04
BlockM4Nv1_ijk_M32Parallel20,194.872 us99.04258.44
BlockM4Nv1_ikj_M32162,750.319 us12.2932.07
BlockM4Nv1_ikj_M32Parallel19,956.415 us100.22261.53
BlockM4Nv1_ikj_M4162,956.796 us12.2732.03
BlockM4Nv1_ikj_M4Parallel22,891.121 us87.37228.00
BlockM4Nv3_ikj_M483,848.352 us23.8562.25
BlockM4Nv3_ikj_M4Parallel15,405.993 us129.82338.78
BlockM4Nv3_ikj_M3280,101.363 us24.9765.16
BlockM4Nv3_ikj_M32Parallel13,063.492 us153.10399.53
UseMathNet62,791.007 us31.8583.12
我们表现最好的算法是 BlockM4Nv3_ikj_M32Parallel，它的性能是基础算法的 399.53 倍，此时的每秒浮点计算次数是 153.10 GFOPS。而 .NET 9.0 是 215.26 GFOPS。.NET Framework 下能达到 71%左右的性能（153.10 / 215.26 ≈ 0.711233），表现已经很不错了。
5.3 Arm 架构 - .NET 9.0

同样的源代码，可以在Arm架构的处理器上运行。
由于 MKL 不支持 Arm架构，且目前nuget上没有Arm架构的OpenBLAS包，故没有测试它们。
当前CPU是 “Apple M2”，操作系统是 “macOS Sequoia 15.6.1”。
5.3.1 Single（单精度浮点型）

Arm架构上、.NET 9.0 程序、Single类型的基准测试结果如下。
  1. BenchmarkDotNet v0.15.2, macOS Sequoia 15.6.1 (24G90) [Darwin 24.6.0] Apple M2, 1 CPU, 8 logical and 8 physical cores .NET SDK 9.0.102 [Host] : .NET 9.0.1 (9.0.124.61010), Arm64 AOT AdvSIMD MediumRun : .NET 9.0.1 (9.0.124.61010), Arm64 RyuJIT AdvSIMD
复制代码
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10243,754,237.106 us2,074.9078 us3,041.3732 us3,753,734.041 us1.0000.00
TileRow1024840,254.040 us142.7877 us200.1684 us840,254.500 us0.2240.00
TileRowSpan1024761,125.180 us498.8317 us715.4096 us760,974.583 us0.2030.00
TileRowRef1024330,818.469 us106.5126 us152.7572 us330,825.885 us0.0880.00
TileRowSimd1024117,413.868 us90.4758 us129.7577 us117,410.512 us0.0310.00
TileRowSimdFmaBcl1024117,210.060 us67.9110 us99.5431 us117,214.725 us0.0310.00
TileRowSimdFma1024117,189.597 us64.3478 us94.3201 us117,187.067 us0.0310.00
TileRowSimdParallel102426,510.611 us185.3270 us277.3887 us26,510.653 us0.0070.00
TileRowSimdFmaParallel102425,593.386 us211.9040 us317.1679 us25,553.712 us0.0070.00
BlockM4Nv1_ijk_M321024152,420.390 us107.9145 us161.5213 us152,420.901 us0.0410.00
BlockM4Nv1_ijk_M32Parallel102438,359.750 us243.9897 us342.0395 us38,300.182 us0.0100.00
BlockM4Nv1_ikj_M321024148,065.516 us119.2124 us178.4315 us148,031.271 us0.0390.00
BlockM4Nv1_ikj_M32Parallel102436,589.371 us249.5473 us349.8305 us36,578.652 us0.0100.00
BlockM4Nv1_ikj_M41024146,542.856 us441.5088 us647.1579 us146,317.281 us0.0390.00
BlockM4Nv1_ikj_M4Parallel102436,993.328 us268.5206 us393.5940 us37,065.863 us0.0100.00
BlockM4Nv3_ikj_M41024489,578.227 us2,903.1298 us4,069.7828 us493,334.646 us0.1300.00
BlockM4Nv3_ikj_M4Parallel102417,373.172 us80.5072 us115.4610 us17,399.023 us0.0050.00
BlockM4Nv3_ikj_M321024487,247.653 us1,268.0786 us1,777.6693 us488,772.104 us0.1300.00
BlockM4Nv3_ikj_M32Parallel102417,908.937 us119.7722 us179.2693 us17,917.776 us0.0050.00
BlockM4Nv3On512_ikj_M41024NANANANA??
BlockM4Nv3On512_ikj_M4Parallel1024NANANANA??
BlockM4Nv3On512_ikj_M321024NANANANA??
BlockM4Nv3On512_ikj_M32Parallel1024NANANANA??
UseMathNet1024154,141.940 us837.3170 us1,253.2568 us153,909.807 us0.0410.00
测试结果对比如下。
Method 耗时 GFOPS 性能倍数
Basic3,754,237.106 us0.531.00
TileRow840,254.040 us2.384.47
TileRowSpan761,125.180 us2.634.93
TileRowRef330,818.469 us6.0511.35
TileRowSimd117,413.868 us17.0331.97
TileRowSimdFmaBcl117,210.060 us17.0632.03
TileRowSimdFma117,189.597 us17.0732.04
TileRowSimdParallel26,510.611 us75.44141.61
TileRowSimdFmaParallel25,593.386 us78.15146.69
BlockM4Nv1_ijk_M32152,420.390 us13.1224.63
BlockM4Nv1_ijk_M32Parallel38,359.750 us52.1497.87
BlockM4Nv1_ikj_M32148,065.516 us13.5125.36
BlockM4Nv1_ikj_M32Parallel36,589.371 us54.66102.60
BlockM4Nv1_ikj_M4146,542.856 us13.6525.62
BlockM4Nv1_ikj_M4Parallel36,993.328 us54.06101.48
BlockM4Nv3_ikj_M4489,578.227 us4.097.67
BlockM4Nv3_ikj_M4Parallel17,373.172 us115.12216.09
BlockM4Nv3_ikj_M32487,247.653 us4.107.70
BlockM4Nv3_ikj_M32Parallel17,908.937 us111.68209.63
UseMathNet154,141.940 us12.9824.36
我们表现最好的算法是 BlockM4Nv3_ikj_M4Parallel，它的性能是基础算法的 216.09 倍，此时的每秒浮点计算次数是 115.12 GFOPS。
AMD Ryzen 7 7840H 的Single运算速度是458.39 GFOPS。故对于单精度浮点数（Single）的运算，现代CPU均能达到 每秒千亿次 （100GFOPS）级别的运算速度。
Apple M2的向量宽度是128位，而AMD Ryzen 7 7840H 的向量宽度是 512位。假设 Apple M2的向量宽度也达到512位的话，那么理论上会有4倍性能提升，即理论上会是 115.12*4=460.48 GFOPS，这与 AMD Ryzen 7 7840H 的成绩非常接近。
5.3.2 Double（双精度浮点型）

Arm架构上、.NET 9.0 程序、Double类型的基准测试结果如下。
  1. BenchmarkDotNet v0.15.2, macOS Sequoia 15.6.1 (24G90) [Darwin 24.6.0] Apple M2, 1 CPU, 8 logical and 8 physical cores .NET SDK 9.0.102 [Host] : .NET 9.0.1 (9.0.124.61010), Arm64 AOT AdvSIMD MediumRun : .NET 9.0.1 (9.0.124.61010), Arm64 RyuJIT AdvSIMD
复制代码
Method N Mean Error StdDev Median Ratio RatioSD Code Size
Basic10243,979,217.591 us3,084.0733 us4,423.0862 us3,979,746.459 us1.0000.00
TileRow1024880,127.965 us571.4544 us855.3261 us880,017.875 us0.2210.00
TileRowSpan1024797,048.246 us486.5361 us713.1583 us797,011.417 us0.2000.00
TileRowRef1024348,452.077 us324.7616 us465.7634 us348,371.812 us0.0880.00
TileRowSimd1024242,866.758 us231.8556 us325.0292 us242,941.903 us0.0610.00
TileRowSimdFmaBcl1024241,862.247 us306.1996 us429.2491 us241,910.333 us0.0610.00
TileRowSimdFma1024241,918.592 us195.5033 us292.6202 us241,910.924 us0.0610.00
TileRowSimdParallel102457,335.473 us813.2521 us1,192.0544 us57,027.893 us0.0140.00
TileRowSimdFmaParallel102455,068.792 us505.7084 us741.2608 us54,973.000 us0.0140.00
BlockM4Nv1_ijk_M321024268,039.331 us388.0223 us580.7736 us267,957.021 us0.0670.00
BlockM4Nv1_ijk_M32Parallel102465,948.754 us894.1736 us1,338.3570 us66,392.143 us0.0170.00
BlockM4Nv1_ikj_M321024265,015.397 us231.9062 us332.5930 us265,013.260 us0.0670.00
BlockM4Nv1_ikj_M32Parallel102465,237.383 us958.7143 us1,434.9585 us65,071.180 us0.0160.00
BlockM4Nv1_ikj_M41024257,808.595 us340.6194 us488.5062 us257,729.927 us0.0650.00
BlockM4Nv1_ikj_M4Parallel102465,421.381 us691.5592 us1,013.6786 us65,436.662 us0.0160.00
BlockM4Nv3_ikj_M41024594,005.257 us137,852.8179 us206,331.6239 us694,300.630 us0.1490.05
BlockM4Nv3_ikj_M4Parallel102435,414.032 us584.9320 us875.4988 us35,398.204 us0.0090.00
BlockM4Nv3_ikj_M321024587,200.067 us134,524.2102 us201,349.5203 us685,058.547 us0.1480.05
BlockM4Nv3_ikj_M32Parallel102435,606.024 us439.7341 us630.6535 us35,467.028 us0.0090.00
BlockM4Nv3On512_ikj_M41024NANANANA??
BlockM4Nv3On512_ikj_M4Parallel1024NANANANA??
BlockM4Nv3On512_ikj_M321024NANANANA??
BlockM4Nv3On512_ikj_M32Parallel1024NANANANA??
UseMathNet1024172,491.316 us1,034.3648 us1,415.8503 us172,643.514 us0.0430.00
测试结果对比如下。
Method 耗时 GFOPS 性能倍数
Basic3,979,217.591 us0.501.00
TileRow880,127.965 us2.274.52
TileRowSpan797,048.246 us2.514.99
TileRowRef348,452.077 us5.7411.42
TileRowSimd242,866.758 us8.2316.38
TileRowSimdFmaBcl241,862.247 us8.2716.45
TileRowSimdFma241,918.592 us8.2716.45
TileRowSimdParallel57,335.473 us34.8869.40
TileRowSimdFmaParallel55,068.792 us36.3272.26
BlockM4Nv1_ijk_M32268,039.331 us7.4614.85
BlockM4Nv1_ijk_M32Parallel65,948.754 us30.3360.34
BlockM4Nv1_ikj_M32265,015.397 us7.5515.02
BlockM4Nv1_ikj_M32Parallel65,237.383 us30.6661.00
BlockM4Nv1_ikj_M4257,808.595 us7.7615.43
BlockM4Nv1_ikj_M4Parallel65,421.381 us30.5760.82
BlockM4Nv3_ikj_M4594,005.257 us3.376.70
BlockM4Nv3_ikj_M4Parallel35,414.032 us56.47112.36
BlockM4Nv3_ikj_M32587,200.067 us3.416.78
BlockM4Nv3_ikj_M32Parallel35,606.024 us56.17111.76
UseMathNet172,491.316 us11.5923.07
我们表现最好的算法是 BlockM4Nv3_ikj_M4Parallel，它的性能是基础算法的 112.36 倍，此时的每秒浮点计算次数是 56.47 GFOPS。
该算法对于 Single 类型是115.12 GFOPS，Double类型是 56.47 GFOPS。Single的性能是Double的2倍左右（115.12 / 56.47 ≈ 2.0386），这与理论值相同。
六、结语

我最开始对 C# 的矩阵乘法程序进行优化时，没料到能取得这么大的进展。前几个版本的算法，与采用手动汇编优化等优化技术的MKL、OpenBLAS比起来，性能差距望尘莫及。当时觉得我优化的算法，能超过老牌的 MathNet就很不错了。
经过不断改进，性能大幅度提升，最终取得了惊人的成果——不仅单线程成绩就能超过并行计算的MathNet，而且并行计算的成绩能小幅超过MKL，甚至与OpenBLAS的差距在一倍以内，处于同一梯队。
即纯 C# 编写的托管代码程序，与采用了手动汇编优化等优化技术的专业矩阵运算库，达到同一梯队的性能水平。
首先，这得益于 .NET 9.0 生成的汇编代码的质量比较高。例如 .NET Framework 生成的汇编代码就差一些，性能只有 .NET 9.0 的 60%左右。
其次，虽然.NET 9.0 生成的汇编代码质量比较高，但比起MKL、OpenBLAS里手工优化的汇编代码，各方面还是差了不少。但得益于现代 CPU的解码、乱序执行、分支预测等功能的日益强大，即使面对质量稍差汇编代码，CPU也会尽可能高效地运算。
其实矩阵乘法程序还有很大的优化空间，例如 使用各种CPU架构的矩阵运算加速专用指令、为各种矩阵尺寸找出最佳的分块尺寸等。现在OpenBLAS在很多场合超过 MKL，就是因为这些原因。
希望我的文章能给读者们带来启发。
附录


  • 源代码地址: https://github.com/zyl910/MatrixBenchmarkCs/tree/v0.1
  • 【深缓中字】为什么6层嵌套循环让这个算法快了120倍？
  • 矩阵乘法优化过程（DGEMM）
  • 关于sgemm_hsw的一点解释说明
  • 通用矩阵乘法（GEMM）原理实现与优化方法总览
  • OpenBLAS项目与矩阵乘法优化 | AI 研习社
  • OpenBLAS
  • Intel® oneAPI Math Kernel Library (oneMKL)
  • C# 使用SIMD向量类型加速浮点数组求和运算（1）：使用Vector4、Vector
  • C# 使用SIMD向量类型加速浮点数组求和运算（4）：用引用代替指针, 摆脱unsafe关键字，兼谈Unsafe类的使用
出处：http://www.cnblogs.com/zyl910/ 版权声明：自由转载-非商用-非衍生-保持署名 | Creative Commons BY-NC-ND 3.0.

来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！
使用NET跨平台SIMD硬件

使用道具 举报

相关推荐
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册
发帖
陆菊
4 天前

0

粉丝关注

26

主题发布

板块介绍填写区域，请于后台编辑
微信扫一扫