Noah Kochavi's blog

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:

  1. Naive counting
  2. Digging into the assembly
  3. Compiler Optimizations
  4. Handwritten Assembly
  5. Intro to Parallelism

Hardware specs:

Software specs:

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.

Whoa! That’s a wide range of results. Before I start looking at the assembly, I will guess what is going on here.

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
	ret

I’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	gettimeofday

Looking 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
	ret

Here, 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	gettimeofday

Looking 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	gettimeofday

At 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	gettimeofday

The 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	printf

First, 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, %esi

Counted 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.

  1. 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.
  2. 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