# HPC kernel for DGEMM: compiler v.s. assembly

This is a correct version, for computing a small matrix multiplication: C += A * B, where C is 12 * 4 in non-transpose setting, A is 12 * k in non-transpose setting and B is k * 4 in transpose setting.

Computation is scheduled k number of matrix rank-1 update. 12 and 4 are choose to best use available register file. On x86-64 with AVX and FMA support, this rank-1 update should take exactly 16 YMM registers. Ideally, each vector loaded from A, i.e., A1_vec, A2_vec, A3_vec below, will be reused 4 times, while each vector loaded from B, i.e., B1_vec, B2_vec, B3_vec, B4_vec are reused 3 times. All vectors C1_vec, ..., C12_vec from C are held in registers for k times's reuse. As a result, the ratio:

flops : data read = [2 * (12/4) * 4]/[(12/4) + 4] = 3.43


Plus the theoretical 4* speedup by AVX, performance should be 10-12 times higher than ordinary scalar implementation. In practice, k is a cache blocking factor. On most x86 with 32KB L1 Dcache, the optimal factor is 60.

Originally everything was written in C, using Intel's intrinsics from <immintrin.h>. But GCC compilation does not yield optimal register allocation. Typically I see that C11_vec and C12_vec are alternately load from and stored to memory, instead of being hold on register file all the time. I suspect this degrades performance. The use of FMA instruction might make compiler feel difficult to make optimal decision, as there are three possible choices, combined with possibility of direct memory load. As a result, I produce an inline assembly implementation that does what I see best.

As a comparison, original C code is posted as well. This is a good opportunity for benchmarking. It would be good if we can:

• see how good different compiler does on x86-64 architecture (otherwise AVX does not work, right?), by checking assembler output for register allocation.
• Benchmark with compiler.

Due to some misunderstanding I posted this question here far too earlier before it is ready. But never mind. Just take it as a special welcome to me. This site is brilliant, anyway.

