4
\$\begingroup\$

Benchmark

An optimised 4x4 double precision matrix multiply using intel AVX intrinsics. Two different variations.

Gist

For quick benchmark (with a compatible system) copy paste the command below. Runs tests on clang and gcc on optimisation levels 0 -> 3. Runs a naive matrix multiplication NORMAL as a reference.

curl -OL https://gist.githubusercontent.com/juliuskoskela/2381831c86041eb2ebf271011db7b2bf/raw/cebbcec5f2f02ac75c36e9346776cd97bf1e68a5/run.sh
bash run.sh >> M4X4D_BENCHMARK_RESULTS.txt
cat M4X4D_TEST_RESULTS.txt

Quick result Linux AMD Ryzen 3 3100 4-Core:

(best)NORMAL
gcc -O1 -march=native -mavx
time: 0.168030

(best)AVX MUL + FMADD
clang -O1 -march=native -mavx
time: 0.033416

Results indicate that there is some discrepancy in optimisation between gcc and clang and on Linux and Mac OS. On Mac OS neither compiler is able to properly optimise the normal version. On Linux clang still didn't know how to optimise further, but gcc Managed to find the same speeds as we hit with the intrinsics version BUT only on -O1 and -O2 and NOT on -O3.

Code

Version of the matrix multiplication which uses one multiply and one multiply + add instruction in the inner loop. Inner loop is unrolled and with comments.

typedef union u_m4d
{
    __m256d m256d[4];
    double  d[4][4];
}   t_m4d;

// Example matrices.
//
// left[0] (1 2 3 4)  | right[0] (4 1 1 1)
// left[1] (2 2 3 4)  | right[1] (3 2 2 2)
// left[2] (3 2 3 4)  | right[2] (2 3 3 3)
// left[3] (4 2 3 4)  | right[3] (1 4 4 4)

static inline void m4x4d_avx_mul(
        t_m4d *restrict dst,
        const t_m4d *restrict left,
        const t_m4d *restrict right)
{
    __m256d ymm0;
    __m256d ymm1;
    __m256d ymm2;
    __m256d ymm3;

    // Fill registers ymm0 -> ymm3 with a single value
    // from the i:th column of the left
    // hand matrix.
    //
    // left[0] (1 2 3 4) -> ymm0 (1 1 1 1)
    // left[1] (2 2 3 4) -> ymm1 (2 2 2 2)
    // left[2] (3 2 3 4) -> ymm2 (3 3 3 3)
    // left[3] (4 2 3 4) -> ymm3 (4 4 4 4)

    ymm0 = _mm256_broadcast_sd(&left->d[0][0]);
    ymm1 = _mm256_broadcast_sd(&left->d[0][1]);
    ymm2 = _mm256_broadcast_sd(&left->d[0][2]);
    ymm3 = _mm256_broadcast_sd(&left->d[0][3]);

    // Multiply vector at register ymm0 with right row[0]
    //
    // 1  1  1  1   <- ymm0
    // *
    // 4  2  3  4   <- right[0]
    // ----------
    // 4  2  3  4   <- ymm0

    ymm0 = _mm256_mul_pd(ymm0, right->m256d[0]);

    // Multiply vector at register ymm1 with right hand
    // row[1] and add at each multiply add the corresponding
    // value at ymm0 tp the result.
    //
    // 2  2  2  2   <- ymm1
    // *
    // 3  2  3  4   <- right[1]
    // +
    // 4  2  3  4   <- ymm0
    // ----------
    // 10 6  9  12  <- ymm0

    ymm0 = _mm256_fmadd_pd(ymm1, right->m256d[1], ymm0);

    // We repeat for ymm2 -> ymm3.
    //
    // 3  3  3  3   <- ymm2
    // *
    // 2  2  3  4   <- right[2]
    // ----------
    // 6  6  9  12  <- ymm2
    //
    // 2  2  2  2   <- ymm3
    // *
    // 3  2  3  4   <- right[3]
    // +
    // 6  6  9  12  <- ymm2
    // ----------
    // 10 14 21 28  <- ymm2

    ymm2 = _mm256_mul_pd(ymm2, right->m256d[2]);
    ymm2 = _mm256_fmadd_pd(ymm3, right->m256d[3], ymm2);

    // Sum accumulated vectors at ymm0 and ymm2.
    //
    // 10  6   9   12   <- ymm0
    // +
    // 10  14  21  28   <- ymm2
    // ----------
    // 20  20  30  40   <- dst[0] First row!

    dst->m256d[0] = _mm256_add_pd(ymm0, ymm2);

    // Calculate dst[1]
    ymm0 = _mm256_broadcast_sd(&left->d[1][0]);
    ymm1 = _mm256_broadcast_sd(&left->d[1][1]);
    ymm2 = _mm256_broadcast_sd(&left->d[1][2]);
    ymm3 = _mm256_broadcast_sd(&left->d[1][3]);
    ymm0 = _mm256_mul_pd(ymm0, right->m256d[0]);
    ymm0 = _mm256_fmadd_pd(ymm1, right->m256d[1], ymm0);
    ymm2 = _mm256_mul_pd(ymm2, right->m256d[2]);
    ymm2 = _mm256_fmadd_pd(ymm3, right->m256d[3], ymm2);
    dst->m256d[1] = _mm256_add_pd(ymm0, ymm2);

    // Calculate dst[2]
    ymm0 = _mm256_broadcast_sd(&left->d[2][0]);
    ymm1 = _mm256_broadcast_sd(&left->d[2][1]);
    ymm2 = _mm256_broadcast_sd(&left->d[2][2]);
    ymm3 = _mm256_broadcast_sd(&left->d[2][3]);
    ymm0 = _mm256_mul_pd(ymm0, right->m256d[0]);
    ymm0 = _mm256_fmadd_pd(ymm1, right->m256d[1], ymm0);
    ymm2 = _mm256_mul_pd(ymm2, right->m256d[2]);
    ymm2 = _mm256_fmadd_pd(ymm3, right->m256d[3], ymm2);
    dst->m256d[2] = _mm256_add_pd(ymm0, ymm2);

    // Calculate dst[3]
    ymm0 = _mm256_broadcast_sd(&left->d[3][0]);
    ymm1 = _mm256_broadcast_sd(&left->d[3][1]);
    ymm2 = _mm256_broadcast_sd(&left->d[3][2]);
    ymm3 = _mm256_broadcast_sd(&left->d[3][3]);
    ymm0 = _mm256_mul_pd(ymm0, right->m256d[0]);
    ymm0 = _mm256_fmadd_pd(ymm1, right->m256d[1], ymm0);
    ymm2 = _mm256_mul_pd(ymm2, right->m256d[2]);
    ymm2 = _mm256_fmadd_pd(ymm3, right->m256d[3], ymm2);
    dst->m256d[3] = _mm256_add_pd(ymm0, ymm2);
}

