Как попросить GCC полностью развернуть этот цикл (т.е. очистить этот цикл)?

Наивная функция для достижения этого -

void slow_4(double       *c,
            const double *a,
            const double *b,
            size_t        n)
{
    size_t row, col, i;

    for (row = 0; row < 4; row++)
        for (col = 0; col < 4; col++)
            for (i = 0; i < n; i++)
                c[4*row+col] += a[4*i+row] * b[4*i+col];
}

GCC генерирует довольно хороший код для SSE2 / SSE3, используя

#if defined(__SSE2__) || defined(__SSE3__)

typedef  double  vec2d  __attribute__((vector_size (2 * sizeof (double))));

void fast_4(vec2d       *c,
            const vec2d *a,
            const vec2d *b,
            size_t       n)
{
    const vec2d *const b_end = b + 2L * n;

    vec2d s00 = c[0];
    vec2d s02 = c[1];
    vec2d s10 = c[2];
    vec2d s12 = c[3];
    vec2d s20 = c[4];
    vec2d s22 = c[5];
    vec2d s30 = c[6];
    vec2d s32 = c[7];

    while (b < b_end) {
        const vec2d b0 = b[0];
        const vec2d b2 = b[1];
        const vec2d a0 = { a[0][0], a[0][0] };
        const vec2d a1 = { a[0][1], a[0][1] };
        const vec2d a2 = { a[1][0], a[1][0] };
        const vec2d a3 = { a[1][1], a[1][1] };
        s00 += a0 * b0;
        s10 += a1 * b0;
        s20 += a2 * b0;
        s30 += a3 * b0;
        s02 += a0 * b2;
        s12 += a1 * b2;
        s22 += a2 * b2;
        s32 += a3 * b2;
        b += 2;
        a += 2;
    }

    c[0] = s00;
    c[1] = s02;
    c[2] = s10;
    c[3] = s12;
    c[4] = s20;
    c[5] = s22;
    c[6] = s30;
    c[7] = s32; 
}

#endif

Для AVX GCC может работать еще лучше с

#if defined(__AVX__) || defined(__AVX2__)

typedef  double  vec4d  __attribute__((vector_size (4 * sizeof (double))));

void fast_4(vec4d       *c,
            const vec4d *a,
            const vec4d *b,
            size_t       n)
{
    const vec4d *const b_end = b + n;

    vec4d s0 = c[0];
    vec4d s1 = c[1];
    vec4d s2 = c[2];
    vec4d s3 = c[3];

    while (b < b_end) {
        const vec4d bc = *(b++);
        const vec4d ac = *(a++);
        const vec4d a0 = { ac[0], ac[0], ac[0], ac[0] };
        const vec4d a1 = { ac[1], ac[1], ac[1], ac[1] };
        const vec4d a2 = { ac[2], ac[2], ac[2], ac[2] };
        const vec4d a3 = { ac[3], ac[3], ac[3], ac[3] };
        s0 += a0 * bc;
        s1 += a1 * bc;
        s2 += a2 * bc;
        s3 += a3 * bc;
    }

    c[0] = s0;
    c[1] = s1;
    c[2] = s2;
    c[3] = s3;
}

#endif

Версия SSE3 сгенерированной сборки с использованием gcc-4.8.4 (-O2 -march=x86-64 -mtune=generic -msse3) по сути

fast_4:
        salq    $5, %rcx
        movapd  (%rdi), %xmm13
        addq    %rdx, %rcx
        cmpq    %rcx, %rdx
        movapd  16(%rdi), %xmm12
        movapd  32(%rdi), %xmm11
        movapd  48(%rdi), %xmm10
        movapd  64(%rdi), %xmm9
        movapd  80(%rdi), %xmm8
        movapd  96(%rdi), %xmm7
        movapd  112(%rdi), %xmm6
        jnb     .L2
