如何使用 C# 中的矢量化计算双精度矩阵的总和？

首先，你应该做一些基准测试和/或分析，问问自己这是否真的重要？求和是一个非常简单的计算，200x200不是很大。我猜它可能会达到一微秒的量级，但这只是一个猜测。你还需要一个基准来决定你是否真的取得了任何改进，或者你是否只是无缘无故地使代码变得更加复杂。

但这真的是应用程序的最大瓶颈吗？优化通常是关于避免做工作，或避免重做工作。任何 SIMD 优化都能为您提供的最好的效果是持续的加速。浪费时间优化对用户没有明显影响的功能是没有意义的。

如果你决定你需要优化，那么我会从摆脱指数计算开始。当你这样做时，框架本质上是做一个 -calculation。这可能需要比实际值求和更长的时间。优化器可能会删除其中的一些内容，但是在未实际确认的情况下，我不会从优化器中假设任何内容。您可以使用 执行不安全的路由，也可以创建一个自定义矩阵类，该类使用 1D 数组进行存储，该数组仅允许您使用单个循环对值求和，如果您出于其他原因需要语法，请自行实现 2D 索引器。matrix[i, j]i * width + jfixed (double* ptr = matrix )[x, y]

如果你真的需要 SIMD 的性能，你可以采取两种方式

1. Vector<T>
2. 内部函数，如 Vector256

请参阅比较。简而言之，内部函数提供了更好的性能，但代价是将其绑定到特定的 CPU 平台。

无论哪种情况，您都需要了解内存布局才能正确加载元素。但是一旦完成，它应该非常简单，只需将所有向量相加，最后对元素求和即可。如果元素计数不能与向量长度相等，则最后可能会有一些标量代码。

这可以通过向量 API 很好地完成，至少在自由使用类的情况下是这样。System.NumericsUnsafe

据我所知，没有好的“标准”方法可以从 2D 矩阵加载向量。正常载荷的过载都不适用，并且没有正常的方法来获取 2D 数组。但是有了，我们无论如何都可以完成它。Span<T>Unsafe

使用带有 8 个独立累加器的 8 展开（请参阅使用多个累加器展开 FP 循环），并通过操作引用将 2D 矩阵视为 1D 数组，我们可以这样做：（未测试，但在 sharplab.io 上编译）Unsafe

static unsafe double Sum(double[,] matrix)
{
    Vector<double> sum0 = Vector<double>.Zero;
    Vector<double> sum1 = Vector<double>.Zero;
    Vector<double> sum2 = Vector<double>.Zero;
    Vector<double> sum3 = Vector<double>.Zero;
    Vector<double> sum4 = Vector<double>.Zero;
    Vector<double> sum5 = Vector<double>.Zero;
    Vector<double> sum6 = Vector<double>.Zero;
    Vector<double> sum7 = Vector<double>.Zero;
    double sum8 = 0;
    uint vlen = (uint)Vector<double>.Count;

    ref double unaligneddata = ref matrix[0, 0];
    uint i = 0;
    uint alignmask = vlen * sizeof(double) - 1;
    for (; i < matrix.Length && ((IntPtr)Unsafe.AsPointer(ref unaligneddata) & alignmask) != 0; i++)
    {
        sum8 += unaligneddata;
        unaligneddata = ref Unsafe.Add(ref unaligneddata, 1);
    }
    uint alignment_skipped = i;
    ref Vector<double> data = ref Unsafe.As<double, Vector<double>>(ref unaligneddata);
    uint bigChunk = ((uint)matrix.Length - alignment_skipped & (0u - (vlen * 8))) + alignment_skipped;
    for (; i < bigChunk; i += vlen * 8)
    {
        sum0 += data;
        sum1 += Unsafe.Add(ref data, 1);
        sum2 += Unsafe.Add(ref data, 2);
        sum3 += Unsafe.Add(ref data, 3);
        sum4 += Unsafe.Add(ref data, 4);
        sum5 += Unsafe.Add(ref data, 5);
        sum6 += Unsafe.Add(ref data, 6);
        sum7 += Unsafe.Add(ref data, 7);
        data = ref Unsafe.Add(ref data, 8);
    }
    uint smallChunk = ((uint)matrix.Length - alignment_skipped & (0u - vlen)) + alignment_skipped;
    for (; i < smallChunk; i += vlen)
    {
        sum0 += data;
        data = ref Unsafe.Add(ref data, 1);
    }
    ref double remainder = ref Unsafe.As<Vector<double>, double>(ref data);
    for (; i < matrix.Length; i++)
    {
        sum8 += remainder;
        remainder = ref Unsafe.Add(ref remainder, 1);
    }

    sum0 += sum1;
    sum2 += sum3;
    sum4 += sum5;
    sum6 += sum7;
    sum0 += sum2;
    sum4 += sum6;
    sum0 += sum4;
    return Vector.Dot(sum0, new Vector<double>(1.0)) + sum8;
}

最后用一个水平求和有点傻，但很短，而且只发生一次。Vector.Dot

开始时尝试使地址对齐的循环主要用于不使用 AVX 时。不幸的是，据我所知，这需要（关键字，而不是类），即使原始指针立即转换为整数并且从未用作指针。unsafe

当 AVX2 可用时（在没有 AVX2 的情况下为 128 位，即使您只使用 float/double），主循环在程序集中可能如下所示：Vector<T>

L008c: vaddpd ymm0, ymm0, [rax]
L0091: vaddpd ymm1, ymm1, [rax+0x20]
L0097: vaddpd ymm2, ymm2, [rax+0x40]
L009d: vaddpd ymm3, ymm3, [rax+0x60]
L00a3: vaddpd ymm4, ymm4, [rax+0x80]
L00ac: vaddpd ymm5, ymm5, [rax+0xa0]
L00b5: vaddpd ymm6, ymm6, [rax+0xc0]
L00be: vaddpd ymm7, ymm7, [rax+0xe0]
L00c7: add rax, 0x100
L00cd: add r8d, 0x20
L00d1: cmp r8d, ecx
L00d4: jb short L008c

对我来说看起来不错。我们可以通过直接比较地址而不是保留冗余索引来保存此处，但这没什么大不了的。add