Version of the matrix multiplication which uses two multiplies and and add in the inner loop.

static inline void m4x4d_avx_mul2(
        t_m4d *restrict dst,
        const t_m4d *restrict left,
        const t_m4d *restrict right)
{
    __m256d ymm[4];

    for (int i = 0; i < 4; i++)
    {
        ymm[0] = _mm256_broadcast_sd(&left->d[i][0]);
        ymm[1] = _mm256_broadcast_sd(&left->d[i][1]);
        ymm[2] = _mm256_broadcast_sd(&left->d[i][2]);
        ymm[3] = _mm256_broadcast_sd(&left->d[i][3]);
        ymm[0] = _mm256_mul_pd(ymm[0], right->m256d[0]);
        ymm[1] = _mm256_mul_pd(ymm[1], right->m256d[1]);
        ymm[0] = _mm256_add_pd(ymm[0], ymm[1]);
        ymm[2] = _mm256_mul_pd(ymm[2], right->m256d[2]);
        ymm[3] = _mm256_mul_pd(ymm[3], right->m256d[3]);
        ymm[2] = _mm256_add_pd(ymm[2], ymm[3]);
        dst->m256d[i] = _mm256_add_pd(ymm[0], ymm[2]);
    }
}

Comparison matrix multiply that doesn't use intrinsics.

static inline void m4x4d_mul(double d[4][4], double l[4][4], double r[4][4])
{
    d[0][0] = l[0][0] * r[0][0] + l[0][1] * r[1][0] + l[0][2] * r[2][0] + l[0][3] * r[3][0];
    d[0][1] = l[0][0] * r[0][1] + l[0][1] * r[1][1] + l[0][2] * r[2][1] + l[0][3] * r[3][1];
    d[0][2] = l[0][0] * r[0][2] + l[0][1] * r[1][2] + l[0][2] * r[2][2] + l[0][3] * r[3][2];
    d[0][3] = l[0][0] * r[0][3] + l[0][1] * r[1][3] + l[0][2] * r[2][3] + l[0][3] * r[3][3];
    d[1][0] = l[1][0] * r[0][0] + l[1][1] * r[1][0] + l[1][2] * r[2][0] + l[1][3] * r[3][0];
    d[1][1] = l[1][0] * r[0][1] + l[1][1] * r[1][1] + l[1][2] * r[2][1] + l[1][3] * r[3][1];
    d[1][2] = l[1][0] * r[0][2] + l[1][1] * r[1][2] + l[1][2] * r[2][2] + l[1][3] * r[3][2];
    d[1][3] = l[1][0] * r[0][3] + l[1][1] * r[1][3] + l[1][2] * r[2][3] + l[1][3] * r[3][3];
    d[2][0] = l[2][0] * r[0][0] + l[2][1] * r[1][0] + l[2][2] * r[2][0] + l[2][3] * r[3][0];
    d[2][1] = l[2][0] * r[0][1] + l[2][1] * r[1][1] + l[2][2] * r[2][1] + l[2][3] * r[3][1];
    d[2][2] = l[2][0] * r[0][2] + l[2][1] * r[1][2] + l[2][2] * r[2][2] + l[2][3] * r[3][2];
    d[2][3] = l[2][0] * r[0][3] + l[2][1] * r[1][3] + l[2][2] * r[2][3] + l[2][3] * r[3][3];
    d[3][0] = l[3][0] * r[0][0] + l[3][1] * r[1][0] + l[3][2] * r[2][0] + l[3][3] * r[3][0];
    d[3][1] = l[3][0] * r[0][1] + l[3][1] * r[1][1] + l[3][2] * r[2][1] + l[3][3] * r[3][1];
    d[3][2] = l[3][0] * r[0][2] + l[3][1] * r[1][2] + l[3][2] * r[2][2] + l[3][3] * r[3][2];
    d[3][3] = l[3][0] * r[0][3] + l[3][1] * r[1][3] + l[3][2] * r[2][3] + l[3][3] * r[3][3];
};