void GEPDOT_AVX__12_4 (int k, double *A, double *B, double *C) {
double *A_end = A + 12 * k - 20;
asm volatile (
"vmovapd  (%[A]), %%ymm13\n\t"
"addq  $128, %[A]\n\t" "vmovapd -96(%[A]), %%ymm14\n\t" "vmovapd -64(%[A]), %%ymm13\n\t" "vmovapd -32(%[A]), %%ymm0\n\t" "vbroadcastsd (%[B]), %%ymm6\n\t" "addq$40, %[B]\n\t"
"vmulpd  %%ymm13, %%ymm6, %%ymm4\n\t"
"vmulpd  %%ymm14, %%ymm6, %%ymm5\n\t"
"vmulpd  %%ymm15, %%ymm6, %%ymm6\n\t"
"vmulpd  %%ymm13, %%ymm9, %%ymm7\n\t"
"vmulpd  %%ymm14, %%ymm9, %%ymm8\n\t"
"vmulpd  %%ymm15, %%ymm9, %%ymm9\n\t"
"vmulpd  %%ymm13, %%ymm12, %%ymm10\n\t"
"vmulpd  %%ymm14, %%ymm12, %%ymm11\n\t"
"vmulpd  %%ymm15, %%ymm12, %%ymm12\n\t"
"vmulpd  %%ymm2, %%ymm13, %%ymm13\n\t"
"vmulpd  %%ymm2, %%ymm14, %%ymm14\n\t"
"vmulpd  %%ymm2, %%ymm15, %%ymm15\n\t"
"cmpl  $1, %[k]\n\t" "je End\n" "K%=:\n\t" /* first iteration */ "vfmadd231pd %%ymm0, %%ymm3, %%ymm4\n\t" "vmovapd (%[A]), %%ymm1\n\t" "addq$192, %[A]\n\t"
"vmovapd  -160(%[A]), %%ymm2\n\t"
"addq  $64, %[B]\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm7\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm8\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm9\n\t" "vbroadcastsd -56(%[B]), %%ymm3\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm10\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm11\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm12\n\t" "vbroadcastsd -48(%[B]), %%ymm3\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm13\n\t" "vmovapd -128(%[A]), %%ymm0\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm14\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm15\n\t" "vbroadcastsd -40(%[B]), %%ymm3\n\t" /* second iteration */ "vfmadd231pd %%ymm0, %%ymm3, %%ymm4\n\t" "vmovapd -96(%[A]), %%ymm1\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm5\n\t" "vmovapd -64(%[A]), %%ymm2\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm6\n\t" "vbroadcastsd -32(%[B]), %%ymm3\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm7\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm8\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm9\n\t" "vbroadcastsd -24(%[B]), %%ymm3\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm10\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm11\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm12\n\t" "vbroadcastsd -16(%[B]), %%ymm3\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm13\n\t" "vmovapd -32(%[A]), %%ymm0\n\t" "cmpq %[A], %[A_end]\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm14\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm15\n\t" "vbroadcastsd -8(%[B]), %%ymm3\n\t" "jne K%=\n" "End:\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm4\n\t" "vmovapd (%[A]), %%ymm1\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm5\n\t" "vmovapd 32(%[A]), %%ymm2\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm6\n\t" "vbroadcastsd 8(%[B]), %%ymm3\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm7\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm8\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm9\n\t" "vbroadcastsd 16(%[B]), %%ymm3\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm10\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm11\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm12\n\t" "vbroadcastsd 24(%[B]), %%ymm3\n\t" "vfmadd231pd %%ymm0, %%ymm3, %%ymm13\n\t" "vfmadd231pd %%ymm1, %%ymm3, %%ymm14\n\t" "vfmadd231pd %%ymm2, %%ymm3, %%ymm15\n\t" /* write-back */ "vaddpd (%[C]), %%ymm4, %%ymm4\n\t" "vmovapd %%ymm4, (%[C])\n\t" "vaddpd 32(%[C]), %%ymm5, %%ymm5\n\t" "vmovapd %%ymm5, 32(%[C])\n\t" "vaddpd 64(%[C]), %%ymm6, %%ymm6\n\t" "vmovapd %%ymm6, 64(%[C])\n\t" "vaddpd 96(%[C]), %%ymm7, %%ymm7\n\t" "vmovapd %%ymm7, 96(%[C])\n\t" "vaddpd 128(%[C]), %%ymm8, %%ymm8\n\t" "vmovapd %%ymm8, 128(%[C])\n\t" "vaddpd 160(%[C]), %%ymm9, %%ymm9\n\t" "vmovapd %%ymm9, 160(%[C])\n\t" "vaddpd 192(%[C]), %%ymm10, %%ymm10\n\t" "vmovapd %%ymm10, 192(%[C])\n\t" "vaddpd 224(%[C]), %%ymm11, %%ymm11\n\t" "vmovapd %%ymm11, 224(%[C])\n\t" "vaddpd 256(%[C]), %%ymm12, %%ymm12\n\t" "vmovapd %%ymm12, 256(%[C])\n\t" "vaddpd 288(%[C]), %%ymm13, %%ymm13\n\t" "vmovapd %%ymm13, 288(%[C])\n\t" "vaddpd 320(%[C]), %%ymm14, %%ymm14\n\t" "vmovapd %%ymm14, 320(%[C])\n\t" "vaddpd 352(%[C]), %%ymm15, %%ymm15\n\t" "vmovapd %%ymm15, 352(%[C])\n\t" /* exit AVX mode */ "vzeroupper\n\t" : [A] "+r" (A), [B] "+r" (B), [C] "+r" (C) : [k] "r" (k), [A_end] "r" (A_end) : "ymm0", "ymm1", "ymm2", "ymm3", "ymm4", "ymm5", "ymm6", "ymm7", "ymm8", "ymm9", "ymm10", "ymm11", "ymm12", "ymm13", "ymm14", "ymm15", "memory"); }  C version #include <immintrin.h> void GEPDOT_AVX__12_4 (int k, double *A, double *B, double *C) { __m256d A1_vec = _mm256_load_pd(A); A += 4; __m256d B_vec = _mm256_broadcast_sd(B); B++; __m256d C1_vec = A1_vec * B_vec; __m256d A2_vec = _mm256_load_pd(A); A += 4; __m256d C2_vec = A2_vec * B_vec; __m256d A3_vec = _mm256_load_pd(A); A += 4; __m256d C3_vec = A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; __m256d C4_vec = A1_vec * B_vec; __m256d C5_vec = A2_vec * B_vec; __m256d C6_vec = A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; __m256d C7_vec = A1_vec * B_vec; __m256d C8_vec = A2_vec * B_vec; __m256d C9_vec = A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; __m256d C10_vec = A1_vec * B_vec; A1_vec = _mm256_load_pd(A); A += 4; __m256d C11_vec = A2_vec * B_vec; __m256d C12_vec = A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; k--; while (k--) { /* first iteration */ C1_vec += A1_vec * B_vec; A2_vec = _mm256_load_pd(A); A += 4; C2_vec += A2_vec * B_vec; A3_vec = _mm256_load_pd(A); A += 4; C3_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; C4_vec += A1_vec * B_vec; C5_vec += A2_vec * B_vec; C6_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; C7_vec += A1_vec * B_vec; C8_vec += A2_vec * B_vec; C9_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; C10_vec += A1_vec * B_vec; A1_vec = _mm256_load_pd(A); A += 4; C11_vec += A2_vec * B_vec; C12_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; /* second iteration */ C1_vec += A1_vec * B_vec; A2_vec = _mm256_load_pd(A); A += 4; C2_vec += A2_vec * B_vec; A3_vec = _mm256_load_pd(A); A += 4; C3_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; C4_vec += A1_vec * B_vec; C5_vec += A2_vec * B_vec; C6_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; C7_vec += A1_vec * B_vec; C8_vec += A2_vec * B_vec; C9_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; C10_vec += A1_vec * B_vec; A1_vec = _mm256_load_pd(A); A += 4; C11_vec += A2_vec * B_vec; C12_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; } C1_vec += A1_vec * B_vec; A2_vec = _mm256_load_pd(A); A += 4; C2_vec += A2_vec * B_vec; A3_vec = _mm256_load_pd(A); C3_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; C4_vec += A1_vec * B_vec; C5_vec += A2_vec * B_vec; C6_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); B++; C7_vec += A1_vec * B_vec; C8_vec += A2_vec * B_vec; C9_vec += A3_vec * B_vec; B_vec = _mm256_broadcast_sd(B); C10_vec += A1_vec * B_vec; C11_vec += A2_vec * B_vec; C12_vec += A3_vec * B_vec; /* [write-back] */ A1_vec = _mm256_load_pd(C); C1_vec += A1_vec; _mm256_store_pd(C, C1_vec); A2_vec = _mm256_load_pd(C + 4); C2_vec += A2_vec; _mm256_store_pd(C + 4, C2_vec); A1_vec = _mm256_load_pd(C + 8); C3_vec += A1_vec; _mm256_store_pd(C + 8, C3_vec); A2_vec = _mm256_load_pd(C + 12); C4_vec += A2_vec; _mm256_store_pd(C + 12, C4_vec); A1_vec = _mm256_load_pd(C + 16); C5_vec += A1_vec; _mm256_store_pd(C + 16, C5_vec); A2_vec = _mm256_load_pd(C + 20); C6_vec += A2_vec; _mm256_store_pd(C + 20, C6_vec); A1_vec = _mm256_load_pd(C + 24); C7_vec += A1_vec; _mm256_store_pd(C + 24, C7_vec); A2_vec = _mm256_load_pd(C + 28); C8_vec += A2_vec; _mm256_store_pd(C + 28, C8_vec); A1_vec = _mm256_load_pd(C + 32); C9_vec += A1_vec; _mm256_store_pd(C + 32, C9_vec); A1_vec = _mm256_load_pd(C + 36); C10_vec += A1_vec; _mm256_store_pd(C + 36, C10_vec); A2_vec = _mm256_load_pd(C + 40); C11_vec += A2_vec; _mm256_store_pd(C + 40, C11_vec); A1_vec = _mm256_load_pd(C + 44); C12_vec += A1_vec; _mm256_store_pd(C + 44, C12_vec); }  Well, as usual, lots of thanks to Peter Cordes! • For future readers: history of this code / question, as summarized by the OP in an answer to one of his early SO question about it. Then discussion on another Q led to this being posted in an incomplete state, but now it's all cleaned up. (thanks, Codereview.SE mods) – Peter Cordes Apr 4 '16 at 18:57 ## 1 Answer I put your code up on the Godbolt Compiler Explorer to check out the asm from the intrinsics version. You're right that even with gcc 5.3 or clang 3.8, there are spills / reloads in the inner loop. So you may actually get a speedup from hand-written asm here, if those store-forwarding round trips aren't hidden by out-of-order execution. The problem with using "ymm0" and so on in the clobber list (instead of "xmm0" is due to using a really old version of gcc. Declaring clobbers on ymm0 works fine even with gcc 4.9.2. (current is gcc 5.3, which you should definitely use. gcc5 made some significant improvements.) ### Comment your code The asm code is very lacking in comments. It's appropriate to describe what's going on in each block. It's also very appropriate for the C and asm versions to include a comment like what you have in your opening paragraph of the question that describes what they do, and what the 12 and 4 mean. ### Consider making the function a stand-alone asm function, instead of inline ASM Upside: you don't have to write a clobber list naming every vector register. Large amounts of asm code in GNU C inline asm looks messy. Downside: you lose portability between non-Windows (x86-64 SysV ABI) and Windows (__vectorcall). With inline asm, gcc will take care of the differences in calling convention, and save/restore the call-preserved xmm regs in the windows ABI. There's no performance advantage to be gained from it. gcc isn't wasting instructions. If you care about portability to Windows while still compiling with GNU C, it's probably easiest to keep using GNU C inline asm. If you need portability to MSVC, you might want to just compile this function stand-alone with gcc for Windows. ### Optimize for macro-fusion for Intel and AMD CPUs: "cmpl$1, %[k]\n\t"
"je  .L2\n\t"