.L3:
        movddup (%rsi), %xmm5
        addq    $32, %rdx
        movapd  -32(%rdx), %xmm1
        addq    $32, %rsi
        movddup -24(%rsi), %xmm4
        movapd  %xmm5, %xmm14
        movddup -16(%rsi), %xmm3
        movddup -8(%rsi), %xmm2
        mulpd   %xmm1, %xmm14
        movapd  -16(%rdx), %xmm0
        cmpq    %rdx, %rcx
        mulpd   %xmm0, %xmm5
        addpd   %xmm14, %xmm13
        movapd  %xmm4, %xmm14
        mulpd   %xmm0, %xmm4
        addpd   %xmm5, %xmm12
        mulpd   %xmm1, %xmm14
        addpd   %xmm4, %xmm10
        addpd   %xmm14, %xmm11
        movapd  %xmm3, %xmm14
        mulpd   %xmm0, %xmm3
        mulpd   %xmm1, %xmm14
        mulpd   %xmm2, %xmm0
        addpd   %xmm3, %xmm8
        mulpd   %xmm2, %xmm1
        addpd   %xmm14, %xmm9
        addpd   %xmm0, %xmm6
        addpd   %xmm1, %xmm7
        ja      .L3
.L2:
        movapd  %xmm13, (%rdi)
        movapd  %xmm12, 16(%rdi)
        movapd  %xmm11, 32(%rdi)
        movapd  %xmm10, 48(%rdi)
        movapd  %xmm9, 64(%rdi)
        movapd  %xmm8, 80(%rdi)
        movapd  %xmm7, 96(%rdi)
        movapd  %xmm6, 112(%rdi)
        ret

Версия AVX созданной сборки (-O2 -march=x86-64 -mtune=generic -mavx) по существу

fast_4:
        salq       $5, %rcx
        vmovapd    (%rdi), %ymm5
        addq       %rdx, %rcx
        vmovapd    32(%rdi), %ymm4
        cmpq       %rcx, %rdx
        vmovapd    64(%rdi), %ymm3
        vmovapd    96(%rdi), %ymm2
        jnb        .L2
.L3:
        addq       $32, %rsi
        vmovapd    -32(%rsi), %ymm1
        addq       $32, %rdx
        vmovapd    -32(%rdx), %ymm0
        cmpq       %rdx, %rcx
        vpermilpd  $0, %ymm1, %ymm6
        vperm2f128 $0, %ymm6, %ymm6, %ymm6
        vmulpd     %ymm0, %ymm6, %ymm6
        vaddpd     %ymm6, %ymm5, %ymm5
        vpermilpd  $15, %ymm1, %ymm6
        vperm2f128 $0, %ymm6, %ymm6, %ymm6
        vmulpd     %ymm0, %ymm6, %ymm6
        vaddpd     %ymm6, %ymm4, %ymm4
        vpermilpd  $0, %ymm1, %ymm6
        vpermilpd  $15, %ymm1, %ymm1
        vperm2f128 $17, %ymm6, %ymm6, %ymm6
        vperm2f128 $17, %ymm1, %ymm1, %ymm1
        vmulpd     %ymm0, %ymm6, %ymm6
        vmulpd     %ymm0, %ymm1, %ymm0
        vaddpd     %ymm6, %ymm3, %ymm3
        vaddpd     %ymm0, %ymm2, %ymm2
        ja         .L3
.L2:
        vmovapd    %ymm5, (%rdi)
        vmovapd    %ymm4, 32(%rdi)
        vmovapd    %ymm3, 64(%rdi)
        vmovapd    %ymm2, 96(%rdi)
        vzeroupper
        ret

Планирование регистров, я полагаю, не оптимально, но и не выглядит ужасным. Я лично доволен вышеизложенным, не пытаясь вручную оптимизировать его на данном этапе.

На процессоре Core i5-4200U (с поддержкой AVX2) быстрые версии вышеуказанных функций вычисляют произведение двух матриц 4 × 256 за 1843 цикла ЦП (медиана) для SSE3 и 1248 циклов для AVX2. Это составляет 1,8 и 1,22 цикла на запись в матрице. Для сравнения, невекторизованная медленная версия занимает около 11 циклов на запись в матрице.

(Счетчики циклов являются средними значениями, то есть половина тестов была быстрее. Я провел лишь приблизительный тест с ~ 100 тыс. Повторений или около того, так что относитесь к этим числам с недоверием.)

На этом процессоре эффекты кеширования таковы, что AVX2 при размере матрицы 4 × 512 по-прежнему составляет 1,2 цикла на запись, но при 4 × 1024 он падает до 1,4, при 4 × 4096 до 1,5, при 4 × 8192 до 1,8 и при 4 × 65536 до 2,2 цикла на запись. Версия SSE3 остается на 1,8 цикла на запись до 4 × 3072, после чего она начинает замедляться; при 4 × 65536 это тоже примерно 2,2 цикла на вход. Я действительно считаю, что этот (ноутбук!) CPU имеет ограниченную пропускную способность кэша на данный момент.