Computers Counting Part 6: SIMD with Vector Instructions
How fast can computers count?
This is the fifth post in a series of writing programs to see how fast computers can count. The series will start with naive methods, then will make use of optimization and later parallelism within the CPU, then will utilize a GPU, and the finale will use a dedicated FPGA for counting.
Previous posts:
- Naive counting
- Digging into the assembly
- Compiler Optimizations
- Handwritten Assembly
- Intro to Parallelism
Hardware specs:
- CPU: AMD Ryzen 9 3900X (12x Zen 2 cores)
- GPU: Nvidia GeForce RTX 5060
- RAM: 32GB DDR4
- FPGA: PYNQ-Z2
Software specs:
- OS: Fedora 43
- Compiler: gcc 15.2.1
Parallelism within a CPU core
All modern CPU instruction set architectures include a vector processing unit. This unit can perform operations on a vector of numbers rather than a scalar number. This is known as SIMD, or single-instruction, multiple-data. Instead of telling the CPU to add a scalar by 1, I can tell the CPU to add 1 to all scalars within a vector.
The concept of SIMD is crucial throughout this series, as the bulk of the instructions that the CPU (or eventual GPU) will be running will be the “add 1” instructions. It’s the same instruction regardless of where the addition is being done in the hardware.
Naive attempt
My computer is equipped with a AMD Ryzen 9 3900X, which contains 12 Zen 2 CPU cores. In this part, I will only use 1 core and will not utilize hyperthreading, as I am taking my first steps towards parallelism in this post. A Zen 2 core has the hardware to execute AVX-2 instructions, which is a modern (but not the latest and greatest that is AVX-512) set of vector instructions that operate on 256-bit registers. AVX-2 is an expansion of the old SSE instruction set.
Since the registers I am operating on are 256-bit, I can fit 8 32-bit ints within a register. Knowing this, I will create the first iteration of the code, with irrelevant parts omitted:
void setToZero(int* vec)
{
for(int i=0; i<8; i++)
{
vec[i] = 0;
}
}
void countByOne(int* vec)
{
for(int i=0; i<8; i++)
{
vec[i]++;
}
}
int main()
{
int countVector[8];
setToZero(countVector);
struct timeval startTime, endTime;
gettimeofday(&startTime, NULL);
while(countVector[0] < 100'000'000)
{
countByOne(countVector);
}
gettimeofday(&endTime, NULL);I initalize an array of 8 ints, then set all values of the array to zero in a crude fashion (this part doesn’t matter anyway), then have each element of the vector increment 100,000,000 times. Let’s see how this does:
Counted to 100 million in 3474890 microseconds.
Yikes, that’s awful. It’s a counting rate of only about 28.8 million per second, much less than all but the “very naive JavaScript” example. Why is it so bad?
Did you notice the bug?
There is a major bug in the code that prevents the paraellism from actually doing anything. I check the first element of the vector against the max of 100,000,000, but don’t sum anything up in the end. Let’s fix that.
gettimeofday(&startTime, NULL);
while(countVector[0] < 12'500'000)
{
countByOne(countVector);
}
int count = 0;
for(int i=0; i<8; i++)
{
count += countVector[i];
}
gettimeofday(&endTime, NULL);Here, I have each vector element count to 12,500,000, and then sum up each element of the vector to 100,000,000.
Counted to 100 million in 393495 microseconds.
That’s much better, but still not good at all, at 254 million increments per second. In fact, it’s almost exactly 8x faster than the last attempt. Is this parallelism thing all a waste?
Compiler optimizations
So far, I have compiled these programs with -O0, which optimizes for compile time, not runtime performance. gcc will not auto-vectorize this code on -O0, so it can’t be faster than running the naive C program. There is also a lot of overhead in calling the countByOne function, which takes in a pointer to an int, then does its own separate incrementing and pointer arithmetic (the vec[i] part), 100 million times.
Let’s see if compiling with optimizations on does anything.
- O1:
Counted to 100 million in 46455 microseconds.(2.15B increments/s) - O2:
Counted to 100 million in 3308 microseconds.(30.2B increments/s) - O3:
Counted to 100 million in 0 microseconds.(infinity increments/s)
Whoa! That’s a wide range of results. Before I start looking at the assembly, I will guess what is going on here.
- O1: The compiler has probably optimized the overhead of the countByOne function, but has probably not parallelized anything yet.
- O2: The compiler has either successfully parallelized the program with AVX-2 instructions, or has possibly cheated with some kind of a shortcut.
- O3: The compiler has definitely cheated by bypassing the counting altogether like it did in part 3 of this series.
While the entire code changes when setting compiler optimization levels, I will only focus on the relevant sections of assembly to see how the counting is done and to see if the counting is legitimate.
O0 Analysis
I purposely put all of the actual counting code in a function for readability of the assembly code. Let’s see what it is doing:
countByOne:
.LFB1:
.cfi_startproc
pushq %rbp
.cfi_def_cfa_offset 16
.cfi_offset 6, -16
movq %rsp, %rbp
.cfi_def_cfa_register 6
movq %rdi, -24(%rbp)
movl $0, -4(%rbp)
jmp .L5
.L6:
movl -4(%rbp), %eax
cltq
leaq 0(,%rax,4), %rdx
movq -24(%rbp), %rax
addq %rdx, %rax
movl (%rax), %edx
addl $1, %edx
movl %edx, (%rax)
addl $1, -4(%rbp)
.L5:
cmpl $7, -4(%rbp)
jle .L6
nop
nop
popq %rbp
.cfi_def_cfa 7, 8
retI’m not gonna analyze this line by line but it is clearly inefficient. There are 30 lines of assembly code with a lot of move instructions and even two “no-operation” instructions. It’s also performing its operations on the normal x86-64 registers rather than the vector registers that should be used, so it is not doing anything in parallel.
call gettimeofday
jmp .L8
.L9:
leaq -48(%rbp), %rax
movq %rax, %rdi
call countByOne
.L8:
movl -48(%rbp), %eax
cmpl $12499999, %eax
jle .L9
movl $0, -4(%rbp)
movl $0, -8(%rbp)
jmp .L10
.L11:
movl -8(%rbp), %eax
cltq
movl -48(%rbp,%rax,4), %eax
addl %eax, -4(%rbp)
addl $1, -8(%rbp)
.L10:
cmpl $7, -8(%rbp)
jle .L11
leaq -80(%rbp), %rax
movl $0, %esi
movq %rax, %rdi
call gettimeofdayLooking at the section between the two gettimeofday calls, the CPU is legitimately counting, including the summation at the end. Not much more to say, as the CPU spends the vast majority of its time in the countByOne function.
O1 Analysis
countByOne:
.LFB12:
.cfi_startproc
leaq 32(%rdi), %rax
.p2align 4
.L5:
addl $1, (%rdi)
addq $4, %rdi
cmpq %rax, %rdi
jne .L5
retHere, the assembly for countByOne is much simpler, so I will analyze it in further detail.
The function starts off with leaq setting register %rax to memory address of the 1st element of the count vector, plus 32. Then, it adds 1 to the current element of the count vector with addl. Next, it adds 4 to the pointer to the memory address of the current element of the count vector so that it points to the next element of the count vector. A comparison then takes place with register %rax, and if they are not equal, the execution jumps to .L5. The function returns when all 8 increments happen.
As you might have noticed, this is all done serially on regular x86-64 registers. This is much more efficient and allows for decent serial execution speed (2.15B increments/s), but does not utilize any kind of parallelism.
call gettimeofday
cmpl $12499999, 32(%rsp)
jg .L8
leaq 32(%rsp), %rbx
.L9:
movq %rbx, %rdi
call countByOne
cmpl $12499999, 32(%rsp)
jle .L9
.L8:
leaq 32(%rsp), %rax
leaq 64(%rsp), %rcx
movl $0, %edx
.p2align 4
.L10:
addl (%rax), %edx
movl %edx, %ebx
addq $4, %rax
cmpq %rcx, %rax
jne .L10
movq %rsp, %rdi
movl $0, %esi
call gettimeofdayLooking at the section between the two gettimeofday calls, the CPU is legitimately counting, including the summation at the end.
O2 Analysis
This is the most interesting and exciting result. Let’s first check if it is legitimate.
call gettimeofday
pxor %xmm0, %xmm0
pcmpeqd %xmm2, %xmm2
movdqa %xmm0, %xmm1
psrld $31, %xmm2
.p2align 5
.p2align 4
.p2align 3
.L5:
paddd %xmm2, %xmm1
paddd %xmm2, %xmm0
movd %xmm1, %eax
cmpl $12499999, %eax
jle .L5
leaq 48(%rsp), %rdi
xorl %esi, %esi
movaps %xmm1, 16(%rsp)
movaps %xmm0, (%rsp)
call gettimeofdayAt first glance, this looks legitimate. There’s a lot of stuff between the gettimeofday calls. However, a call to CountByOne is missing. Let’s analyze further.
Before I analyze further, I should explain what the %xmm registers are. These are 128-bit vector registers, meaning that the compiler is finally making use of the wider registers I alluded to in the beginning of this post. I expected to see the 256-bit %ymm registers here, but the compiler is likely using %xmm for compatability across all x86-64 CPUs. While %xmm registers are essentially universal across x86-64 CPUs, %ymm registers only were introduced on x86-64 CPUs in 2011, and weren’t universal on new CPUs until circa 2020.
The first several instructions before the .L5 tag set things up. The pxor is a zeroing idiom for %xmm0, the movdqa instruction copies the zero into %xmm1. The pcmpeqd and psrld instructions set %xmm2 to [1,1,1,1]. Note that the instructions beginning with p are SIMD instructions.
The first two paddd instructions add %xmm2 to %xmm1 and %xmm0. paddd is an opcode to do a vector addition between a source and destination vector of 32-bit ints. Since %xmm2 contains 1 in each of its elements, this effectively adds 1 to each element in %xmm0 and %xmm1. This means that there is finally an example of parallel counting in the code!
The movd instruction moves the low 32-bit integer element from %xmm1 to register %eax, which is legitimate because all elements in %xmm0 and %xmm1 are the same after all of the additions. The cmpl instruction compares the per-element count to 12,499,999, and jumps back to .L5 if 12,499,999 is less than the per-element count. The jumping will stop once the per-element count reaches 12,500,000.
Unfortunately, the next 4 instructions don’t do anything useful for summing the elements of %xmm0 and %xmm1 into an accumulator. There’s some useless pointer arithmetic, a useless zeroing idiom, and spilling of the values of %xmm1 and %xmm0 into the stack.
So no, this is not legitimate. It got most of the way there but didn’t sum everything up in the end. I can only say that this code counts to 12,500,000.
O3 Analysis
call gettimeofday
leaq 16(%rsp), %rdi
xorl %esi, %esi
call gettimeofdayThe CPU is clearly cheating here. End of story.
Fixing the O2 version
If I take a look at the lines of assembly immediately following the second gettimeofday call, I get this:
movdqa (%rsp), %xmm0
paddd 16(%rsp), %xmm0
xorl %eax, %eax
movq 48(%rsp), %rbx
subq 32(%rsp), %rbx
movl $.LC1, %edi
movdqa %xmm0, %xmm1
imulq $1000000, %rbx, %rbx
addq 56(%rsp), %rbx
subq 40(%rsp), %rbx
psrldq $8, %xmm1
paddd %xmm1, %xmm0
movdqa %xmm0, %xmm1
psrldq $4, %xmm1
paddd %xmm1, %xmm0
movd %xmm0, %esi
call printfFirst, notice that the instructions using the regular x86-64 registers are interspersed with the instructions using the vector registers. The former set of instructions are used to prepare the printf call and are unimportant in this analysis. The instructions are interspersed because this allows the CPU to more efficiently issue instructions to run the whole code faster overall. However, I’m interested in benchmarking the counting ability of the CPU, not how well it can mix that with function calls to other things.
The first two instructions reload %xmm0 from the stack, and then add %xmm1 to %xmm0. Since each element in %xmm0 and %xmm1 were 12,500,000, this doubles each element of %xmm0 to 25,000,000. The next movdqa instruction copies the new %xmm0 vector to %xmm1.
Then, the psrldq instruction shifts %xmm1 right by 8 bytes, effectively clearing the first two 32-bit ints in the %xmm1 vector to 0 but keeping the other two 32-bit ints the same. At this point, %xmm1 is equal to [0, 0, 25M, 25M] Another vector addition is done, leaving %xmm0 with a value of [25M, 25M, 50M, 50M]. This vector is copied onto %xmm1 again and is this time shifted right by 4 bytes, leaving %xmm1 with a value of [0, 25M, 25M, 50M]. This is added to %xmm0 to get [25M, 50M, 75M, 100M]. Finally, the last value of 100M is moved to register %esi. (All vectors are displayed as “little-endian” in this paragraph.)
So, the code does eventually get around to summing the vectors up, but it does so after the call to gettimeofday. Let’s fix it by moving the gettimeofday call to after the summation happens.
.L5:
paddd %xmm2, %xmm1
paddd %xmm2, %xmm0
movd %xmm1, %eax
cmpl $12499999, %eax
jle .L5
leaq 48(%rsp), %rdi
xorl %esi, %esi
movaps %xmm1, 16(%rsp)
movaps %xmm0, (%rsp)
movdqa (%rsp), %xmm0
paddd 16(%rsp), %xmm0
movdqa %xmm0, %xmm1
psrldq $8, %xmm1
paddd %xmm1, %xmm0
movdqa %xmm0, %xmm1
psrldq $4, %xmm1
paddd %xmm1, %xmm0
call gettimeofday
movd %xmm0, %esiCounted to 100 million in 3343 microseconds.
All I did here was swap a bunch of lines here so that the counting happens before the call to gettimeofday, while the other stuff happens after the gettimeofday call “offscreen”.
Note that the final move to the %esi register does happen after the gettimeofday call because the execution ran into a segfault and apparently %esi needs to be 0 in such a call according to Claude. Moving the call up one line is still legitimate in my books because the final summation to 100 million happens on the final paddd instruction, which is still before the call to gettimeofday.
Conclusion
Wow, that was a lot! Properly utilizing SIMD instructions take more work than properly utilizing regular instructions due to the added complexity, but it was worth it. The fixed -O2 version uses SIMD to count at a rate of almost 30 billion increments per second. This is about 8 times as fast as the best previous attempt, and also about 8 increments per clock cycle of the CPU!
Can we do better? Yes.
- Notice how the 128-bit %xmm registers were being used instead of the newer 256-bit %ymm registers. This could theoretially double throughput, although the CPU might not be able to dual-issue the adds like it did with the %xmm registers.
- Using 32-bit integers is actually not very efficient. The x86-64 SIMD instructions allow operations on not just 32-bit integers; I can go down to 8-bit integers to count up to 4 times as fast. However, I will need to periodically accumulate before each 8-bit integer rolls over, creating some overhead.
I will focus on these optimizations in the next part. Stay tuned!
Summary table of counting speed
| Method | Increments / second |
|---|---|
| Very naive JavaScript | 12.6M |
| Naive JavaScript | 369M |
| Naive C | 1.55B |
| Optimized C (-O0) | 2.02B |
| Optimized C (-O3) | 3.00B |
| Hand-optimized assembly | 3.77B |
| Basic SIMD (1 CPU core) | 29.9B |