.NET BLAS库开发中FMA指令性能低于AVX的问题排查

核心实现代码 #

using System;
using System.Runtime.Intrinsics;
using System.Runtime.Intrinsics.X86;

namespace LinearAlgebra
{
    public static class Level1
    {
        const int Vector256Length = 8;
        const int Vector256Block0 = Vector256Length * 0;
        const int Vector256Block1 = Vector256Length * 1;
        const int Vector256Block2 = Vector256Length * 2;

        public static unsafe void ApplyAffineTransformFma(int N, float* X, float* Y, float m11, float m12, float m21, float m22, float dx, float dy)
        {
            var vm11 = Avx.BroadcastScalarToVector256(&m11);
            var vm12 = Avx.BroadcastScalarToVector256(&m12);
            var vm21 = Avx.BroadcastScalarToVector256(&m21);
            var vm22 = Avx.BroadcastScalarToVector256(&m22);
            var vdx = Avx.BroadcastScalarToVector256(&dx);
            var vdy = Avx.BroadcastScalarToVector256(&dy);

            int n = 0;
            for (; n <= N - Vector256Block2; n += Vector256Block2)
            {
                float* xn = X + n;
                float* yn = Y + n;

                var vx0 = Avx.LoadVector256(xn + Vector256Block0);
                var vx1 = Avx.LoadVector256(xn + Vector256Block1);
                var vy0 = Avx.LoadVector256(yn + Vector256Block0);
                var vy1 = Avx.LoadVector256(yn + Vector256Block1);

                var sx0 = Fma.MultiplyAdd(vm21, vy0, vdx);
                var sx1 = Fma.MultiplyAdd(vm21, vy1, vdx);
                var sy0 = Fma.MultiplyAdd(vm12, vx0, vdy);
                var sy1 = Fma.MultiplyAdd(vm12, vx1, vdy);

                vx0 = Fma.MultiplyAdd(vm11, vx0, sx0);
                vx1 = Fma.MultiplyAdd(vm11, vx1, sx1);
                vy0 = Fma.MultiplyAdd(vm22, vy0, sy0);
                vy1 = Fma.MultiplyAdd(vm22, vy1, sy1);

                Avx.Store(xn + Vector256Block0, vx0);
                Avx.Store(xn + Vector256Block1, vx1);
                Avx.Store(yn + Vector256Block0, vy0);
                Avx.Store(yn + Vector256Block1, vy1);
            }

            for (; n <= N - Vector256Length; n += Vector256Length)
            {
                float* xn = X + n;
                float* yn = Y + n;

                var vx = Avx.LoadVector256(xn);
                var vy = Avx.LoadVector256(yn);

                var sx = Fma.MultiplyAdd(vm21, vy, vdx);
                var sy = Fma.MultiplyAdd(vm12, vx, vdy);

                vx = Fma.MultiplyAdd(vm11, vx, sx);
                vy = Fma.MultiplyAdd(vm22, vy, sy);

                Avx.Store(xn, vx);
                Avx.Store(yn, vy);
            }

            ApplyAffineTransformOdd(N - n, X + n, Y + n, m11, m12, m21, m22, dx, dy);
        }

        public static unsafe void ApplyAffineTransformAvx(int N, float* X, float* Y, float m11, float m12, float m21, float m22, float dx, float dy)
        {
            var vm11 = Avx.BroadcastScalarToVector256(&m11);
            var vm12 = Avx.BroadcastScalarToVector256(&m12);
            var vm21 = Avx.BroadcastScalarToVector256(&m21);
            var vm22 = Avx.BroadcastScalarToVector256(&m22);
            var vdx = Avx.BroadcastScalarToVector256(&dx);
            var vdy = Avx.BroadcastScalarToVector256(&dy);

            int n = 0;
            for (; n <= N - Vector256Length; n += Vector256Length)
            {
                float* xn = X + n;
                float* yn = Y + n;

                var vx = Avx.LoadVector256(xn);
                var vy = Avx.LoadVector256(yn);

                var x1 = Avx.Multiply(vm11, vx);
                var x2 = Avx.Multiply(vm21, vy);
                var y1 = Avx.Multiply(vm12, vx);
                var y2 = Avx.Multiply(vm22, vy);

                x1 = Avx.Add(x1, x2);
                y1 = Avx.Add(y1, y2);

                vx = Avx.Add(x1, vdx);
                vy = Avx.Add(y1, vdy);

                Avx.Store(xn, vx);
                Avx.Store(yn, vy);
            }

            ApplyAffineTransformOdd(N - n, X + n, Y + n, m11, m12, m21, m22, dx, dy);
        }

        public static unsafe void ApplyAffineTransformOdd(int N, float* X, float* Y, float m11, float m12, float m21, float m22, float dx, float dy)
        {
            for (int n = 0; n < N; n++)
            {
                float xn = X[n];
                float yn = Y[n];

                X[n] = m11 * xn + m21 * yn + dx;
                Y[n] = m12 * xn + m22 * yn + dy;
            }
        }
    }
}

初始测试环境与结果 #

测试运行在Intel Core i5-7200U（Skylake架构）设备，代码采用.NET 5 x64模式编译，初始基准测试结果如下：

基准测试代码 #

using BenchmarkDotNet.Attributes;
using LinearAlgebra;
using LinearAlgebra.Tests;

namespace Benchmarks.LinearAlgebra
{
    [MinColumn, MaxColumn, MedianColumn]
    public unsafe class TransformBenchmark
    {
        [Params(68, 1159, 4101, 16389, 32773)]
        public int N { get; set; }

        public float[] X;
        public float[] Y;

        const float m11 = 1f;
        const float m12 = 0f;
        const float m21 = 1f;
        const float m22 = 1f;
        const float dx = 0f;
        const float dy = 0f;

        [GlobalSetup]
        public void GlobalSetup()
        {
            X = LinAlgGenerator.GetRandomVector(N);
            Y = LinAlgGenerator.GetRandomVector(N);
        }

        [Benchmark(Baseline = true)]
        public void AffineTransformFma()
        {
            fixed (float* ptrX = X, ptrY = Y)
            {
                Level1.ApplyAffineTransformFma(N, ptrX, ptrY, m11, m12, m21, m22, dx, dy);
            }
        }
        
        [Benchmark]
        public void AffineTransformAvx()
        {
            fixed (float* ptrX = X, ptrY = Y)
            {
                Level1.ApplyAffineTransformAvx(N, ptrX, ptrY, m11, m12, m21, m22, dx, dy);
            }
        }
    }
}

问题分析 #

初始结果完全不符合理论预期：FMA代码性能几乎仅为AVX版本的一半。根据英特尔Intrinsics指南，Fma.MultiplyAdd延迟为4、吞吐量为0.5，理论上处理16个浮点数FMA仅需50周期，AVX则需要66周期，FMA理应更快。开发人员已确认代码256位寄存器占用未超过16个，关闭流水线测试结果仍一致，未找到问题根源。

排查结论 #

后续排查发现异常为BenchmarkDotNet的bug导致，使用自研基准测试代码多次运行得到如下结果：

修正测试后多数场景下FMA版本性能比AVX高30%，结果不受测试执行顺序影响。

内容的提问来源于stack exchange，提问作者Станислав Герасименко