[C#] 复数乘法的跨平台SIMD硬件加速向量算法（不但支持X86的Sse、Avx、Avx5
将复数乘法改造为SIMD向量算法，是稍微有一些的难度的。首个难点是需要重新调整元素的位置，才能满足复数乘法公式。而“调整元素的位置”与内存中数据布局有关，不同办法的性能不同。还需考虑优化内存访问等细节。最近知乎有个 帖子 讨论了该话题，且 hez2010 给出了修正后的基于Avx指令集 HorizontalAdd（程度加法）的向量算法。
于是我来说说基于 Shuffle（换位）的向量算法吧。且这些算法是跨平台的，同一份源代码，能在 X86（Sse、Avx等指令集）及Arm（AdvSimd指令集）等架构上运行，且均享有SIMD硬件加速。本文还介绍了512位向量算法，用于对比测试Avx512指令集的硬件加速。
一、简朴算法
先回顾一下 .NET 里怎么做复数乘法。
从 .NET 4.0 开始，提供了 Complex范例。于是不必手工地根据数学知识来编写复数乘法函数，而是可以使用Complex范例来做复数运算，简化了不少。
为了便于进行基准测试，可以将一个数组作为输入参数，随后对各个元素自乘并进行累加。最后返回累加的结果。
public static Complex mul(Complex[] a) {
Complex c = 0;
for (int i = 0; i < a.Length; i++) {
c += a * a;
}
return c;
}二、向量算法
2.1 算法思路
2.1.1 复数乘法的数学定义
向量范例并未提供复数乘法的方法，于是需要手工地根据数学知识来编写复数乘法函数了。
先来回顾一下复数乘法的数学定义。
(a + bi)*(c + di) = (ac – bd) + (ad + bc)ihttps://img2024.cnblogs.com/blog/169212/202412/169212-20241228005637218-1865452110.png
a + bi 是乘法的左侧参数， c + di 是乘法的右侧参数。
数学表达式一样平常喜欢省略“乘法运算符”。例如上式里，“ac”其实表示“a*c”。为了能使乘法运算更清晰，上式可以写成下面的形式。
(a + bi)*(c + di) = (a*c – b*d) + (a*d + b*c)i可以看出，复数乘法是由4次实数乘法，以及若干个加减运算所构成。
2.1.2 复数的数据布局
Complex 是一个布局体，此中次序存放了2个Double范例的成员，分别为“实部”与“虚部”。Double范例是64位的，2个Double就是128位，于是 Complex范例是128位的，即16字节。
由于大多数架构都支持128位向量的SIMD硬件加速。所以 C# 步伐在使用Complex范例时，其实已经享受到SIMD硬件加速了。这就是为什么手写的向量算法，偶然还比不过 Complex的原因。于是需要更细致的进行优化。
对于 Complex数组，数据是连续存放的。于是使用向量范例从Complex数组加载数据时，能加载到整数个Complex。
[*]对于128位向量，能存放1个Complex。以Double元素的视角来看，向量中的元素依次是 。注：“a”代表复数的“实部”，“b”代表复数的“虚部”。
[*]对于256位向量，能存放2个Complex。以Double元素的视角来看，向量中的元素依次是 。
[*]对于512位向量，能存放4个Complex。以Double元素的视角来看，向量中的元素依次是 。
由于数据都是这样连续存放的，对于上面这种数据布局，本文将它简称为 a + bi 形式。后面将会经常使用这种简称。
例如将实部与虚部交换，变为 b + ai 形式。那么代表的是下面这样的数据布局。
[*]128位向量：。
[*]256位向量：。
[*]512位向量：。
2.1.3 第1步：计算 (a*c) + (-b*d)i
先来观察一下 向量乘法（Vector.Multiply）对复数数据的运算结果。
假设已经将复数数据分别加载到向量a（a + bi）、向量c（c + di）之中，随后对这2个向量进行向量乘法运算。由于向量乘法依然是对逐个地对相同位置的元素做乘法处理，所以计算结果为 (a*c) + (b*d)i。可以观察到，它正好包含了“复数乘法内部的4个实数乘法”中的2项——a*c、b*d。
表示这个思路是精确的。但是可以注意到，复数乘法需要使用 -b*d，而上面的 b*d是不带“负号”的。于是我们需要对数据进行进一步的处理，将奇数位置（虚部）的元素做一次“求负”处理。Vector范例里没有单独提供这种处理的方法，于是需要自行编写算法。
对于“求负”处理，性能最高的办法是直接翻转浮点范例的符号位。在IEEE754浮点数标准里，规定了最高位是符号位，该位为0时表示正数，为1时表示负数。于是IEEE754浮点数标准里还能表达 -0.0（负零），它的最高位符号位为1，剩余数据位（阶码s、尾数m）均为0。
二进制的Xor（异或运算）可以使相干二进制进行翻转。于是，将某个浮点数，与 -0.0（负零）执行Xor运算，就能将符号位翻转，这正好是“求负”运算。
由于仅需对奇数位置求负，于是我们可以预先准备一个向量mask，它的偶数位置为 0（正零），奇数位置为 -0.0（负零）。向量c（c + di）与mask进行Xor运算后，结果是 c + (-d)i。随后再与向量a（a + bi）执行乘法，结果是 (a*c) + (-b*d)i。满足所需。
Vector<double> mask = Vectors.CreateRotate(0.0, -0.0);
Vector<double> a = *p; // a + bi
var c = a; // c + di
var e = Vector.Multiply(a, Vector.Xor(c, mask)); // (a*c) + (-b*d)i由于 Vector 的长度不是固定的，手工给mask赋值不太方便。于是这里使用 VectorTraits 库的 Vectors.CreateRotate 方法来简化赋值，它会循环将各个参数依次放置在向量的各个元素里。从而满足了“偶数为正零，奇数为负零”的要求。
Vector a = *p 表示根据指针p，将数据加载到向量a。由于是复数数据，于是a的数学意义为 a + bi，即复数乘法的左侧参数。
为了与前面的“简朴算法”保持同等，于是将a复制给了c。此时c的数学意义为 c + di，即复数乘法的右侧参数。
随后根据上面的履历，算出 (a*c) + (-b*d)i，赋值给向量e。
2.1.4 第2步：计算 (a*d) + (b*c)i
根据上一节的履历，使用2个步骤就能计算出 (a*d) + (b*c)i——
[*]将c的实部与虚部交换，得到 (d) + (c)i，赋值给向量f。
[*]对 a、f 执行向量乘法（Vector.Multiply），便能得到(a*d) + (b*c)i。
第2步容易实现，但第1步遇到了困难——.NET 的向量方法里，没有合适的方法来做这一工作。
从.NET 7.0开始，Vector128等向量范例增加了 Shuffle 方法。用该方法，可以给向量内的元素进行换位。但是直至 .NET 8.0，Shuffle都没有硬件加速，而是使用了标量回退代码。
为相识决 Shuffle 方法没有硬件加速的问题，我开发了VectorTraits 库。它使用了各个架构的shuffle类别的指令，从而使 Shuffle 方法具有硬件加速。具体来说，它分别使用了以下指令。
[*]X86: 使用 _mm_shuffle_epi8等指令.
[*]Arm: 使用 vqvtbl1q_u8 指令.
[*]Wasm: 使用 i8x16.swizzle 指令.
VectorTraits 不但为固定巨细的向量范例（如 Vector128）提供了Shuffle方法，它还为主动巨细的向量范例（Vector）也提供了Vector方法。
虽然VectorTraits的Shuffle方法是有硬件加速的，但它不是最佳办法。VectorTraits还提供了YShuffleG2方法，专门用来处理2元素组的换位。它用起来更简朴，且大多数时间的性能更高。
YShuffleG2方法通过ShuffleControlG2参数来控制换位模式，例如 ShuffleControlG2.YX 表示将“XY”次序给换成“YX”次序，即我们所需的“实部与虚部交换”。
根据这些信息，便能完成第2步的代码了。
var f = Vectors.YShuffleG2(c, ShuffleControlG2.YX); // (d) + (c)i
f = Vector.Multiply(a, f); // (a*d) + (b*c)i2.1.5 第3步：计算结果归并
颠末上面的步骤，已经算出了 (a*c) + (-b*d)i、(a*d) + (b*c)i。复数乘法规则是 (a*c – b*d) + (a*d + b*c)i，实数乘法的步骤均处理完了，还剩下加法与数据变换的处理。
向量加法与向量乘法一样，是对逐个地对相同位置的元素做相加处理。于是我们需要将上面的那2组数据，变换为 (a*c) + (a*d)i、(-b*d) + (b*c)i，这样才好交给向量加法去处理。
可以观察到，这种变换，就是 2x2矩阵的转置操作。将数据写成矩阵形式，便能清晰看出它是转置操作。
[(a*c) (-b*d)] --> [(a*c) (a*d)]
[(a*d) (b*c)] --> [(-b*d) (b*c)]VectorTraits库提供了YGroup2Transpose方法，用它便能实现2x2矩阵的转置操作。
var g = Vectors.YGroup2Transpose(e, f, out var h); // g is {(a*c) + (a*d)i}; h is {(-b*d) + (b*c)i}
g += h; // (a*c - b*d) + (a*d + b*c)i自此，完成了复数乘法的计算。
2.2 算法实现（UseVectors）
有了上面的思路，便能编写出对数组计算向量乘法累加的算法了。源代码如下。
public static unsafe Complex UseVectorsDo(Complex[] numbers) {
int blockWidth = Vector<double>.Count / 2; // Complex is double*2
int cntBlock = numbers.Length / blockWidth;
int cntRem = numbers.Length - (cntBlock * blockWidth);
// -- Processs body.
Complex result;
Vector<double> acc = Vector<double>.Zero;
Vector<double> mask = Vectors.CreateRotate(0.0, -0.0);
fixed (Complex* pnumbers = numbers) {
Vector<double>* p = (Vector<double>*)pnumbers; // Set pointer to numbers.
Vector<double>* pEnd = p + cntBlock;
while (p < pEnd) {
// -- Complex multiply: (a + bi)*(c + di) = (ac – bd) + (ad + bc)i
Vector<double> a = *p; // a + bi
var c = a; // c + di
var e = Vector.Multiply(a, Vector.Xor(c, mask)); // (a*c) + (-b*d)i
var f = Vectors.YShuffleG2(c, ShuffleControlG2.YX); // (d) + (c)i
f = Vector.Multiply(a, f); // (a*d) + (b*c)i
var g = Vectors.YGroup2Transpose(e, f, out var h); // g is {(a*c) + (a*d)i}; h is {(-b*d) + (b*c)i}
g += h; // (a*c - b*d) + (a*d + b*c)i
// Sum
acc += g;
// Next
++p;
}
// Vector to a Complex.
double re = 0.0, im = 0.0;
for (int i = 0; i < Vector<double>.Count; i += 2) {
re += acc;
im += acc;
}
result = new Complex(re, im);
// -- Processs remainder.
Complex* q = (Complex*)pEnd;
for (int i = 0; i < cntRem; i++) {
result += (*q) * (*q);
// Next
++q;
}
}
return result;
}该方法是用指针来处理数据的。代码分为4个部门。
[*]最前面的部门，是变量初始化。例如计算 cntBlock（有多少个块）、cntRem（剩余部门有多少个复数）等。
[*]“Processs body”是向量处理的主体代码，它构造好mask，并计算好指针地址，随后循环处理主要数据（向量范例的整数倍的数据）。用的就是上面提到的算法。
[*]“Vector to a Complex”是将向量里存放的复数数据，转换为单个复数。
[*]“Processs remainder”是剩余部门的处理。
2.3 深入探讨
2.3.1 YGroup2Transpose的深入探讨
此时，有些读者可能会产生疑问——Vector范例有可能是256位（CPU支持Avx2指令集时），这时向量里不只是2个Double，而是4个Double，YGroup2Transpose还能正常工作吗？
答案是——YGroup2Transpose当然能正常工作。
当 Vector 为256位时，4个Double被分为了2组。既可以看作有2组2x2矩阵，随后分别进行了矩阵转置处理。正好Complex范例是128位的，由2个Double所构成，与2x2矩阵相匹配。
从中可以看出规律——假设向量里可以存放 N*2 个元素（如Double元素），那么 YGroup2Transpose可以对 N组2x2矩阵进行转置。且“N个Complex”做复数乘法时，能使用YGroup2Transpose方法。
对于X86架构，YGroup2Transpose是用Shuffle类别的指令来实现的。对于Arm架构，它是使用 AdvSimd中转置类指令 TRN1、TRN2 来实现的，服从很高。
另注：为了使方法名的可读性更高，未来筹划将 YGroup2Transpose 改名为 YMatrix2Transpose。初步筹划将在VectorTraits库下一个大版本，做这个改名。
2.3.2 HorizontalAdd（程度加法）算法与本算法的区别
HorizontalAdd（程度加法）算法与本算法是非常相似的，仅“第3步：计算结果归并”不同。
先来看看128位向量时的运算环境。HorizontalAdd 会将 相邻2个元素相加，并把结果放在1个元素里，于是2个向量在处理后被归并成了1个向量。
颠末前面2步，已经算出了 (a*c) + (-b*d)i、(a*d) + (b*c)i。此时进行 HorizontalAdd，便会算出 (a*c – b*d) + (a*d + b*c)i，正好是复数乘法的运算结果！
随后改为用Avx指令集的256位向量来实现，虽然也能精确算出结果，但其实此时 HorizontalAdd 是Avx指令集的特别版本——它不是对整个向量进行程度加法的，而是按每128位小道（lane）来做程度加法的。
例如Double范例的源数据是 、，那么这2种程度加法的结果是不同的——
[*]整256位向量：
[*]每128位小道：
第1种方式（整256位向量）比力符合常规思路，但第2种方式（ 每128位小道）在很多时间更有用。例如复数范例Complex是128位，要使Complex的运算结果精确，于是需要第2种方式（ 每128位小道）的程度加法。幸好Avx指令集，提供的就是第2种方式的程度加法，从而方便了计算。
虽然仅需要1条程度加法指令就能实现“第3步：计算结果归并”，但由于该指令牵涉2个向量的程度操作，所以处理器会稍微多花一些时间。下面摘录了Intel手册对程度加法指令（vhaddpd）的说明，“Latency and Throughput”就是介绍延迟与吞吐率的，可以发现它的值比力高。
__m256d _mm256_hadd_pd (__m256d a, __m256d b)
#include <immintrin.h>
Instruction: vhaddpd ymm, ymm, ymm
CPUID Flags: AVX
Description
Horizontally add adjacent pairs of double-precision (64-bit) floating-point elements in a and b, and pack the results in dst.
Operation
dst := a + a
dst := b + b
dst := a + a
dst := b + b
dst := 0
Latency and Throughput
Architecture Latency Throughput (CPI)
Alderlake 5 2
Icelake Intel Core 6 2
Icelake Xeon 6 2
Sapphire Rapids 5 2
Skylake 7 2对于今世处理器，程度加法算法与本算法的性能，一样平常是是差不多的。详见“三、基准测试结果”。
而且Avx512系列指令集里，尚未提供512位的程度加法指令，故更保举使用本算法来处理复数加法。
2.4 用引用代替指针, 摆脱unsafe关键字（UseVectorsSafeDo）
从 C# 7.3开始，可以用引用代替指针, 摆脱unsafe关键字与fixed语句。其实编程思路与指针差不多的，只是一些地方需要换一种写法。具体办法可参考《C# 使用SIMD向量范例加速浮点数组求和运算（4）：用引用代替指针, 摆脱unsafe关键字，兼谈Unsafe类的使用》。
改造后的代码如下。
public static Complex UseVectorsSafeDo(Complex[] numbers) {
int blockWidth = Vector<double>.Count / 2; // Complex is double*2
int cntBlock = numbers.Length / blockWidth;
int cntRem = numbers.Length - (cntBlock * blockWidth);
// -- Processs body.
Vector<double> acc = Vector<double>.Zero;
Vector<double> mask = Vectors.CreateRotate(0.0, -0.0);
ref Vector<double> p = ref Unsafe.As<Complex, Vector<double>>(ref numbers); // Set pointer to numbers.
ref Vector<double> pEnd = ref Unsafe.Add(ref p, cntBlock);
while (Unsafe.IsAddressLessThan(ref p, ref pEnd)) {
// -- Complex multiply: (a + bi)*(c + di) = (ac – bd) + (ad + bc)i
Vector<double> a = p; // a + bi
var c = a; // c + di
var e = Vector.Multiply(a, Vector.Xor(c, mask)); // (a*c) + (-b*d)i
var f = Vectors.YShuffleG2(c, ShuffleControlG2.YX); // (d) + (c)i
f = Vector.Multiply(a, f); // (a*d) + (b*c)i
var g = Vectors.YGroup2Transpose(e, f, out var h); // g is {(a*c) + (a*d)i}; h is {(-b*d) + (b*c)i}
g += h; // (a*c - b*d) + (a*d + b*c)i
// Sum
acc += g;
// Next
p = ref Unsafe.Add(ref p, 1);
}
// Vector to a Complex.
double re = 0.0, im = 0.0;
for (int i = 0; i < Vector<double>.Count; i += 2) {
re += acc;
im += acc;
}
Complex result = new Complex(re, im);
// -- Processs remainder.
ref Complex q = ref Unsafe.As<Vector<double>, Complex>(ref pEnd);
for (int i = 0; i < cntRem; i++) {
result += q * q;
// Next
q = ref Unsafe.Add(ref q, 1);
}
return result;
}上面的代码，与指针版代码（UseVectorsDo）是等价的。步伐的性能，也是差不多的。
2.5 循环展开（UseVectorsX2Do）
对于小循环，循环时的跳转处理也是一笔不小的开销。此时可以使用循环展开的办法，在循环内处理多条数据，从而使循环开销的占比减低。优化了性能。
这里选用了2倍循环展开。它还能带来一个利益，就是能将2元素组换位（YShuffleG2），改为4元素组换位（YShuffleG4X2）。因为一些架构上有“4元素组换位”的专业指令（例如 Avx2.Permute4x64），性能很高。
当初因为主动巨细向量偶然只有128位，只能存放2个Double，无法满足“4元素组换位”的要求，故审慎的使用了YShuffleG2方法。而如今有了2倍数据，即使主动巨细向量只有128位，也能保证至少共有4个元素，故可以使用 YShuffleG4X2方法。
YShuffleG4X2方法通过ShuffleControlG4参数来控制换位模式，例如 ShuffleControlG4.YXWZ 表示将“XYZW”次序给换成“YXWZ”次序，即我们所需的“实部与虚部交换”。（其实， ShuffleControlG4.YXWZ 就是Avx2.Permute4x64 指令的参数 0b10110001，如今用枚举来描述，代码的可读性更好）
由于如今所用的ShuffleControlG4参数是固定的常数，于是还可以使用 YShuffleG4X2_Const，它的性能一样平常更好。
根据上面的履历，便能编写出2倍循环展开时的算法了。源代码如下。
public static Complex UseVectorsX2Do(Complex[] numbers) {
const int batchWidth = 2; // X2
int blockWidth = Vector<double>.Count * batchWidth / 2; // Complex is double*2
int cntBlock = numbers.Length / blockWidth;
int cntRem = numbers.Length - (cntBlock * blockWidth);
// -- Processs body.
Vector<double> acc = Vector<double>.Zero;
Vector<double> acc1 = Vector<double>.Zero;
Vector<double> mask = Vectors.CreateRotate(0.0, -0.0);
ref Vector<double> p = ref Unsafe.As<Complex, Vector<double>>(ref numbers); // Set pointer to numbers.
ref Vector<double> pEnd = ref Unsafe.Add(ref p, cntBlock * batchWidth);
while (Unsafe.IsAddressLessThan(ref p, ref pEnd)) {
// -- Complex multiply: (a + bi)*(c + di) = (ac – bd) + (ad + bc)i
Vector<double> a0 = p; // a + bi
var a1 = Unsafe.Add(ref p, 1);
var c0 = a0; // c + di
var c1 = a1;
var e0 = Vector.Multiply(a0, Vector.Xor(c0, mask)); // (a*c) + (-b*d)i
var e1 = Vector.Multiply(a1, Vector.Xor(c1, mask));
var f0 = Vectors.YShuffleG4X2_Const(c0, c1, ShuffleControlG4.YXWZ, out var f1); // (d) + (c)i
f0 = Vector.Multiply(a0, f0); // (a*d) + (b*c)i
f1 = Vector.Multiply(a1, f1);
var g0 = Vectors.YGroup2Transpose(e0, f0, out var h0); // g is {(a*c) + (a*d)i}; h is {(-b*d) + (b*c)i}
var g1 = Vectors.YGroup2Transpose(e1, f1, out var h1);
g0 += h0; // (a*c - b*d) + (a*d + b*c)i
g1 += h1;
// Sum
acc += g0;
acc1 += g1;
// Next
p = ref Unsafe.Add(ref p, batchWidth);
}
acc += acc1;
// Vector to a Complex.
double re = 0.0, im = 0.0;
for (int i = 0; i < Vector<double>.Count; i += 2) {
re += acc;
im += acc;
}
Complex result = new Complex(re, im);
// -- Processs remainder.
ref Complex q = ref Unsafe.As<Vector<double>, Complex>(ref pEnd);
for (int i = 0; i < cntRem; i++) {
result += q * q;
// Next
q = ref Unsafe.Add(ref q, 1);
}
return result;
}2.6 512位算法（UseVector512sX2Do）
.NET 8.0 新增了对 X86架构的Avx512系列指令集的支持，且新增了 Vector512范例。
VectorTraits 3.0版已经支持了Avx512系列指令集，于是能够很容易将主动巨细向量的算法，改造为512位向量的算法。只需要做文本替换，将“Vector”替换为“Vector512”，便基本完成了改造，顶多有少量报错需修正。而不用学习复杂的Avx512指令集，大大降低了门槛。源代码如下。
#if NET8_0_OR_GREATER
public void UseVector512sX2() {
if (!Vector512s.IsHardwareAccelerated) throw new NotSupportedException("Vector512 does not have hardware acceleration!");
_destination = UseVector512sX2Do(_array);
}
public static Complex UseVector512sX2Do(Complex[] numbers) {
const int batchWidth = 2; // X2
int blockWidth = Vector512<double>.Count * batchWidth / 2; // Complex is double*2
int cntBlock = numbers.Length / blockWidth;
int cntRem = numbers.Length - (cntBlock * blockWidth);
// -- Processs body.
Vector512<double> acc = Vector512<double>.Zero;
Vector512<double> acc1 = Vector512<double>.Zero;
Vector512<double> mask = Vector512s.CreateRotate(0.0, -0.0);
ref Vector512<double> p = ref Unsafe.As<Complex, Vector512<double>>(ref numbers); // Set pointer to numbers.
ref Vector512<double> pEnd = ref Unsafe.Add(ref p, cntBlock * batchWidth);
while (Unsafe.IsAddressLessThan(ref p, ref pEnd)) {
// -- Complex multiply: (a + bi)*(c + di) = (ac – bd) + (ad + bc)i
Vector512<double> a0 = p; // a + bi
var a1 = Unsafe.Add(ref p, 1);
var c0 = a0; // c + di
var c1 = a1;
var e0 = Vector512s.Multiply(a0, Vector512s.Xor(c0, mask)); // (a*c) + (-b*d)i
var e1 = Vector512s.Multiply(a1, Vector512s.Xor(c1, mask));
var f0 = Vector512s.YShuffleG4X2_Const(c0, c1, ShuffleControlG4.YXWZ, out var f1); // (d) + (c)i
f0 = Vector512s.Multiply(a0, f0); // (a*d) + (b*c)i
f1 = Vector512s.Multiply(a1, f1);
var g0 = Vector512s.YGroup2Transpose(e0, f0, out var h0); // g is {(a*c) + (a*d)i}; h is {(-b*d) + (b*c)i}
var g1 = Vector512s.YGroup2Transpose(e1, f1, out var h1);
g0 += h0; // (a*c - b*d) + (a*d + b*c)i
g1 += h1;
// Sum
acc += g0;
acc1 += g1;
// Next
p = ref Unsafe.Add(ref p, batchWidth);
}
acc += acc1;
// Vector to a Complex.
double re = 0.0, im = 0.0;
for (int i = 0; i < Vector512<double>.Count; i += 2) {
re += acc;
im += acc;
}
Complex result = new Complex(re, im);
// -- Processs remainder.
ref Complex q = ref Unsafe.As<Vector512<double>, Complex>(ref pEnd);
for (int i = 0; i < cntRem; i++) {
result += q * q;
// Next
q = ref Unsafe.Add(ref q, 1);
}
return result;
}
#endif // NET8_0_OR_GREATER注意，从.NET 8.0才开始支持 Vector512，故需要使用条件编译符号NET8_0_OR_GREATER进行判断。
由于如今有不少处理器尚未支持 Avx512系列指令集，于是需要用if语句判断一下“Vector512s.IsHardwareAccelerated”是否为真。否则就不要执行了。
三、基准测试结果
3.1 X86 架构
X86架构下的基准测试结果如下。
BenchmarkDotNet v0.14.0, Windows 11 (10.0.26100.2605)
AMD Ryzen 7 7840H w/ Radeon 780M Graphics, 1 CPU, 16 logical and 8 physical cores
.NET SDK 9.0.101
: .NET 8.0.11 (8.0.1124.51707), X64 RyuJIT AVX-512F+CD+BW+DQ+VL+VBMI
DefaultJob : .NET 8.0.11 (8.0.1124.51707), X64 RyuJIT AVX-512F+CD+BW+DQ+VL+VBMI
| Method | Count | Mean | Error | StdDev | Ratio | Code Size |
|----------------- |------ |---------:|---------:|---------:|------:|----------:|
| CallMul | 65536 | 44.45 us | 0.329 us | 0.308 us |1.00 | 128 B |
| CallMul2 | 65536 | 92.50 us | 0.104 us | 0.087 us |2.08 | NA |
| UseVectors | 65536 | 22.48 us | 0.068 us | 0.061 us |0.51 | NA |
| UseVectorsSafe | 65536 | 22.47 us | 0.084 us | 0.070 us |0.51 | NA |
| UseVectorsX2 | 65536 | 17.95 us | 0.080 us | 0.075 us |0.40 | NA |
| UseVector512sX2| 65536 | 17.26 us | 0.179 us | 0.167 us |0.39 | NA |
| Hez2010Simd_Mul2 | 65536 | 23.69 us | 0.206 us | 0.193 us |0.53 | NA |
| Hez2010Simd | 65536 | 23.01 us | 0.151 us | 0.134 us |0.52 | 298 B |方法说明——
[*]CallMul: 简朴算法。
[*]CallMul2: 知乎帖子提出给出的有问题的Avx向量算法。
[*]UseVectors: 指针写法的向量算法（2.2 算法实现）。
[*]UseVectorsSafe: 引用写法的安全向量算法（2.4 用引用代替指针, 摆脱unsafe关键字）。
[*]UseVectorsX2: 2倍循环展开的向量算法（2.5 循环展开）。
[*]UseVector512sX2: 2倍循环展开的512位向量算法（2.6 512位算法）。
[*]Hez2010Simd_Mul2: hez2010将CallMul2修正后的向量算法。
[*]Hez2010Simd: hez2010的安全向量算法。
如今来对比一下各个方法的性能。
[*]CallMul2: 44.45/92.50 ≈ 0.4805（倍）。
[*]UseVectors: 44.45/22.48 ≈ 1.9773（倍）。
[*]UseVectorsSafe: 44.45/22.47 ≈ 1.9782（倍）。
[*]UseVectorsX2: 44.45/17.95 ≈ 2.4763（倍）。性能是UseVectorsSafe的 22.47/17.95 ≈ 1.2518（倍）
[*]UseVector512sX2: 44.45/17.26 ≈ 2.5753（倍）。性能是UseVectorsSafe的 22.47/17.26 ≈ 1.3019（倍）
[*]Hez2010Simd_Mul2: 44.45/23.69 ≈ 1.8763（倍）。
[*]Hez2010Simd: 44.45/23.01 ≈ 1.9318（倍）。
首先可以注意到，UseVectors、UseVectorsSafe的性能几乎是一样的。这是因为不论是指针写法，还是引用写法，它们的算法是相同的，所以性能是一样的。而且若观察JIT生成汇编代码，你会发现它们基本是一样的。
其次，可以发现 Hez2010Simd_Mul2、Hez2010Simd、UseVectors、UseVectorsSafe 这4种方法的耗时很接近，都是23us左右。差距很小，可看作测试误差范围内。故可以得出结论，程度加法算法与本算法的性能是几乎是一样的。Avx是256位向量宽度，比Sse的128位向量宽度翻了一倍，理论性能是2倍。如今的测试结果，很接近这个理论值。
再来看 UseVectorsX2，会发现它的性能又有提升，比起普通向量算法（UseVectorsSafe）快了20%左右。看来此时做循环展开，确实有效。
最厥后看 UseVector512sX2。Avx512是512位向量宽度，是Sse的128位向量宽度的4倍，理论性能应该是4倍。但是实际测试时，它仅比 UseVectorsX2 稍微快了一点点。
这是因为当前处理器是 Zen4微架构的，它并没有专门的512位运算单元，而是通过2个256位运算单元组合而成的，还不能完全发挥Avx512指令集的潜力。若换成 Zen5、Sapphire Rapids等具有专门512位运算单元的微架构的处理器，能获得进一步性能提升。
3.1.1 更深入的说明
细致观察一下上面的测试结果，会发现本算法（UseVectorsSafe）比起程度加法算法（Hez2010Simd），稍微快一点点。
其实这跟CPU微架构有关。在AMD的处理器上，很多时间是本算法稍微快一点；而在Intel的处理器上，很多时间是程度加法算法稍微快一点。
而且在 Intel 处理器上测试 UseVectorsX2 算法时，偶然它的性能与 UseVectorsSafe 相差不大。这也是CPU微架构的差别。
3.2 Arm 架构
同样的源代码可以在 Arm 架构上运行。基准测试结果如下。
BenchmarkDotNet v0.14.0, macOS Sequoia 15.1.1 (24B91)
Apple M2, 1 CPU, 8 logical and 8 physical cores
.NET SDK 8.0.204
: .NET 8.0.4 (8.0.424.16909), Arm64 RyuJIT AdvSIMD
DefaultJob : .NET 8.0.4 (8.0.424.16909), Arm64 RyuJIT AdvSIMD
| Method | Count | Mean | Error | StdDev | Ratio | RatioSD |
|----------------- |------ |---------:|---------:|---------:|------:|--------:|
| CallMul | 65536 | 56.30 us | 0.051 us | 0.045 us |1.00 | 0.00 |
| CallMul2 | 65536 | NA | NA | NA | ? | ? |
| UseVectors | 65536 | 56.90 us | 0.468 us | 0.415 us |1.01 | 0.01 |
| UseVectorsSafe | 65536 | 56.32 us | 0.019 us | 0.017 us |1.00 | 0.00 |
| UseVectorsX2 | 65536 | 47.18 us | 0.025 us | 0.024 us |0.84 | 0.00 |
| UseVector512sX2| 65536 | NA | NA | NA | ? | ? |
| Hez2010Simd_Mul2 | 65536 | NA | NA | NA | ? | ? |
| Hez2010Simd | 65536 | NA | NA | NA | ? | ? |首先可以注意到，CallMul2、Hez2010Simd_Mul2、Hez2010Simd 都无法执行基准测试。这是因为它们都依靠Avx指令集，但如今是Arm架构的处理器，没有Avx指令集，于是报错了。
此时便能体现出本文介绍的算法的优势了——支持跨平台。同一份源代码，能在 X86（Sse、Avx等指令集）及Arm（AdvSimd指令集）等架构上运行，且均享有SIMD硬件加速。
UseVector512sX2是我们主动抛出了非常。虽然同一份源代码可以运行，但由于此时没有硬件加速，测试起来没有意义。故不必测试。
随厥后对比一下各个方法的性能。
[*]UseVectors: 56.30/56.90 ≈ 0.9895（倍）。
[*]UseVectorsSafe: 56.30/56.32 ≈ 0.9996（倍）。
[*]UseVectorsX2: 56.30/47.18 ≈ 1.1933（倍）。
UseVectors、UseVectorsSafe的耗时，与CallMul的结果几乎相同。前面提到过，Complex范例内部是已经使用了128位向量加速的。由于Arm架构的AdvSimd指令集是 128位的，于是此时 Vector范例的宽度也是128位的。所以此时用Vector范例实现的算法，理论上跟Complex范例的性能是一样的。
而UseVectorsX2比CallMul快20%左右。这表示此时做循环展开，确实有效。
附录
[*]完备源代码: https://github.com/zyl910/VectorTraits.Sample.Benchmarks/blob/main/VectorTraits.Sample.Benchmarks.Inc/Complexes/ComplexMultiplySumBenchmark.cs
[*]YGroup2Transpose 的文档: https://zyl910.github.io/VectorTraits_doc/api/Zyl.VectorTraits.Vectors.YGroup2Transpose.html
[*]YShuffleG2 的文档: https://zyl910.github.io/VectorTraits_doc/api/Zyl.VectorTraits.Vectors.YShuffleG2.html
[*]YShuffleG4X2_Const 的文档: https://zyl910.github.io/VectorTraits_doc/api/Zyl.VectorTraits.Vectors.YShuffleG4X2_Const.html
[*]VectorTraits 的NuGet包: https://www.nuget.org/packages/VectorTraits
[*]VectorTraits 的在线文档: https://zyl910.github.io/VectorTraits_doc/
[*]VectorTraits 源代码: https://github.com/zyl910/VectorTraits
[*]C#simd使用Avx类的代码比普通的for循环代码慢，什么原因呢？ - hez2010的回复 - 知乎
[*]C# 使用SIMD向量范例加速浮点数组求和运算（4）：用引用代替指针, 摆脱unsafe关键字，兼谈Unsafe类的使用
[*]《 24位图像程度翻转的跨平台SIMD硬件加速向量算法的关键——YShuffleX3Kernel源码解读（如Avx2解决shuffle的跨lane问题、Avx512优化等）》
出处：http://www.cnblogs.com/zyl910/ 版权声明：自由转载-非商用-非衍生-保持署名 | Creative Commons BY-NC-ND 3.0.
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。
页:
[1]