should be next to each other so they can macro-fuse into a single compare-and-branch uop. See also Agner Fog's guides, including the microarchitecture pdf, and other stuff on the SO x86 tag wiki.

### Use local labels

instead of "K%=:\n\t", use ".LK%=:\n\t". Labels that start with .L don't appear in the object file, but other labels do. The disassembly output looks like (objdump -dwr gepdot_avx_12_4.o)

00000000000002b0 <GEPDOT_AVX__12_4_updatedasm>:
2b0:   48 8d 86 80 f8 ff ff    lea    -0x780(%rsi),%rax
...
327:   c5 05 59 fa             vmulpd %ymm2,%ymm15,%ymm15
32b:   83 ff 01                cmp    $0x1,%edi 32e: 0f 84 db 00 00 00 je 40f <End> 0000000000000334 <K177>: 334: c4 e2 e5 b8 e0 vfmadd231pd %ymm0,%ymm3,%ymm4 339: c5 fd 28 0e vmovapd (%rsi),%ymm1 33d: 48 81 c6 c0 00 00 00 add$0xc0,%rsi
...


because the auto-numbered (by %=) K177 label isn't a .L local label. That's weird, since normally you'd do that for a static function or something, not for labels within a global function.

It's hard to come up with good names for loop labels, but K is not a very good one. Even something like .Lmainloop%= would be ok. (To make it clear it's a loop, not an if()).

