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

Question

我正在尝试使用 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;
        }

Answer 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 倍，这可以通过额外的单向量每次迭代循环来避免。