c# 使用 system.numerics 将数组元素相乘

Posted 2023-02-16

【中文标题】c# 使用 system.numerics 将数组元素相乘

【英文标题】：c# multiplying array elements using system.numerics

【发布时间】：2020-02-09 11:30:33

【问题描述】：

我正在尝试使用 System.Numerics 来处理多个数组元素。 有没有更快的方法将结果向量 (accVector) 的元素相乘？ 目前 accVector 需要转换为使用 LINQ 将元素相乘的数组。

        private double VectorMultiplication(double[] array)
        
            int vectorSize = Vector<double>.Count;
            var accVector = Vector<double>.One;
            int i;

            for (i = 0; i <= array.Length - vectorSize; i += vectorSize)
            
                var v = new Vector<double>(array, i);
                accVector = Vector.Multiply(accVector, v);
            

            var tempArray = new double[Vector<double>.Count];
            accVector.CopyTo(tempArray);
            var result = tempArray.Aggregate(1d, (p, d) => p * d);

            for (; i < array.Length; i++)
            
                result *= array[i];
            
            return result;
        

【问题讨论】：

在最后一个 for 循环中 i 从 array.Length 开始，因为您没有将 i 重置为零。你能解释一下这些乘法背后的逻辑吗？

您通常希望 CPU 执行的操作是随机播放和垂直 SIMD 相乘，直到您减少到 1 个标量结果。与水平总和相同，但使用 mul 而不是 add：Fastest way to do horizontal float vector sum on x86。 IDK 如果 C# 甚至可以作为第一步从 256 缩小到 128，但 log2(N) shuffle / mul 步骤通常是您理想的结果，而不是 store -> N-1 标量乘法，尤其是没有一个 long依赖链。

【参考方案1】：

有没有更快的方法将结果向量 (accVector) 的元素相乘？

在 Sytem.Numerics 内，没有。正如 Peter 在 cmets 中提到的，通常您首先将 256 位向量分成两个 128 位的一半并将它们相乘，然后使用 shuffle 来处理 128 位的部分。但是 System.Numerics 不提供随机播放，并且它不允许您选择正在使用的向量的大小。

通常的方法可以与 System.Runtime.Intrinsics.X86 API 一起使用，这需要 .NET Core 3.0 或更高版本。

例如：

static double product(Vector256<double> vec)

    var t = Sse2.Multiply(vec.GetLower(), vec.GetUpper());
    return t.GetElement(0) * t.GetElement(1);

这看起来可能不好，留下一个神秘的GetElement 让JIT引擎去弄清楚，但实际上codegen真的很合理：

 vmovupd     ymm0,ymmword ptr [rcx] 
 vextractf128 xmm0,ymm0,1  
 vmovupd     ymm1,ymmword ptr [rcx]  
 vmulpd      xmm0,xmm1,xmm0  
 vmovaps     xmm1,xmm0  
 vpshufd     xmm0,xmm0,0EEh  
 vmulsd      xmm0,xmm0,xmm1

所以看起来GetElement(0) 是隐式的，而GetElement(1) 的结果是vpshufd，这很好。将xmm0 复制到xmm1 而不是使用非破坏性vpshufd 有点神秘，但还不错，总体上比我通常对.NET 的预期要好。我测试了这个函数非内联，通常应该是内联，负载应该消失。

---

主循环可以改进，因为乘法的吞吐量比它的延迟要好得多。现在，乘法一次完成一次（即，一次 vector 乘法）之间有延迟（Haswell 上 5 个周期，Broadwell 上 4 个周期和更新版本）以等待前一个乘法完成，但例如英特尔 Haswell 可能会在每个周期开始两次乘法，这是 10 倍。实际上，改进不会那么大，但为指令级并行性创造一些机会会有所帮助。

例如（未测试）：

var acc0 = Vector<double>.One;
var acc1 = Vector<double>.One;
var acc2 = Vector<double>.One;
var acc3 = Vector<double>.One;
var acc4 = Vector<double>.One;
var acc5 = Vector<double>.One;
var acc6 = Vector<double>.One;
var acc7 = Vector<double>.One;
int i;

for (i = 0; i <= array.Length - vectorSize * 8; i += vectorSize * 8)

    acc0 = Vector.Multiply(acc0, new Vector<double>(array, i));
    acc1 = Vector.Multiply(acc1, new Vector<double>(array, i + vectorSize));
    acc2 = Vector.Multiply(acc2, new Vector<double>(array, i + vectorSize * 2));
    acc3 = Vector.Multiply(acc3, new Vector<double>(array, i + vectorSize * 3));
    acc4 = Vector.Multiply(acc4, new Vector<double>(array, i + vectorSize * 4));
    acc5 = Vector.Multiply(acc5, new Vector<double>(array, i + vectorSize * 5));
    acc6 = Vector.Multiply(acc6, new Vector<double>(array, i + vectorSize * 6));
    acc7 = Vector.Multiply(acc7, new Vector<double>(array, i + vectorSize * 7));

acc0 = Vector.Multiply(acc0, acc1);
acc2 = Vector.Multiply(acc2, acc3);
acc4 = Vector.Multiply(acc4, acc5);
acc6 = Vector.Multiply(acc6, acc7);
acc0 = Vector.Multiply(acc0, acc2);
acc4 = Vector.Multiply(acc4, acc6);
acc0 = Vector.Multiply(acc0, acc4);
// from here on it's the same
var tempArray = new double[Vector<double>.Count];
acc0.CopyTo(tempArray);
var result = tempArray.Aggregate(1d, (p, d) => p * d);
for (; i < array.Length; i++)
    result *= array[i];

这使得最后一个循环的运行时间可能是过去的 8 倍，这可以通过额外的单向量每次迭代循环来避免。

【讨论】：

谢谢哈罗德。你有我可以查看的使用 System.Runtime.Intrinsics.X86 API 的示例代码吗？

@little_stone_05 我加了一点

奇怪的是 JIT 引擎在 FP 操作之间使用 vpshufd。 （而不是vshufpd 或vunpckhpd）。我认为在这种情况下，某些 AMD CPU 可能会绕过转发。 （英特尔没有。）

@harold 你将如何将数组加载到 Vector256？

@little_stone_05 以LoadVector256 为例