Main method and utils

///////////////////////////////////////////////////////////////////////////////
//
// Main and utils for testing.

t_v4d   v4d_set(double n0, double n1, double n2, double n3)
{
    t_v4d   v;

    v.d[0] = n0;
    v.d[1] = n1;
    v.d[2] = n2;
    v.d[3] = n3;
    return (v);
}

t_m4d   m4d_set(t_v4d v0, t_v4d v1, t_v4d v2, t_v4d v3)
{
    t_m4d   m;

    m.m256d[0] = v0.m256d;
    m.m256d[1] = v1.m256d;
    m.m256d[2] = v2.m256d;
    m.m256d[3] = v3.m256d;
    return (m);
}

int main(int argc, char **argv)
{
    t_m4d   left;
    t_m4d   right;
    t_m4d   res;
    t_m4d   ctr;

    if (argc != 2)
        return (printf("usage: avx4x4 [iters]"));

    left = m4d_set(
        v4d_set(1, 2, 3, 4),
        v4d_set(2, 2, 3, 4),
        v4d_set(3, 2, 3, 4),
        v4d_set(4, 2, 3, 4));

    right = m4d_set(
        v4d_set(4, 2, 3, 4),
        v4d_set(3, 2, 3, 4),
        v4d_set(2, 2, 3, 4),
        v4d_set(1, 2, 3, 4));

    size_t  iters;
    clock_t begin;
    clock_t end;
    double  time_spent;

    // Test 1
    m4x4d_mul(ctr.d, left.d, right.d);
    iters = atoi(argv[1]);

    begin = clock();
    for (size_t i = 0; i < iters; i++)
    {
        m4x4d_mul(res.d, left.d, right.d);
        
        // To prevent loop unrolling with optimisation flags.
        __asm__ volatile ("" : "+g" (i));
    }
    end = clock();

    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    printf("\nNORMAL\n\ntime: %lf\n", time_spent);

    // Test 2
    m4x4d_avx_mul(&ctr, &left, &right);
    iters = atoi(argv[1]);

    begin = clock();
    for (size_t i = 0; i < iters; i++)
    {
        m4x4d_avx_mul(&res, &left, &right);
        __asm__ volatile ("" : "+g" (i));
    }
    end = clock();

    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    printf("\nAVX MUL + FMADD\n\ntime: %lf\n", time_spent);

    // Test 3
    m4x4d_avx_mul2(&ctr, &left, &right);
    iters = atoi(argv[1]);

    begin = clock();
    for (size_t i = 0; i < iters; i++)
    {
        m4x4d_avx_mul2(&res, &left, &right);
        __asm__ volatile ("" : "+g" (i));
    }
    end = clock();

    time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
    printf("\nAVX MUL + MUL + ADD\n\ntime: %lf\n", time_spent);
}
```
\$\endgroup\$
2
  • \$\begingroup\$ You've clearly done your homework, but can you explain a bit about the expected usage patterns? For example, it makes a difference whether the operands/result are always expected to be in memory (not from the point of view of C, but in reality), or might have preferred to be in registers (eg when "chaining" multiplications). It's fine if it's supposed to be general purpose too, but if it isn't then I can avoid making recommendations that are inappropriate for your use case. \$\endgroup\$
    – harold
    Jul 26 at 8:20
  • 1
    \$\begingroup\$ I don’t have a clear use case for this. The project was born out of just curiosity. I was thinking it could be one step in calculating a bigger matrix for example, possibly in a multithreaded way. Any speculation as to the strenghts and weaknessess of the approach is welcome. I should also update the benchmarks. Need to tag the functions with __attribute__((noinline)) to get proper results. When the compiler isn’t able to inline test loop I get a more reasonable result with the AVX methdo being clearly superior. \$\endgroup\$
    – jkoskela
    Jul 26 at 8:37
3
\$\begingroup\$

It looks pretty good to me. I still have some notes to make.

Let's start with, there are 5 arithmetic operations (per row) to do the work that 4 FMAs could do. How good or bad is that? Well basically not bad at all in the context of just a simple 4x4 multiplication, because there are also 8 loads for those 5 arithmetic operations, so the problem is the loads. OTOH there is also no advantage of "splitting" the calculation in that case, because having the FMAs of a single row all in a dependent chain wouldn't matter in that context. That picture can change though. For example in next example with 3 matrices, the second half of the code still has 5 arithmetic operations per row but only 4 loads in the good case. But chaining also means that the latency of calculating a row matters, and therefore splitting up the dependency chain into more of a tree the way you did becomes increasingly more interesting as you chain more multiplications .. but how interesting exactly, given that the calculations for different rows are independent, I couldn't say. So how does it balance? Who knows, try it.

Here's an odd thing about this code: let's compare which order to multiply 3 matrices in, (AB)C or A(BC)

Like (AB)C:

void m4x4d_avx_mul3(
        t_m4d *restrict dst,
        const t_m4d *restrict A,
        const t_m4d *restrict B,
        const t_m4d *restrict C)
{
    t_m4d tmp;
    m4x4d_avx_mul(&tmp, A, B);
    m4x4d_avx_mul(dst, &tmp, C);
}

(on godbolt)

Or like A(BC):

void m4x4d_avx_mul3(
        t_m4d *restrict dst,
        const t_m4d *restrict A,
        const t_m4d *restrict B,
        const t_m4d *restrict C)
{
    t_m4d tmp;
    m4x4d_avx_mul(&tmp, B, C);
    m4x4d_avx_mul(dst, A, &tmp);
}

(on godbolt)

Algebraically that's the same thing, but perhaps surprisingly, the second way works out better: that way, tmp is not actually written to scratch space on the stack, it is created in vector registers and those registers are then used as input for the second multiply. That cannot happen in the first implementation, because then tmp's individual entries are broadcasted from. OK actually it could happen, but GCC didn't do it and it would have to use shuffles to make it happen, which aren't that cheap.

It's not necessarily bad that it works out like this, but something to keep in mind.

You've mentioned in the comments that you could use this as part of multiplying larger matrices. That is possible, but it wouldn't be as good as it could be, by a significant margin. The problem with that is that 4x4 is too small, which causes two issues. First, it makes the ratio of loads to arithmetic too high. There are 8 loads to 4 FMAs here, but the ratio should be 1:1 or less, to avoid bottlenecking on loads (most modern hardware can execute two of each per cycle, except that loads tend to be slower on average due to various factors so you should even below 1:1). For actual 4x4 matrices there is no hope of achieving that, but for larger matrices there is, namely by using a bigger "small matrix".

For example using a 6x8 matrix as the base (6 rows deep, 2 vectors wide) would give you an inner loop with:

  • 2 vector loads
  • 6 broadcasted scalar loads
  • 12 FMAs

That still does not exceed the vector register budget, which is 16 on x64. Unfortunately a nice round 8x8 would exceed the register budget and introduce spilling, which rules out that strategy.

By the way I counted only FMAs and no muls, not only because they are equivalent in cost anyway, but also because if you add tiling (which is essentially mandatory for multiplying big matrices) then you will be starting with a non-zero destination matrix and adding more products into it.

Another issue with basing large matrix multiplication on the multiplication of 4x4 matrices is that it would either (when multiplying a bunch of independent 4x4 matrices, ie the products aren't summed) force a lot of stores (which is not great) or (when accumulating several products before storing them) it would run into the high latency-throughput product of FMA on typical hardware (commonly ranging from 8 to even 10). If the latency-throughput product is 8 (for example a latency of 4 and a throughput of 2/cycle) then that means there must be at least 8 independent FMAs simultaneously in the pipeline to keep it full, and a 4x4 matrix normally only has 4 independent FMAs at any time. It does have 8 if you decouple it like you did, but that also costs more work, so at a large scale it's better to get those independent FMAs by Going Big and not doing that decoupling.

\$\endgroup\$
1
  • \$\begingroup\$ Thanks for the great insights! Good point that I should look into rectangular divisions of the matrix as well when scaling the algorithm. Btw. I updated the tests so that the functions wouldn't be inlined into the main loop on higher O levels. Now I'm seeing a more consistent gain across the board and also the MUL + FMADD edging slightly ahead of MUL + MUL + ADD. \$\endgroup\$
    – jkoskela
    Jul 26 at 10:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.