使用system.numerics对数组元素进行c#乘法

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在评论中提到的，通常你会将一个256位的向量拆分为两个128位的部分，并将它们相乘，然后使用shuffles来处理128位的部分。

通常的方法可以用于需要.NET Core3.0或更高版本的System.Runtime.Intrinsics.X86 API。

例如：

static double product(Vector256<double> vec)
{
    var t = Sse2.Multiply(vec.GetLower(), vec.GetUpper());
    return t.GetElement(0) * t.GetElement(1);
}

这看起来可能很糟糕，给JIT引擎留下了一个神秘的GetElement来解决，但实际上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的预期要好。我以非内联的方式测试了这个函数，通常它应该是内联的，负载应该会消失。

主循环可以改进，因为乘法的吞吐量比它的延迟要好得多。现在，乘法一次完成一个(即一次完成一个向量乘法)，其间有延迟( Haswell上有5个周期，Broadwell上有4个周期及更高版本)，以等待上一个乘法完成，但例如，Intel 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倍，这可以通过使用额外的每次迭代单个向量的循环来避免。