### One extra insn could save code size in the epilogue

Displacements outside [-128..127] take on instruction byte, but outside that range requires 4 instruction bytes.

 482:   c5 c5 58 79 60          vaddpd 0x60(%rcx),%ymm7,%ymm7
487:   c5 fd 29 79 60          vmovapd %ymm7,0x60(%rcx)
// insert an add here and change the displacements of all the following insns
48c:   c5 3d 58 81 80 00 00 00         vaddpd 0x80(%rcx),%ymm8,%ymm8
494:   c5 7d 29 81 80 00 00 00         vmovapd %ymm8,0x80(%rcx)
49c:   c5 35 58 89 a0 00 00 00         vaddpd 0xa0(%rcx),%ymm9,%ymm9
4a4:   c5 7d 29 89 a0 00 00 00         vmovapd %ymm9,0xa0(%rcx)
... 16 insns with 4B displacements.


You can put in an

"add     \$256, %[C]\n\t"   // spend 7 code bytes for this insn to save 16*3 bytes for displacements in following instructions


Actually 16*3 + 2, because one of the addresses can use no displacement at all. This increase in code density is probably worth it even considering the uop cache. With giant instructions, you're probably packing fewer than 6 fused-domain uops per cache line, since every 32B boundary has to start a new uop cache line. So the extra uop for the add instruction is probably balanced out by using fewer uop cache lines to store the rest of the uops.

It's outside the loop, so it's probably not uncommon for it to be executing from the decoders, rather than the uop-cache, in which case increasing code density is a big win.

### clang ends up with repeated vzeroupper insns

Pulling vzeroupper out of the inline asm and using _mm256_zeroupper() like I did in the version I linked on godbolt (above) didn't help. Even when the compiler can "see" that there's already a vzeroupper, it still adds an extra one. It's not a big problem, only 4 uops per function call (and they don't use any execution units; just issue/retire bandwidth).

nvm this part, didn't realize you were reusing the loaded vector in future instructions. In this case it's probably not worth it to write several instructions that all load the same data. The 4 fused-domain uops per clock pipeline width is not going to be a problem when your throughput is limited by doing mostly FMAs (throughout of only 2 per clock on current Intel HW). Thus, throwing in separate load instructions is fine.

You could maybe save code size and fused-domain uops by folding some loads into FMAs. You could probably replace:

    "vmovapd  32(%[A]), %%ymm2\n\t"


with a vfmadd with a memory operand. The 132, 213, and 231 versions give you some choice of which is the output and which is the memory source operand.

Like I said, this is only a good idea if you're not reusing the loaded data much / at all.

• Comments are not for extended discussion; this conversation has been moved to chat. – Mathieu Guindon Apr 5 '16 at 15:34