[C#] 使用 .NET 的跨平台SIMD硬件加速功能，将 GEMM（通用矩阵乘法）算法速度提升1080倍，比肩 MKL、OpenBLAS

[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 列的点积。

MultiplyMatrix_basic

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)，那么它每秒浮点运算次数为——

每秒浮点运算次数 = 浮点运算总次数 / 耗时 = 2 / t (GFOPS)

1.2 C++ 开发的矩阵乘法最基础实现

根据上面的矩阵乘法定义，可以编程实现它。例如用 C++ 开发的矩阵乘法最基础实现如下。

void gemmVanilla(const std::vector<float> &matA, const std::vector<float> &matB,
  std::vector<float> &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<float> 类型来传递矩阵数据。虽然 std::vector<float> 本身仅能存储一维数据，而这段代码隐式约定了其存储了 matSize * matSize 个数据，正好是二维的矩阵数据。由于 std::vector<float> 是线性存储数据的，所以此时存储的矩阵，可看作行主序的矩阵。

这段代码的开头做了“清除矩阵C”（Clear matC）的工作。随后按照矩阵乘法定义，编写3层循环，实现了矩阵乘法的功能。

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

1.3.1 矩阵乘法的实现

将上述代码由 C++ 转为 C# 语言，便能得到 C# 版矩阵乘法。见下面的 StaticBasic 方法。

// My type.
using TMy = Single;

partial class MultiplyMatrixStatic {

    /// <summary>
    /// Basic on Array.
    /// </summary>
    /// <param name="M">The number of rows in matrix A (矩阵A的行数).</param>
    /// <param name="N">The number of columns in matrix B (矩阵B的列数).</param>
    /// <param name="K">The number of columns in matrix A, or the number of rows in matrix B (矩阵A的列数, 或矩阵B的行数).</param>
    /// <param name="A">Matrix A.</param>
    /// <param name="strideA">Stride of A.</param>
    /// <param name="B">Matrix B.</param>
    /// <param name="strideB">Stride of B.</param>
    /// <param name="C">Matrix C.</param>
    /// <param name="strideC">Stride of C.</param>
    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方法，它的基准测试代码如下。

[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 里编写了条件分支来引入这些包。

<PropertyGroup Condition=" '$(ExUseNativeDll)' == '' AND '$(TargetFramework)' == 'net9.0' AND '$(OS)' == 'Windows_NT' ">
    <ExUseNativeDll>true</ExUseNativeDll>
</PropertyGroup>
<PropertyGroup Condition="'$(ExUseNativeDll)' == ''">
    <ExUseNativeDll>false</ExUseNativeDll>
</PropertyGroup>
<PropertyGroup Condition="'$(ExUseNativeDll)' == 'true'">
    <DefineConstants>$(DefineConstants);USE_NATIVE_DLL</DefineConstants>
</PropertyGroup>
<ItemGroup Condition="'$(ExUseNativeDll)' == 'true'">
    <PackageReference Include="MKL.NET" Version="$(MKLNETVersion)" />
    <PackageReference Include="OpenBlasSharp" Version="$(OpenBlasSharpVersion)" />
</ItemGroup>
<Choose>
    <When Condition=" '$(ExUseNativeDll)' == 'true' AND '$(OS)' == 'Windows_NT' ">
        <ItemGroup>
            <PackageReference Include="MKL.NET.win-x64" Version="$(MKLNETwinx64Version)" />
            <PackageReference Include="OpenBlasSharp.Windows" Version="$(OpenBlasSharpWindowsVersion)" />
            <!--
            <PackageReference Include="MKL.NET.win-86" Version="$(MKLNETwinx64Version)" />
            -->
        </ItemGroup>
    </When>
    <When Condition=" '$(ExUseNativeDll)' == 'true' AND '$([System.Runtime.InteropServices.RuntimeInformation]::IsOSPlatform($([System.Runtime.InteropServices.OSPlatform]::Linux)))' ">
        <ItemGroup>
            <PackageReference Include="MKL.NET.linux-x64" Version="$(MKLNETlinuxx64Version)" />
        </ItemGroup>
    </When>
    <When Condition=" '$(ExUseNativeDll)' == 'true' AND '$([System.Runtime.InteropServices.RuntimeInformation]::IsOSPlatform($([System.Runtime.InteropServices.OSPlatform]::OSX)))' ">
        <ItemGroup>
            <PackageReference Include="MKL.NET.osx-x64" Version="$(MKLNETosxx64Version)" />
        </ItemGroup>
    </When>
    <Otherwise>
    </Otherwise>
</Choose>

上面配置为——仅 .NET 9.0、Windows操作系统时，才使用原生动态库。若你想在其他平台上测试，可以修改 <ExUseNativeDll>false</ExUseNativeDll> 中的值。将它改为“true”。

版本号是在 Directory.Build.props 里定义的。

<MathNetNumericsVersion>6.0.0-beta2</MathNetNumericsVersion>

<MKLNETVersion>1.6.0</MKLNETVersion>
<MKLNETwinx64Version>2022.0.0.115</MKLNETwinx64Version>
<MKLNETlinuxx64Version>2022.0.1.117</MKLNETlinuxx64Version>
<MKLNETosxx64Version>2022.0.0.105</MKLNETosxx64Version>

<OpenBlasSharpVersion>0.3.4</OpenBlasSharpVersion>
<OpenBlasSharpWindowsVersion>2025.1.13</OpenBlasSharpWindowsVersion>

1.4.2 这些库的基准测试方法

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

由于 MathNet 的矩阵运算函数不直接支持数组，而需要用它提供的矩阵类型。于是按它要求，定义了相关变量，并进行了初始化。

protected MathNet.Numerics.LinearAlgebra.Matrix<TMy>? matA;
protected MathNet.Numerics.LinearAlgebra.Matrix<TMy>? matB;
protected MathNet.Numerics.LinearAlgebra.Matrix<TMy>? matC;

protected override void GlobalSetupCore() {
    base.GlobalSetupCore();
    try {
        matA = MathNet.Numerics.LinearAlgebra.Matrix<TMy>.Build.DenseOfRowMajor(MatrixM, MatrixK, arrayA);
        matB = MathNet.Numerics.LinearAlgebra.Matrix<TMy>.Build.DenseOfRowMajor(MatrixK, MatrixN, arrayB);
        matC = MathNet.Numerics.LinearAlgebra.Matrix<TMy>.Build.Dense(MatrixM, MatrixN);
    } catch (Exception ex) {
        BenchmarkUtil.WriteItem("# Setup-error", string.Format("{0}", ex.Message));
    }
}

MathNet的基准测试的方法如下。

[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的基准测试的方法如下。

[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的基准测试的方法如下。

[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，目前的基准测试结果如下。

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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
UseMathNet	1024	53,205,723.8 ns	836,527.21 ns	1,226,170.84 ns	52,803,610.0 ns	0.011	0.00	7,575 B
UseOpenBLAS	1024	2,932,070.9 ns	117,359.68 ns	160,643.26 ns	2,864,316.2 ns	0.001	0.00	NA
UseMKL	1024	4,379,927.2 ns	58,880.02 ns	88,128.84 ns	4,340,348.8 ns	0.001	0.00	508 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 算法的实现

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”来实现该功能，源代码如下。

/// <summary>
/// Fill value (填充值).
/// </summary>
/// <typeparam name="T">The element type (元素的类型).</typeparam>
/// <param name="value">The value (值).</param>
/// <param name="rows">The number of rows in matrix (矩阵的行数).</param>
/// <param name="cols">The number of columns in matrix (矩阵的列数).</param>
/// <param name="matrix">The matrix (矩阵).</param>
/// <param name="stride">The stride of matrix. When it is 0, use cols (矩阵的跨距. 为 0 时 使用 cols).</param>
/// <param name="start">The start index of matrix (矩阵的开始索引).</param>
public static void Fill<T>(T value, int rows, int cols, Span<T> matrix, int stride = 0, int start = 0) {
    if (0 == stride) {
        stride = cols;
    }
    int idx0 = start;
    for (int i = 0; i < rows; i++) {
        Span<T> span = matrix.Slice(idx0, cols);
        span.Fill(value);
        // Next.
        idx0 += stride;
    }
}

这个工具函数，是逐行对数组进行处理的。对于每一行的内容，转为Span做清零处理。当一行处理完毕时，会根据stride参数移动到下一行的起始索引。

2.1.2 基准测试方法

对于上面的方法，它的基准测试代码如下。

[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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
TileRow	1024	780,626,575.0 ns	5,970,550.33 ns	8,562,785.59 ns	778,917,250.0 ns	0.166	0.00	1,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时矩阵乘法的内循环。

multiplymatrix4_ijk

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时矩阵乘法的内循环。

multiplymatrix4_ikj

于是可以直观的看到——绿色的长方体在B矩阵、C矩阵上，均仅跨越了虚线，表示它们都是连续内存访问。

对比后可以发现，ikj比ijk次序在内存访问方面具有优势，这就是它性能更好的原因。

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

观察上一个算法，会发现在内循环里，除了矩阵乘法运算所需的乘法，还用到2次乘法——用乘法来计算索引（bIdx、cIdx）。也就是说上一个算法至少有 3 * N^3次乘法运算，而标准矩阵乘法是 N^3次乘法运算，它多了 2倍的 N^3次乘法运算。

消除索引计算时的乘法开销，能够提高性能。具体办法是根据循环的特点，在各级循环里做合适的累加操作，便能算出各索引的值。

另外，还可以用Span代替数组，使该方法不仅能支持数据等托管数据，还能支持非托管的数据。

2.3.1 算法的实现

public static void StaticTileRowSpan(int M, int N, int K, Span<TMy> A, int strideA, Span<TMy> B, int strideB, Span<TMy> 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 基准测试方法

对于上面的方法，它的基准测试代码如下。

[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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
TileRowSpan	1024	613,236,262.1 ns	6,326,567.15 ns	9,273,400.86 ns	610,494,600.0 ns	0.130	0.00	1,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 算法的实现

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<TMy*>(A)。

2.4.1.1 矩阵清零

因现在使用了托管指针，这就没法使用先前的Span版参数的“MatrixUtil.Fill”方法。于是新增了一个托管指针版参数的 “MatrixUtil.Fill”方法。

public static void Fill<T>(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<T> 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 基准测试方法

对于上面的方法，它的基准测试代码如下。

[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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
TileRowRef	1024	332,063,973.2 ns	4,375,391.05 ns	6,275,055.62 ns	329,953,500.0 ns	0.070	0.00	1,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<TMy>.Count 的值，那么SIMD向量算法的关键思路如下。

“A的索引保持不变”：使用 Vector<TMy> vA = new Vector<TMy>(pA)，就能根据当前A矩阵中的元素，广播到向量变量里（vA[i] = pA）。
“B、C的索引仅需+1”：例如使用 ref Vector<TMy> pB = ref Unsafe.As<TMy, Vector<TMy>>(ref pB0)，就能直接将 TMy 的引用，重新解释为 Vector<TMy>的引用。继而能够从 pB0 这个地址开始，依次将“W个连续值”加载到向量变量里（(pB[i]) = pB0[i]）。同理可这样加载pC。

随后执行 pC = Vector.Add(Vector.Multiply(vA, pB), pC)，便能对这“W个连续值”中的每个元素执行 pC += vA * pB 的操作。

标量算法每次仅能处理1个值，而SIMD向量算法每次能处理 W个值，故它的理论计算性能是标量算法的 W倍。

若读者对 .NET SIMD硬件加速技术不熟悉，可查阅 “C# 使用SIMD向量类型加速浮点数组求和运算（1）：使用Vector4、Vector<T>” 等文章。

multiplymatrix4_ikj

假设现在的向量类型是128位（如X86的SSE指令集、Arm的AdvSimd指令集），那么向量类型的字节尺寸为 128 / 8 = 16。因Single类型的尺寸是4个字节，那么向量类型里可以存储 16 / 4 = 4个元素。即 W（Vector<TMy>.Count）的值为4，很适合处理N（B、C矩阵的列数）为4时的矩阵。例如上图，它就是4时的矩阵，绿色长方体在标量算法时需要写循环来实现，而用了SIMD向量算法后能一次性算出这4个值。

当N为W的整数倍时，可以多次执行上述的SIMD向量算法，使整行的数据处理完毕。若不是整数倍，则需专门处理尾部数据，下面的源代码里使用“尾部覆盖”的策略。

3.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<TMy>.Count || !Vector.IsHardwareAccelerated) {
        StaticTileRowRef(M, N, K, in A, strideA, in B, strideB, ref C, strideC);
        return;
    }
    int cntRem = N % Vector<TMy>.Count; // Remainder count.
    int cntBlockRaw = N / Vector<TMy>.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<TMy> vA = new Vector<TMy>(aValue);
            // Last.
            int pos = N - Vector<TMy>.Count;
            ref Vector<TMy> pBLast = ref Unsafe.As<TMy, Vector<TMy>>(ref Unsafe.Add(ref pB0, pos));
            ref Vector<TMy> pCLast = ref Unsafe.As<TMy, Vector<TMy>>(ref Unsafe.Add(ref pC0, pos));
            Vector<TMy> vCLast = Vector.Add(Vector.Multiply(vA, pBLast), pCLast);
            // SIMD for.
            if (cntBlock >= 0) {
                ref Vector<TMy> pB = ref Unsafe.As<TMy, Vector<TMy>>(ref pB0);
                ref Vector<TMy> pC = ref Unsafe.As<TMy, Vector<TMy>>(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<TMy>.Count（N小于向量宽度）或没有硬件加速时，便进行回退，调用标量算法 StaticTileRowRef。

注意 pB、pC 都是 Vector<TMy> 类型的托管指针，所以在 Unsafe.Add(ref pB, 1)时，并不是对单个 TMy 元素移动指针（增加 sizeof(TMy)字节），而是对整块 Vector<TMy> 移动指针（增加 sizeof(Vector<TMy>)字节）。

3.1.2 基准测试方法

对于上面的方法，它的基准测试代码如下。

[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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
TileRowSimd	1024	71,736,249.5 ns	2,768,453.89 ns	3,880,986.03 ns	69,775,500.0 ns	0.015	0.00	1,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算法，仅需修改基准测试方法，便能增加它的并行加速的基准测试方法。

[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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
TileRowSimdParallel	1024	11,576,096.3 ns	447,403.26 ns	669,652.19 ns	11,727,672.3 ns	0.002	0.00	8,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 算法的实现

#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<TMy>.Count || !Vector.IsHardwareAccelerated) {
        StaticTileRowRef(M, N, K, in A, strideA, in B, strideB, ref C, strideC);
        return;
    }
    int cntRem = N % Vector<TMy>.Count; // Remainder count.
    int cntBlockRaw = N / Vector<TMy>.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<TMy> vA = new Vector<TMy>(aValue);
            // Last.
            int pos = N - Vector<TMy>.Count;
            ref Vector<TMy> pBLast = ref Unsafe.As<TMy, Vector<TMy>>(ref Unsafe.Add(ref pB0, pos));
            ref Vector<TMy> pCLast = ref Unsafe.As<TMy, Vector<TMy>>(ref Unsafe.Add(ref pC0, pos));
            Vector<TMy> vCLast = Vector.FusedMultiplyAdd(vA, pBLast, pCLast);
            // SIMD for.
            if (cntBlock >= 0) {
                ref Vector<TMy> pB = ref Unsafe.As<TMy, Vector<TMy>>(ref pB0);
                ref Vector<TMy> pC = ref Unsafe.As<TMy, Vector<TMy>>(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 基准测试方法

对于上面的方法，它的基准测试代码如下。

#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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
TileRowSimd	1024	71,736,249.5 ns	2,768,453.89 ns	3,880,986.03 ns	69,775,500.0 ns	0.015	0.00	1,258 B
TileRowSimdFmaBcl	1024	71,320,997.3 ns	1,540,394.38 ns	2,056,382.42 ns	70,551,614.3 ns	0.015	0.00	1,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 指令，使其不再受到编译质量的影响。源代码如下。

#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<TMy>.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<TMy>.Count; // Remainder count.
    int cntBlockRaw = N / Vector<TMy>.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<TMy> vA = new Vector<TMy>(aValue);
            // Last.
            int pos = N - Vector<TMy>.Count;
            ref Vector<TMy> pBLast = ref Unsafe.As<TMy, Vector<TMy>>(ref Unsafe.Add(ref pB0, pos));
            ref Vector<TMy> pCLast = ref Unsafe.As<TMy, Vector<TMy>>(ref Unsafe.Add(ref pC0, pos));
            Vector<TMy> vCLast = Vector.Add(Vector.Multiply(vA, pBLast), pCLast);
            // SIMD for.
            if (cntBlock >= 0) {
                ref Vector<TMy> pB = ref Unsafe.As<TMy, Vector<TMy>>(ref pB0);
                ref Vector<TMy> pC = ref Unsafe.As<TMy, Vector<TMy>>(ref pC0);
                for (int j = 0; j < cntBlock; ++j) {
                    // pC += vA * pB;
                    if (Vector<byte>.Count == Vector256<byte>.Count) {
                        pC = Fma.MultiplyAdd(vA.AsVector256(), pB.AsVector256(), pC.AsVector256()).AsVector();
                    } else if (Vector<byte>.Count == Vector256<byte>.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”方法。

[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static Vector<float> MultiplyAdd(Vector<float> left, Vector<float> right, Vector<float> 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一样简化矩阵乘法的代码。

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<TMy>.Count || !Vector.IsHardwareAccelerated) {
        StaticTileRowRef(M, N, K, in A, strideA, in B, strideB, ref C, strideC);
        return;
    }
    int cntRem = N % Vector<TMy>.Count; // Remainder count.
    int cntBlockRaw = N / Vector<TMy>.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<TMy> vA = new Vector<TMy>(aValue);
            // Last.
            int pos = N - Vector<TMy>.Count;
            ref Vector<TMy> pBLast = ref Unsafe.As<TMy, Vector<TMy>>(ref Unsafe.Add(ref pB0, pos));
            ref Vector<TMy> pCLast = ref Unsafe.As<TMy, Vector<TMy>>(ref Unsafe.Add(ref pC0, pos));
            //Vector<TMy> vCLast = Vector.FusedMultiplyAdd(vA, pBLast, pCLast);
            Vector<TMy> vCLast = VectorHelper.MultiplyAdd(vA, pBLast, pCLast);
            // SIMD for.
            if (cntBlock >= 0) {
                ref Vector<TMy> pB = ref Unsafe.As<TMy, Vector<TMy>>(ref pB0);
                ref Vector<TMy> pC = ref Unsafe.As<TMy, Vector<TMy>>(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 基准测试方法

对于上面的方法，它的基准测试代码如下。

[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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
TileRowSimd	1024	71,736,249.5 ns	2,768,453.89 ns	3,880,986.03 ns	69,775,500.0 ns	0.015	0.00	1,258 B
TileRowSimdFmaBcl	1024	71,320,997.3 ns	1,540,394.38 ns	2,056,382.42 ns	70,551,614.3 ns	0.015	0.00	1,262 B
TileRowSimdFma	1024	72,957,369.2 ns	3,088,652.97 ns	4,329,860.46 ns	71,072,300.0 ns	0.015	0.00	1,262 B
TileRowSimdParallel	1024	11,576,096.3 ns	447,403.26 ns	669,652.19 ns	11,727,672.3 ns	0.002	0.00	8,549 B
TileRowSimdFmaParallel	1024	12,234,503.5 ns	337,579.48 ns	505,273.12 ns	12,306,358.6 ns	0.003	0.00	8,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）)

multiplymatrix4_ikj

设向量宽度 W 为 Vector<TMy>.Count 的值，此时向量算法内循环是处理了 1行、W个列的数据。

若向量算法不再仅处理1行数据，而是处理4行数据，便能进一步提高性能。将它称作“4行、1倍向量宽度”算法。

multiplymatrix4_m4n4

上面的立方图图片展示了 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 类中增加了这些函数

#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<T> LoadUnsafe<T>(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<T, Vector<T>>(ref Unsafe.AsRef(in source));
#endif // NET8_0_OR_GREATER
}

[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void StoreUnsafe<T>(
#if !NET8_0_OR_GREATER
    this
#endif // NET8_0_OR_GREATER
    Vector<T> 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<T, Vector<T>>(ref destination) = source;
#endif // NET8_0_OR_GREATER
}

从 .NET 8.0 开始，向量类型取消了 struct 的约束。于是本类也遵从了这个规则。

4.1.2 矩阵乘法内核算法

将矩阵进行分块后，对于每个分块，专门用一个方法来处比较好。这种方法，不再时处理整个矩阵，而仅是单个分块，故一般称它为内核算法。

内核算法不用管理矩阵的清零操作，且一般仅需支持固定尺寸的数据。

为了便于测试各种矩阵分块策略，我这边会传递固定尺寸的整数倍的数据。于是m、n、k这3个参数也需做循环处理，只是不用考虑非整数倍数据。

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<TMy>.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<TMy> c0;
            Vector<TMy> c1;
            Vector<TMy> c2;
            Vector<TMy> 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<TMy> b = VectorHelper.LoadUnsafe(ref pB);
                ref TMy pA = ref pALine;
                c0 = VectorHelper.MultiplyAdd(new Vector<TMy>(pA), b, c0); pA = ref Unsafe.Add(ref pA, strideA);
                c1 = VectorHelper.MultiplyAdd(new Vector<TMy>(pA), b, c1); pA = ref Unsafe.Add(ref pA, strideA);
                c2 = VectorHelper.MultiplyAdd(new Vector<TMy>(pA), b, c2); pA = ref Unsafe.Add(ref pA, strideA);
                c3 = VectorHelper.MultiplyAdd(new Vector<TMy>(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 参数。

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<TMy>((void*)addressA);
                ref TMy refB = ref Unsafe.AsRef<TMy>((void*)addressB);
                ref TMy refC = ref Unsafe.AsRef<TMy>((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 基准测试方法

对于上面的方法，它的基准测试代码如下。

[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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
BlockM4Nv1_ijk_M32	1024	155,479,981.2 ns	925,593.09 ns	1,327,458.06 ns	155,197,087.5 ns	0.033	0.00	NA
BlockM4Nv1_ijk_M32Parallel	1024	18,389,778.0 ns	391,156.50 ns	573,352.18 ns	18,200,718.8 ns	0.004	0.00	NA

对测试结果进行对比——

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 顺序。

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<TMy>((void*)addressA);
                ref TMy refB = ref Unsafe.AsRef<TMy>((void*)addressB);
                ref TMy refC = ref Unsafe.AsRef<TMy>((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<TMy>.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<TMy>.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 基准测试方法

对于上面的方法，它的基准测试代码如下。

[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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
BlockM4Nv1_ikj_M32	1024	158,203,231.7 ns	1,098,377.31 ns	1,643,999.58 ns	157,511,900.0 ns	0.034	0.00	2,303 B
BlockM4Nv1_ikj_M32Parallel	1024	15,417,779.4 ns	528,972.18 ns	775,360.62 ns	15,101,640.6 ns	0.003	0.00	9,453 B
BlockM4Nv1_ikj_M4	1024	159,175,464.8 ns	2,913,579.61 ns	4,084,432.04 ns	161,698,000.0 ns	0.034	0.00	2,188 B
BlockM4Nv1_ikj_M4Parallel	1024	17,944,453.9 ns	584,658.22 ns	856,984.51 ns	17,830,587.5 ns	0.004	0.00	9,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 增加了这些公共方法。

[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void LoadX3Unsafe<T>(ref readonly T source, out Vector<T> data0, out Vector<T> data1, out Vector<T> data2)
#if VECTOR_WHERE_STRUCT
        where T : struct
#endif // VECTOR_WHERE_STRUCT
        {
    ref Vector<T> p = ref Unsafe.As<T, Vector<T>>(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<T>(ref T destination, Vector<T> data0, Vector<T> data1, Vector<T> data2)
#if VECTOR_WHERE_STRUCT
        where T : struct
#endif // VECTOR_WHERE_STRUCT
        {
    ref Vector<T> p = ref Unsafe.As<T, Vector<T>>(ref destination);
    p = data0;
    Unsafe.Add(ref p, 1) = data1;
    Unsafe.Add(ref p, 2) = data2;
}

4.3.2 矩阵乘法内核算法

随后改进内核算法，让它一次性处理3倍向量宽度。

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<TMy>.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<TMy> c0_0, c0_1, c0_2;
            Vector<TMy> c1_0, c1_1, c1_2;
            Vector<TMy> c2_0, c2_1, c2_2;
            Vector<TMy> 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<TMy> b0, b1, b2;
                Vector<TMy> va;
                VectorHelper.LoadX3Unsafe(ref pB, out b0, out b1, out b2);
                ref TMy pA = ref pALine;
                va = new Vector<TMy>(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<TMy>(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<TMy>(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<TMy>(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循环，且实现计算好了各自的循环次数。

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<TMy>.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<TMy>((void*)addressA);
                ref TMy refB = ref Unsafe.AsRef<TMy>((void*)addressB);
                ref TMy refC = ref Unsafe.AsRef<TMy>((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<TMy>.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<TMy>.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 基准测试方法

对于上面的方法，它的基准测试代码如下。

[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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
BlockM4Nv3_ikj_M4	1024	52,756,227.9 ns	1,133,938.16 ns	1,626,260.36 ns	53,519,345.5 ns	0.011	0.00	3,766 B
BlockM4Nv3_ikj_M4Parallel	1024	6,234,993.0 ns	450,007.94 ns	645,388.00 ns	5,976,738.7 ns	0.001	0.00	10,419 B
BlockM4Nv3_ikj_M32	1024	51,007,023.1 ns	1,010,156.07 ns	1,448,735.77 ns	51,688,581.8 ns	0.011	0.00	3,819 B
BlockM4Nv3_ikj_M32Parallel	1024	5,851,851.6 ns	266,588.14 ns	390,761.47 ns	5,628,710.9 ns	0.001	0.00	10,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指令集了。大多数情况时，使用 Vector256 静态类的方法，会更加方便——对于大多数代码，仅需将 Vector 替换为 Vector512。

4.4.1 公共方法的改进

Vector512 也需要 LoadX3Unsafe等公共方法。于是新增了Vector512Helper类，增加了这些公共方法。

[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static void LoadX3Unsafe<T>(ref readonly T source, out Vector512<T> data0, out Vector512<T> data1, out Vector512<T> data2)
#if VECTOR_WHERE_STRUCT
        where T : struct
#endif // VECTOR_WHERE_STRUCT
        {
    ref Vector512<T> p = ref Unsafe.As<T, Vector512<T>>(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<T>(ref T destination, Vector512<T> data0, Vector512<T> data1, Vector512<T> data2)
#if VECTOR_WHERE_STRUCT
        where T : struct
#endif // VECTOR_WHERE_STRUCT
        {
    ref Vector512<T> p = ref Unsafe.As<T, Vector512<T>>(ref destination);
    p = data0;
    Unsafe.Add(ref p, 1) = data1;
    Unsafe.Add(ref p, 2) = data2;
}

[MethodImpl(MethodImplOptions.AggressiveInlining)]
public static Vector512<float> MultiplyAdd(Vector512<float> left, Vector512<float> right, Vector512<float> 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。

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<TMy>.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<TMy> c0_0, c0_1, c0_2;
            Vector512<TMy> c1_0, c1_1, c1_2;
            Vector512<TMy> c2_0, c2_1, c2_2;
            Vector512<TMy> 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<TMy> b0, b1, b2;
                Vector512<TMy> 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的幂次”的矩阵。

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<TMy>.IsSupported) {
        throw new NotSupportedException("Vector512.IsHardwareAccelerated is false!");
    }
    int vectorWidth = Vector<TMy>.Count;
    int vectorWidth3 = Vector512<TMy>.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<TMy>((void*)addressA);
                ref TMy refB = ref Unsafe.AsRef<TMy>((void*)addressB);
                ref TMy refC = ref Unsafe.AsRef<TMy>((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<TMy>.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<TMy>.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 基准测试方法

对于上面的方法，它的基准测试代码如下。

#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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
BlockM4Nv3On512_ikj_M4	1024	44,685,858.1 ns	4,745,584.17 ns	6,652,646.79 ns	41,372,192.9 ns	0.009	0.00	3,835 B
BlockM4Nv3On512_ikj_M4Parallel	1024	4,842,750.9 ns	306,368.06 ns	429,485.26 ns	4,669,524.2 ns	0.001	0.00	10,499 B
BlockM4Nv3On512_ikj_M32	1024	32,179,510.0 ns	206,091.41 ns	275,126.14 ns	32,157,606.2 ns	0.007	0.00	3,666 B
BlockM4Nv3On512_ikj_M32Parallel	1024	4,363,112.2 ns	73,021.25 ns	97,481.28 ns	4,339,956.2 ns	0.001	0.00	10,848 B
UseMathNet	1024	53,205,723.8 ns	836,527.21 ns	1,226,170.84 ns	52,803,610.0 ns	0.011	0.00	7,575 B
UseOpenBLAS	1024	2,932,070.9 ns	117,359.68 ns	160,643.26 ns	2,864,316.2 ns	0.001	0.00	NA
UseMKL	1024	4,379,927.2 ns	58,880.02 ns	88,128.84 ns	4,340,348.8 ns	0.001	0.00	508 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类型的基准测试结果如下。

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
Basic	1024	4,712,905,103.3 ns	35,792,824.68 ns	53,573,019.05 ns	4,722,569,850.0 ns	1.000	0.02	954 B
TileRow	1024	780,626,575.0 ns	5,970,550.33 ns	8,562,785.59 ns	778,917,250.0 ns	0.166	0.00	1,238 B
TileRowSpan	1024	613,236,262.1 ns	6,326,567.15 ns	9,273,400.86 ns	610,494,600.0 ns	0.130	0.00	1,358 B
TileRowRef	1024	332,063,973.2 ns	4,375,391.05 ns	6,275,055.62 ns	329,953,500.0 ns	0.070	0.00	1,512 B
TileRowSimd	1024	71,736,249.5 ns	2,768,453.89 ns	3,880,986.03 ns	69,775,500.0 ns	0.015	0.00	1,258 B
TileRowSimdFmaBcl	1024	71,320,997.3 ns	1,540,394.38 ns	2,056,382.42 ns	70,551,614.3 ns	0.015	0.00	1,262 B
TileRowSimdFma	1024	72,957,369.2 ns	3,088,652.97 ns	4,329,860.46 ns	71,072,300.0 ns	0.015	0.00	1,262 B
TileRowSimdParallel	1024	11,576,096.3 ns	447,403.26 ns	669,652.19 ns	11,727,672.3 ns	0.002	0.00	8,549 B
TileRowSimdFmaParallel	1024	12,234,503.5 ns	337,579.48 ns	505,273.12 ns	12,306,358.6 ns	0.003	0.00	8,565 B
BlockM4Nv1_ijk_M32	1024	155,479,981.2 ns	925,593.09 ns	1,327,458.06 ns	155,197,087.5 ns	0.033	0.00	NA
BlockM4Nv1_ijk_M32Parallel	1024	18,389,778.0 ns	391,156.50 ns	573,352.18 ns	18,200,718.8 ns	0.004	0.00	NA
BlockM4Nv1_ikj_M32	1024	158,203,231.7 ns	1,098,377.31 ns	1,643,999.58 ns	157,511,900.0 ns	0.034	0.00	2,303 B
BlockM4Nv1_ikj_M32Parallel	1024	15,417,779.4 ns	528,972.18 ns	775,360.62 ns	15,101,640.6 ns	0.003	0.00	9,453 B
BlockM4Nv1_ikj_M4	1024	159,175,464.8 ns	2,913,579.61 ns	4,084,432.04 ns	161,698,000.0 ns	0.034	0.00	2,188 B
BlockM4Nv1_ikj_M4Parallel	1024	17,944,453.9 ns	584,658.22 ns	856,984.51 ns	17,830,587.5 ns	0.004	0.00	9,311 B
BlockM4Nv3_ikj_M4	1024	52,756,227.9 ns	1,133,938.16 ns	1,626,260.36 ns	53,519,345.5 ns	0.011	0.00	3,766 B
BlockM4Nv3_ikj_M4Parallel	1024	6,234,993.0 ns	450,007.94 ns	645,388.00 ns	5,976,738.7 ns	0.001	0.00	10,419 B
BlockM4Nv3_ikj_M32	1024	51,007,023.1 ns	1,010,156.07 ns	1,448,735.77 ns	51,688,581.8 ns	0.011	0.00	3,819 B
BlockM4Nv3_ikj_M32Parallel	1024	5,851,851.6 ns	266,588.14 ns	390,761.47 ns	5,628,710.9 ns	0.001	0.00	10,563 B
BlockM4Nv3On512_ikj_M4	1024	44,685,858.1 ns	4,745,584.17 ns	6,652,646.79 ns	41,372,192.9 ns	0.009	0.00	3,835 B
BlockM4Nv3On512_ikj_M4Parallel	1024	4,842,750.9 ns	306,368.06 ns	429,485.26 ns	4,669,524.2 ns	0.001	0.00	10,499 B
BlockM4Nv3On512_ikj_M32	1024	32,179,510.0 ns	206,091.41 ns	275,126.14 ns	32,157,606.2 ns	0.007	0.00	3,666 B
BlockM4Nv3On512_ikj_M32Parallel	1024	4,363,112.2 ns	73,021.25 ns	97,481.28 ns	4,339,956.2 ns	0.001	0.00	10,848 B
UseMathNet	1024	53,205,723.8 ns	836,527.21 ns	1,226,170.84 ns	52,803,610.0 ns	0.011	0.00	7,575 B
UseOpenBLAS	1024	2,932,070.9 ns	117,359.68 ns	160,643.26 ns	2,864,316.2 ns	0.001	0.00	NA
UseMKL	1024	4,379,927.2 ns	58,880.02 ns	88,128.84 ns	4,340,348.8 ns	0.001	0.00	508 B

测试结果对比如下。

Method	耗时	GFOPS	性能倍数
Basic	4,712,905,103.3 ns	0.42	1.00
TileRow	780,626,575.0 ns	2.56	6.04
TileRowSpan	613,236,262.1 ns	3.26	7.69
TileRowRef	332,063,973.2 ns	6.02	14.19
TileRowSimd	71,736,249.5 ns	27.88	65.70
TileRowSimdFmaBcl	71,320,997.3 ns	28.04	66.08
TileRowSimdFma	72,957,369.2 ns	27.41	64.60
TileRowSimdParallel	11,576,096.3 ns	172.77	407.12
TileRowSimdFmaParallel	12,234,503.5 ns	163.47	385.21
BlockM4Nv1_ijk_M32	155,479,981.2 ns	12.86	30.31
BlockM4Nv1_ijk_M32Parallel	18,389,778.0 ns	108.76	256.28
BlockM4Nv1_ikj_M32	158,203,231.7 ns	12.64	29.79
BlockM4Nv1_ikj_M32Parallel	15,417,779.4 ns	129.72	305.68
BlockM4Nv1_ikj_M4	159,175,464.8 ns	12.56	29.61
BlockM4Nv1_ikj_M4Parallel	17,944,453.9 ns	111.46	262.64
BlockM4Nv3_ikj_M4	52,756,227.9 ns	37.91	89.33
BlockM4Nv3_ikj_M4Parallel	6,234,993.0 ns	320.77	755.88
BlockM4Nv3_ikj_M32	51,007,023.1 ns	39.21	92.40
BlockM4Nv3_ikj_M32Parallel	5,851,851.6 ns	341.77	805.37
BlockM4Nv3On512_ikj_M4	44,685,858.1 ns	44.76	105.47
BlockM4Nv3On512_ikj_M4Parallel	4,842,750.9 ns	412.99	973.19
BlockM4Nv3On512_ikj_M32	32,179,510.0 ns	62.15	146.46
BlockM4Nv3On512_ikj_M32Parallel	4,363,112.2 ns	458.39	1080.17
UseMathNet	53,205,723.8 ns	37.59	88.58
UseOpenBLAS	2,932,070.9 ns	682.11	1607.36
UseMKL	4,379,927.2 ns	456.63	1076.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类型的基准测试结果如下。

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
Basic	1024	4,667,707,243.3 ns	101,856,765.08 ns	152,454,422.51 ns	4,722,324,950.0 ns	1.001	0.05	955 B
TileRow	1024	785,539,882.8 ns	6,951,632.14 ns	10,189,613.10 ns	781,816,000.0 ns	0.168	0.01	1,238 B
TileRowSpan	1024	612,766,200.0 ns	5,098,750.48 ns	7,312,476.18 ns	611,489,700.0 ns	0.131	0.00	1,358 B
TileRowRef	1024	331,644,776.0 ns	7,293,618.64 ns	9,736,772.22 ns	329,782,100.0 ns	0.071	0.00	1,632 B
TileRowSimd	1024	175,592,405.9 ns	27,847,032.37 ns	39,037,653.46 ns	163,639,500.0 ns	0.038	0.01	946 B
TileRowSimdFmaBcl	1024	165,678,173.1 ns	5,407,093.39 ns	7,401,291.09 ns	163,372,412.5 ns	0.036	0.00	946 B
TileRowSimdFma	1024	164,456,937.7 ns	5,192,151.95 ns	6,931,374.30 ns	162,801,575.0 ns	0.035	0.00	946 B
TileRowSimdParallel	1024	27,533,111.2 ns	351,339.32 ns	525,868.19 ns	27,299,954.7 ns	0.006	0.00	8,146 B
TileRowSimdFmaParallel	1024	27,163,961.1 ns	474,347.60 ns	709,981.21 ns	26,996,106.2 ns	0.006	0.00	8,279 B
BlockM4Nv1_ijk_M32	1024	314,843,932.1 ns	13,838,266.66 ns	19,846,430.18 ns	329,845,200.0 ns	0.068	0.00	NA
BlockM4Nv1_ijk_M32Parallel	1024	34,141,308.7 ns	578,416.74 ns	865,747.01 ns	34,121,580.0 ns	0.007	0.00	NA
BlockM4Nv1_ikj_M32	1024	308,316,596.2 ns	1,245,009.33 ns	1,704,182.97 ns	308,077,750.0 ns	0.066	0.00	NA
BlockM4Nv1_ikj_M32Parallel	1024	29,421,452.7 ns	563,599.63 ns	843,569.46 ns	29,088,087.5 ns	0.006	0.00	9,392 B
BlockM4Nv1_ikj_M4	1024	332,366,150.0 ns	1,675,601.19 ns	2,293,582.02 ns	332,275,700.0 ns	0.071	0.00	NA
BlockM4Nv1_ikj_M4Parallel	1024	35,847,960.5 ns	352,075.27 ns	493,560.39 ns	35,929,240.0 ns	0.008	0.00	9,268 B
BlockM4Nv3_ikj_M4	1024	97,701,496.6 ns	4,138,544.95 ns	5,801,662.51 ns	98,073,216.7 ns	0.021	0.00	3,830 B
BlockM4Nv3_ikj_M4Parallel	1024	11,446,186.9 ns	446,374.98 ns	640,177.71 ns	11,182,593.0 ns	0.002	0.00	10,587 B
BlockM4Nv3_ikj_M32	1024	91,633,206.1 ns	1,885,472.12 ns	2,704,088.00 ns	91,787,078.6 ns	0.020	0.00	3,857 B
BlockM4Nv3_ikj_M32Parallel	1024	10,177,112.8 ns	351,523.74 ns	492,787.22 ns	9,954,175.0 ns	0.002	0.00	10,737 B
BlockM4Nv3On512_ikj_M4	1024	103,019,217.0 ns	1,234,728.32 ns	1,648,326.98 ns	102,917,933.3 ns	0.022	0.00	3,831 B
BlockM4Nv3On512_ikj_M4Parallel	1024	9,567,307.4 ns	314,506.73 ns	430,500.40 ns	9,411,418.8 ns	0.002	0.00	10,654 B
BlockM4Nv3On512_ikj_M32	1024	66,571,730.6 ns	722,460.47 ns	1,058,973.28 ns	66,147,812.5 ns	0.014	0.00	3,689 B
BlockM4Nv3On512_ikj_M32Parallel	1024	9,291,286.4 ns	267,082.18 ns	365,585.13 ns	9,083,106.2 ns	0.002	0.00	10,906 B
UseMathNet	1024	57,825,146.5 ns	823,123.70 ns	1,153,904.57 ns	57,571,922.2 ns	0.012	0.00	7,554 B
UseOpenBLAS	1024	5,834,910.0 ns	419,034.34 ns	600,966.58 ns	5,468,444.9 ns	0.001	0.00	NA
UseMKL	1024	13,417,242.3 ns	585,101.67 ns	857,634.51 ns	13,104,798.4 ns	0.003	0.00	508 B

测试结果对比如下。

Method	耗时	GFOPS	性能倍数
Basic	4,667,707,243.3 ns	0.43	1.00
TileRow	785,539,882.8 ns	2.55	5.94
TileRowSpan	612,766,200.0 ns	3.26	7.62
TileRowRef	331,644,776.0 ns	6.03	14.07
TileRowSimd	175,592,405.9 ns	11.39	26.58
TileRowSimdFmaBcl	165,678,173.1 ns	12.07	28.17
TileRowSimdFma	164,456,937.7 ns	12.16	28.38
TileRowSimdParallel	27,533,111.2 ns	72.64	169.53
TileRowSimdFmaParallel	27,163,961.1 ns	73.63	171.83
BlockM4Nv1_ijk_M32	314,843,932.1 ns	6.35	14.83
BlockM4Nv1_ijk_M32Parallel	34,141,308.7 ns	58.58	136.72
BlockM4Nv1_ikj_M32	308,316,596.2 ns	6.49	15.14
BlockM4Nv1_ikj_M32Parallel	29,421,452.7 ns	67.98	158.65
BlockM4Nv1_ikj_M4	332,366,150.0 ns	6.02	14.04
BlockM4Nv1_ikj_M4Parallel	35,847,960.5 ns	55.79	130.21
BlockM4Nv3_ikj_M4	97,701,496.6 ns	20.47	47.78
BlockM4Nv3_ikj_M4Parallel	11,446,186.9 ns	174.73	407.80
BlockM4Nv3_ikj_M32	91,633,206.1 ns	21.83	50.94
BlockM4Nv3_ikj_M32Parallel	10,177,112.8 ns	196.52	458.65
BlockM4Nv3On512_ikj_M4	103,019,217.0 ns	19.41	45.31
BlockM4Nv3On512_ikj_M4Parallel	9,567,307.4 ns	209.05	487.88
BlockM4Nv3On512_ikj_M32	66,571,730.6 ns	30.04	70.12
BlockM4Nv3On512_ikj_M32Parallel	9,291,286.4 ns	215.26	502.37
UseMathNet	57,825,146.5 ns	34.59	80.72
UseOpenBLAS	5,834,910.0 ns	342.76	799.96
UseMKL	13,417,242.3 ns	149.06	347.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类型的基准测试结果如下。

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
Basic	1024	4,454,702.840 us	96,843.2650 us	144,950.4510 us	4,487,360.900 us	1.001	0.05	5,687 B
TileRow	1024	966,327.050 us	11,170.0099 us	16,019.6958 us	962,294.050 us	0.217	0.01	1,145 B
TileRowSpan	1024	1,661,727.444 us	16,190.1043 us	22,696.2670 us	1,657,810.200 us	0.373	0.01	1,793 B
TileRowRef	1024	313,248.250 us	4,727.8263 us	6,929.9871 us	310,946.100 us	0.070	0.00	512 B
TileRowSimd	1024	83,757.697 us	5,095.7616 us	7,308.1896 us	81,437.767 us	0.019	0.00	736 B
TileRowSimdFma	1024	90,579.311 us	1,557.4982 us	2,183.3951 us	90,123.900 us	0.020	0.00	746 B
TileRowSimdParallel	1024	11,748.538 us	192.8809 us	288.6951 us	11,706.321 us	0.003	0.00	7,799 B
TileRowSimdFmaParallel	1024	12,838.369 us	227.1818 us	340.0350 us	12,874.963 us	0.003	0.00	7,809 B
BlockM4Nv1_ijk_M32	1024	79,924.653 us	523.0694 us	733.2703 us	79,925.500 us	0.018	0.00	NA
BlockM4Nv1_ijk_M32Parallel	1024	10,929.294 us	317.3551 us	475.0022 us	10,924.080 us	0.002	0.00	NA
BlockM4Nv1_ikj_M32	1024	79,275.543 us	393.1614 us	563.8604 us	79,307.186 us	0.018	0.00	NA
BlockM4Nv1_ikj_M32Parallel	1024	9,784.844 us	304.0108 us	436.0032 us	9,644.648 us	0.002	0.00	NA
BlockM4Nv1_ikj_M4	1024	79,669.252 us	665.6703 us	975.7309 us	79,380.671 us	0.018	0.00	NA
BlockM4Nv1_ikj_M4Parallel	1024	11,542.263 us	193.4875 us	283.6115 us	11,452.883 us	0.003	0.00	NA
BlockM4Nv3_ikj_M4	1024	41,454.059 us	544.8342 us	745.7753 us	41,140.408 us	0.009	0.00	NA
BlockM4Nv3_ikj_M4Parallel	1024	7,906.391 us	427.4508 us	639.7883 us	7,604.690 us	0.002	0.00	NA
BlockM4Nv3_ikj_M32	1024	41,720.557 us	659.2914 us	924.2345 us	41,472.992 us	0.009	0.00	NA
BlockM4Nv3_ikj_M32Parallel	1024	6,674.869 us	290.9010 us	407.8025 us	6,521.323 us	0.002	0.00	NA
UseMathNet	1024	54,745.115 us	783.4180 us	1,172.5833 us	54,466.025 us	0.012	0.00	279 B

测试结果对比如下。

Method	耗时	GFOPS	性能倍数
Basic	4,454,702.840 us	0.45	1.00
TileRow	966,327.050 us	2.07	4.61
TileRowSpan	1,661,727.444 us	1.20	2.68
TileRowRef	313,248.250 us	6.38	14.22
TileRowSimd	83,757.697 us	23.88	53.19
TileRowSimdFma	90,579.311 us	22.08	49.18
TileRowSimdParallel	11,748.538 us	170.23	379.17
TileRowSimdFmaParallel	12,838.369 us	155.78	346.98
BlockM4Nv1_ijk_M32	79,924.653 us	25.02	55.74
BlockM4Nv1_ijk_M32Parallel	10,929.294 us	182.99	407.59
BlockM4Nv1_ikj_M32	79,275.543 us	25.23	56.19
BlockM4Nv1_ikj_M32Parallel	9,784.844 us	204.40	455.27
BlockM4Nv1_ikj_M4	79,669.252 us	25.10	55.91
BlockM4Nv1_ikj_M4Parallel	11,542.263 us	173.28	385.95
BlockM4Nv3_ikj_M4	41,454.059 us	48.25	107.46
BlockM4Nv3_ikj_M4Parallel	7,906.391 us	252.96	563.43
BlockM4Nv3_ikj_M32	41,720.557 us	47.94	106.77
BlockM4Nv3_ikj_M32Parallel	6,674.869 us	299.63	667.38
UseMathNet	54,745.115 us	36.53	81.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类型的基准测试结果如下。

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
Basic	1024	5,219,216.603 us	555,505.4296 us	814,252.7230 us	4,593,251.700 us	1.024	0.22	5,687 B
TileRow	1024	976,825.167 us	8,050.6241 us	11,285.8515 us	973,636.700 us	0.192	0.03	1,145 B
TileRowSpan	1024	1,710,756.415 us	15,104.7271 us	20,675.5227 us	1,712,805.450 us	0.336	0.05	1,793 B
TileRowRef	1024	329,519.486 us	4,073.9472 us	5,842.7338 us	326,634.675 us	0.065	0.01	512 B
TileRowSimd	1024	226,359.507 us	8,571.1670 us	12,292.5126 us	221,183.317 us	0.044	0.01	736 B
TileRowSimdFma	1024	233,129.988 us	9,004.9541 us	12,623.6891 us	230,386.867 us	0.046	0.01	746 B
TileRowSimdParallel	1024	28,560.033 us	251.5022 us	376.4366 us	28,418.494 us	0.006	0.00	7,799 B
TileRowSimdFmaParallel	1024	29,554.618 us	340.4174 us	509.5208 us	29,640.494 us	0.006	0.00	7,809 B
BlockM4Nv1_ijk_M32	1024	162,907.984 us	570.6867 us	781.1624 us	162,883.600 us	0.032	0.00	NA
BlockM4Nv1_ijk_M32Parallel	1024	20,194.872 us	437.9025 us	628.0267 us	19,967.569 us	0.004	0.00	NA
BlockM4Nv1_ikj_M32	1024	162,750.319 us	1,913.2804 us	2,618.9201 us	162,468.392 us	0.032	0.00	NA
BlockM4Nv1_ikj_M32Parallel	1024	19,956.415 us	356.5837 us	533.7177 us	19,780.269 us	0.004	0.00	NA
BlockM4Nv1_ikj_M4	1024	162,956.796 us	783.5110 us	1,018.7856 us	163,084.538 us	0.032	0.00	NA
BlockM4Nv1_ikj_M4Parallel	1024	22,891.121 us	381.4015 us	570.8639 us	22,858.144 us	0.004	0.00	NA
BlockM4Nv3_ikj_M4	1024	83,848.352 us	3,621.1132 us	5,419.9123 us	80,793.929 us	0.016	0.00	NA
BlockM4Nv3_ikj_M4Parallel	1024	15,405.993 us	592.3353 us	886.5797 us	15,284.463 us	0.003	0.00	NA
BlockM4Nv3_ikj_M32	1024	80,101.363 us	2,405.1717 us	3,449.4257 us	79,586.886 us	0.016	0.00	NA
BlockM4Nv3_ikj_M32Parallel	1024	13,063.492 us	528.6187 us	758.1291 us	12,659.439 us	0.003	0.00	NA
UseMathNet	1024	62,791.007 us	962.8950 us	1,411.3992 us	62,336.978 us	0.012	0.00	279 B

测试结果对比如下。

Method	耗时	GFOPS	性能倍数
Basic	5,219,216.603 us	0.38	1.00
TileRow	976,825.167 us	2.05	5.34
TileRowSpan	1,710,756.415 us	1.17	3.05
TileRowRef	329,519.486 us	6.07	15.84
TileRowSimd	226,359.507 us	8.84	23.06
TileRowSimdFma	233,129.988 us	8.58	22.39
TileRowSimdParallel	28,560.033 us	70.03	182.75
TileRowSimdFmaParallel	29,554.618 us	67.67	176.60
BlockM4Nv1_ijk_M32	162,907.984 us	12.28	32.04
BlockM4Nv1_ijk_M32Parallel	20,194.872 us	99.04	258.44
BlockM4Nv1_ikj_M32	162,750.319 us	12.29	32.07
BlockM4Nv1_ikj_M32Parallel	19,956.415 us	100.22	261.53
BlockM4Nv1_ikj_M4	162,956.796 us	12.27	32.03
BlockM4Nv1_ikj_M4Parallel	22,891.121 us	87.37	228.00
BlockM4Nv3_ikj_M4	83,848.352 us	23.85	62.25
BlockM4Nv3_ikj_M4Parallel	15,405.993 us	129.82	338.78
BlockM4Nv3_ikj_M32	80,101.363 us	24.97	65.16
BlockM4Nv3_ikj_M32Parallel	13,063.492 us	153.10	399.53
UseMathNet	62,791.007 us	31.85	83.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类型的基准测试结果如下。

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
Basic	1024	3,754,237.106 us	2,074.9078 us	3,041.3732 us	3,753,734.041 us	1.000	0.00
TileRow	1024	840,254.040 us	142.7877 us	200.1684 us	840,254.500 us	0.224	0.00
TileRowSpan	1024	761,125.180 us	498.8317 us	715.4096 us	760,974.583 us	0.203	0.00
TileRowRef	1024	330,818.469 us	106.5126 us	152.7572 us	330,825.885 us	0.088	0.00
TileRowSimd	1024	117,413.868 us	90.4758 us	129.7577 us	117,410.512 us	0.031	0.00
TileRowSimdFmaBcl	1024	117,210.060 us	67.9110 us	99.5431 us	117,214.725 us	0.031	0.00
TileRowSimdFma	1024	117,189.597 us	64.3478 us	94.3201 us	117,187.067 us	0.031	0.00
TileRowSimdParallel	1024	26,510.611 us	185.3270 us	277.3887 us	26,510.653 us	0.007	0.00
TileRowSimdFmaParallel	1024	25,593.386 us	211.9040 us	317.1679 us	25,553.712 us	0.007	0.00
BlockM4Nv1_ijk_M32	1024	152,420.390 us	107.9145 us	161.5213 us	152,420.901 us	0.041	0.00
BlockM4Nv1_ijk_M32Parallel	1024	38,359.750 us	243.9897 us	342.0395 us	38,300.182 us	0.010	0.00
BlockM4Nv1_ikj_M32	1024	148,065.516 us	119.2124 us	178.4315 us	148,031.271 us	0.039	0.00
BlockM4Nv1_ikj_M32Parallel	1024	36,589.371 us	249.5473 us	349.8305 us	36,578.652 us	0.010	0.00
BlockM4Nv1_ikj_M4	1024	146,542.856 us	441.5088 us	647.1579 us	146,317.281 us	0.039	0.00
BlockM4Nv1_ikj_M4Parallel	1024	36,993.328 us	268.5206 us	393.5940 us	37,065.863 us	0.010	0.00
BlockM4Nv3_ikj_M4	1024	489,578.227 us	2,903.1298 us	4,069.7828 us	493,334.646 us	0.130	0.00
BlockM4Nv3_ikj_M4Parallel	1024	17,373.172 us	80.5072 us	115.4610 us	17,399.023 us	0.005	0.00
BlockM4Nv3_ikj_M32	1024	487,247.653 us	1,268.0786 us	1,777.6693 us	488,772.104 us	0.130	0.00
BlockM4Nv3_ikj_M32Parallel	1024	17,908.937 us	119.7722 us	179.2693 us	17,917.776 us	0.005	0.00
BlockM4Nv3On512_ikj_M4	1024	NA	NA	NA	NA	?	?
BlockM4Nv3On512_ikj_M4Parallel	1024	NA	NA	NA	NA	?	?
BlockM4Nv3On512_ikj_M32	1024	NA	NA	NA	NA	?	?
BlockM4Nv3On512_ikj_M32Parallel	1024	NA	NA	NA	NA	?	?
UseMathNet	1024	154,141.940 us	837.3170 us	1,253.2568 us	153,909.807 us	0.041	0.00

测试结果对比如下。

Method	耗时	GFOPS	性能倍数
Basic	3,754,237.106 us	0.53	1.00
TileRow	840,254.040 us	2.38	4.47
TileRowSpan	761,125.180 us	2.63	4.93
TileRowRef	330,818.469 us	6.05	11.35
TileRowSimd	117,413.868 us	17.03	31.97
TileRowSimdFmaBcl	117,210.060 us	17.06	32.03
TileRowSimdFma	117,189.597 us	17.07	32.04
TileRowSimdParallel	26,510.611 us	75.44	141.61
TileRowSimdFmaParallel	25,593.386 us	78.15	146.69
BlockM4Nv1_ijk_M32	152,420.390 us	13.12	24.63
BlockM4Nv1_ijk_M32Parallel	38,359.750 us	52.14	97.87
BlockM4Nv1_ikj_M32	148,065.516 us	13.51	25.36
BlockM4Nv1_ikj_M32Parallel	36,589.371 us	54.66	102.60
BlockM4Nv1_ikj_M4	146,542.856 us	13.65	25.62
BlockM4Nv1_ikj_M4Parallel	36,993.328 us	54.06	101.48
BlockM4Nv3_ikj_M4	489,578.227 us	4.09	7.67
BlockM4Nv3_ikj_M4Parallel	17,373.172 us	115.12	216.09
BlockM4Nv3_ikj_M32	487,247.653 us	4.10	7.70
BlockM4Nv3_ikj_M32Parallel	17,908.937 us	111.68	209.63
UseMathNet	154,141.940 us	12.98	24.36

我们表现最好的算法是 BlockM4Nv3_ikj_M4Parallel，它的性能是基础算法的 216.09 倍，此时的每秒浮点计算次数是 115.12 GFOPS。

AMD Ryzen 7 7840H 的Single运算速度是458.39 GFOPS。故对于单精度浮点数（Single）的运算，现代CPU均能达到每秒千亿次（100GFOPS）级别的运算速度。

5.3.2 Double（双精度浮点型）

Arm架构上、.NET 9.0 程序、Double类型的基准测试结果如下。

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
Basic	1024	3,979,217.591 us	3,084.0733 us	4,423.0862 us	3,979,746.459 us	1.000	0.00
TileRow	1024	880,127.965 us	571.4544 us	855.3261 us	880,017.875 us	0.221	0.00
TileRowSpan	1024	797,048.246 us	486.5361 us	713.1583 us	797,011.417 us	0.200	0.00
TileRowRef	1024	348,452.077 us	324.7616 us	465.7634 us	348,371.812 us	0.088	0.00
TileRowSimd	1024	242,866.758 us	231.8556 us	325.0292 us	242,941.903 us	0.061	0.00
TileRowSimdFmaBcl	1024	241,862.247 us	306.1996 us	429.2491 us	241,910.333 us	0.061	0.00
TileRowSimdFma	1024	241,918.592 us	195.5033 us	292.6202 us	241,910.924 us	0.061	0.00
TileRowSimdParallel	1024	57,335.473 us	813.2521 us	1,192.0544 us	57,027.893 us	0.014	0.00
TileRowSimdFmaParallel	1024	55,068.792 us	505.7084 us	741.2608 us	54,973.000 us	0.014	0.00
BlockM4Nv1_ijk_M32	1024	268,039.331 us	388.0223 us	580.7736 us	267,957.021 us	0.067	0.00
BlockM4Nv1_ijk_M32Parallel	1024	65,948.754 us	894.1736 us	1,338.3570 us	66,392.143 us	0.017	0.00
BlockM4Nv1_ikj_M32	1024	265,015.397 us	231.9062 us	332.5930 us	265,013.260 us	0.067	0.00
BlockM4Nv1_ikj_M32Parallel	1024	65,237.383 us	958.7143 us	1,434.9585 us	65,071.180 us	0.016	0.00
BlockM4Nv1_ikj_M4	1024	257,808.595 us	340.6194 us	488.5062 us	257,729.927 us	0.065	0.00
BlockM4Nv1_ikj_M4Parallel	1024	65,421.381 us	691.5592 us	1,013.6786 us	65,436.662 us	0.016	0.00
BlockM4Nv3_ikj_M4	1024	594,005.257 us	137,852.8179 us	206,331.6239 us	694,300.630 us	0.149	0.05
BlockM4Nv3_ikj_M4Parallel	1024	35,414.032 us	584.9320 us	875.4988 us	35,398.204 us	0.009	0.00
BlockM4Nv3_ikj_M32	1024	587,200.067 us	134,524.2102 us	201,349.5203 us	685,058.547 us	0.148	0.05
BlockM4Nv3_ikj_M32Parallel	1024	35,606.024 us	439.7341 us	630.6535 us	35,467.028 us	0.009	0.00
BlockM4Nv3On512_ikj_M4	1024	NA	NA	NA	NA	?	?
BlockM4Nv3On512_ikj_M4Parallel	1024	NA	NA	NA	NA	?	?
BlockM4Nv3On512_ikj_M32	1024	NA	NA	NA	NA	?	?
BlockM4Nv3On512_ikj_M32Parallel	1024	NA	NA	NA	NA	?	?
UseMathNet	1024	172,491.316 us	1,034.3648 us	1,415.8503 us	172,643.514 us	0.043	0.00

测试结果对比如下。

Method	耗时	GFOPS	性能倍数
Basic	3,979,217.591 us	0.50	1.00
TileRow	880,127.965 us	2.27	4.52
TileRowSpan	797,048.246 us	2.51	4.99
TileRowRef	348,452.077 us	5.74	11.42
TileRowSimd	242,866.758 us	8.23	16.38
TileRowSimdFmaBcl	241,862.247 us	8.27	16.45
TileRowSimdFma	241,918.592 us	8.27	16.45
TileRowSimdParallel	57,335.473 us	34.88	69.40
TileRowSimdFmaParallel	55,068.792 us	36.32	72.26
BlockM4Nv1_ijk_M32	268,039.331 us	7.46	14.85
BlockM4Nv1_ijk_M32Parallel	65,948.754 us	30.33	60.34
BlockM4Nv1_ikj_M32	265,015.397 us	7.55	15.02
BlockM4Nv1_ikj_M32Parallel	65,237.383 us	30.66	61.00
BlockM4Nv1_ikj_M4	257,808.595 us	7.76	15.43
BlockM4Nv1_ikj_M4Parallel	65,421.381 us	30.57	60.82
BlockM4Nv3_ikj_M4	594,005.257 us	3.37	6.70
BlockM4Nv3_ikj_M4Parallel	35,414.032 us	56.47	112.36
BlockM4Nv3_ikj_M32	587,200.067 us	3.41	6.78
BlockM4Nv3_ikj_M32Parallel	35,606.024 us	56.17	111.76
UseMathNet	172,491.316 us	11.59	23.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，就是因为这些原因。

希望我的文章能给读者们带来启发。

附录

posted on 2025-08-31 17:35 zyl910 阅读(41) 评论(1) 收藏举报

刷新页面返回顶部

zyl910

公告

一、基本算法与测试办法

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（单精度浮点型）

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

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

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

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