Visual Studio Code – How to calculate 2×2 matrix multiplied by 2D vector using SSE intrinsics (32 bit floating points)? (C++, Mac and Windows)

It’s much better to let compilers vectorize over the whole arrays of LLLL and RRRR samples, not for one left, right sample pair at once.

With the same mixing matrix for a whole array of audio samples, you get nice asm with no shuffles. Borrowing code from Nole’s answer just to illustrate the auto-vectorization (you might want to simplify the

struct Vector {
    std::array<float,2> coef;
};

struct Matrix {
    std::array<float,4> coef;
};

    static Vector multiply( const Matrix& m, const Vector& v ) {
        Vector r;
        r.coef[0]  = v.coef[0]*m.coef[0] + v.coef[1]*m.coef[2];
        r.coef[1] = v.coef[0]*m.coef[1] + v.coef[1]*m.coef[3];
        return r;
    }


// The per-element functions need to inline into the loop,
// so target attributes need to match or be a superset.  
// Or better, just don't use target options on the per-sample function
__attribute__ ((target ("avx,fma")))
void intermix(float *__restrict left, float *__restrict right, const Matrix &m)
{
    for (int i=0 ; i<10240 ; i++){
        Vector v = {left[i], right[i]};
        v = multiply(m, v);
        left[i] = v.coef[0];
        right[i] = v.coef[1];
    }
}

GCC -O3 (without any target options) compiles this to nice AVX1 + FMA code, as per the __attribute__((target("avx,fma"))) (Similar to -march=x86-64-v3). (Godbolt)

# GCC (trunk) -O3
intermix(float*, float*, Matrix const&):
        vbroadcastss    (%rdx), %ymm5
        vbroadcastss    8(%rdx), %ymm4
        xorl    %eax, %eax
        vbroadcastss    4(%rdx), %ymm3
        vbroadcastss    12(%rdx), %ymm2       # broadcast each matrix element separately
.L2:
        vmulps  (%rsi,%rax), %ymm4, %ymm1     # a whole vector of 8 R*B
        vmulps  (%rsi,%rax), %ymm2, %ymm0
        vfmadd231ps     (%rdi,%rax), %ymm5, %ymm1
        vfmadd231ps     (%rdi,%rax), %ymm3, %ymm0
        vmovups %ymm1, (%rdi,%rax)
        vmovups %ymm0, (%rsi,%rax)
        addq    $32, %rax
        cmpq    $40960, %rax
        jne     .L2
        vzeroupper
        ret
main:
        movl    $13, %eax
        ret

Note how there are zero shuffle instructions because the matrix coefficients each get broadcast to a separate vector, so one vmulps can do R*B for 8 R samples in parallel, and so on.

Unfortunately GCC and clang both use indexed addressing modes, so the memory source vmulps and vfma instructions un-laminate into 2 uops for the back-end on Intel CPUs. And the stores can’t use the port 7 AGU on HSW/SKL. -march=skylake or any other specific Intel SnB-family uarch doesn’t fix that for either of them. Clang unrolls by default, so the extra pointer increments to avoid indexed addressing modes would be amortized. (It would actually just be 1 extra add instruction, since we’re modifying L and R in-place. You could of course change the function to copy-and-mix.)

If the data is hot in L1d cache, it’ll bottleneck on the front-end rather than load+FP throughput, but it still comes relatively close to 2 loads and 2 FMAs per clock.

Hmm, GCC is saving instructions but costing extra loads by loading the same L and R data twice, as memory-source operands for vmulps and vfmadd...ps. With -march=skylake it doesn’t do that, instead using a separate vmovups (which has no problem with an indexed addressing mode, but the later store still does.)

I haven’t looked at tuning choices from other GCC versions.

# GCC (trunk) -O3 -march=skylake   (which implies -mtune=skylake)
.L2:
        vmovups (%rsi,%rax), %ymm1
        vmovups (%rdi,%rax), %ymm0     # separate loads
        vmulps  %ymm1, %ymm5, %ymm2    # FP instructions using only registers
        vmulps  %ymm1, %ymm3, %ymm1
        vfmadd231ps     %ymm0, %ymm6, %ymm2
        vfmadd132ps     %ymm4, %ymm1, %ymm0
        vmovups %ymm2, (%rdi,%rax)
        vmovups %ymm0, (%rsi,%rax)
        addq    $32, %rax
        cmpq    $40960, %rax
        jne     .L2

This is 10 uops, so can issue 2 cycles per iteration on Ice Lake, 2.5c on Skylake. On Ice Lake, it will sustain 2x 256-bit mul/FMA per clock cycle.

On Skylake, it doesn’t bottleneck on AGU throughput since it’s 4 uops for ports 2,3 every 2.5 cycles. So that’s fine. No need for indexed addressing modes.