SSE: Vector Floats for x86

• Single, meaning just one.
  
• Instruction, as in a machine code instruction, executed by hardware.
  
• Multiple, as in more than one--from 2 to a thousand or so.
  
• Data, as in floats or ints.
  

• x86, where they are called "SSE", or "Streaming SIMD Extensions"
• PowerPC, where they are called "AltiVec"
• Cell Broadband Engine, where they are part of the "Synergistic Processing Elements"
• Graphics cards, where they are part of the pixel processors

SSE Assembly

	movups xmm1,[thing1]; <- copy the four floats into xmm1
	addps xmm1,xmm1; <- add floats to themselves
	movups [retval],xmm1; <- move that constant into the global "retval"
	
; Print out retval
	extern farray_print
	push 4 ;<- number of floats to print
	push retval ;<- points to array of floats
	call farray_print
	add esp,8 ; <- pop off arguments
	ret

section .data
thing1: dd 10.2, 100.2, 1000.2, 10000.2;<- source constant
retval: dd 0.0, 0.0, 0.0, 0.0 ;<- our return value

(Try this in NetRun now!)

Here's a version that loads up two separate float constants, and adds them:

	movups xmm1,[thing1]; <- copy the four floats into xmm1
	movups xmm6,[thing2]; <- copy the four floats into xmm1
	addps xmm1,xmm6; <- add floats
	movups [retval],xmm1; <- move that constant into the global "retval"
	
; Print out retval
	extern farray_print
	push 4 ;<- number of floats to print
	push retval ;<- points to array of floats
	call farray_print
	add esp,8 ; <- pop off arguments
	ret

section .data
thing1: dd 10.2, 100.2, 1000.2, 10000.2;<- source constant
thing2: dd 1.2, 2.2, 3.2, 4.2;<- source constant
retval: dd 0.0, 0.0, 0.0, 0.0 ;<- our return value

(Try this in NetRun now!)

movss xmm3,[pi]; load up constant
addss xmm3,xmm3 ; add pi to itself
cvtss2si eax,xmm3 ; round to integer
ret
section .data
pi: dd 3.14159265358979 ; constant

(Try this in NetRun now!)

movss xmm3,[pi]; load up constant
addss xmm3,xmm3 ; add pi to itself
movss [output],xmm3; write register out to memory

; Print floating-point output
mov rdi,output ; first parameter: pointer to floats
mov rsi,1 ; second parameter: number of floats
sub rsp,8 ; keep stack 16-byte aligned (else get crash!)
extern farray_print
call farray_print
add rsp,8

ret

section .data
pi: dd 3.14159265358979 ; constant
output: dd 0.0 ; overwritten at runtime

(Try this in NetRun now!)

SSE in C/C++

• GCC 3.x under linux or MacOS X.  Be sure to compile with "-msse", or the header will whine.
  
• Visual C++ 2003 or better under Windows.
  
• Intel C/C++ compiler under Linux or Windows. Makes bogus complaint about "no EMMS instruction before return", although the annoying emms() call is thankfully only needed with MMX, not SSE.
  

	for (int i=0;i<n_vals;i++) { 
	        vals[i]=vals[i]*a+b;
	}

	for (int i=0;i<n_vals;i+=4) { 
	        vals[i+0]=vals[i+0]*a+b;
	        vals[i+1]=vals[i+1]*a+b;
	        vals[i+2]=vals[i+2]*a+b;
	        vals[i+3]=vals[i+3]*a+b;
	}

	__m128 SSEa=_mm_load1_ps(&a);
	__m128 SSEb=_mm_load1_ps(&b);
	__m128 v=_mm_load_ps(&vec[i]);
	v=_mm_add_ps(_mm_mul_ps(v,SSEa),SSEb);
	_mm_store_ps(&vec[i],v);

(Try this in NetRun now!)

• Advantage: it's easy to branch, or call a function, on just one 3-vector.
• Disadvantage: data may need to be padded to occupy all four floats.
• Disadvantage: non-parallel operations like dot product may need some fiddly code.
  

• Advantage: better packing of data into four-float vectors.
• Advantage: non-parallel operations like dot product are easy to code (SSEx+SSEy+SSEz)
  
• Advantage: clearly can scale to longer SSE vector lengths, like with AVX.
  
• Disadvantage: branching per-float is much trickier (see below).

AoS vs SoA: SSE Matrix-Vector Multiply

enum {n=4};
float mat[n][n];
float vec[n];
float outvector[n];

int foo(void) {
 for (int row=0;row<4;row++) {
	 float sum=0.0;
	 for (int col=0;col<n;col++) {
		   float m=mat[row][col];
		   float v=vec[col];
		   sum+=m*v;
	 }
	 outvector[row]=sum;
 }
 return 0;
}

(Try this in NetRun now!)

enum {n=4};
float mat[n][n];
float vec[n];
float outvector[n];

int foo(void) {
	for (int row=0;row<4;row++) {
		float sum=0.0, m,v;
		m=mat[row][0];
		v=vec[0];
		sum+=m*v;
		m=mat[row][1];
		v=vec[1];
		sum+=m*v;
		m=mat[row][2];
		v=vec[2];
		sum+=m*v;
		m=mat[row][3];
		v=vec[3];
		sum+=m*v;
		outvector[row]=sum;
	}
	return 0;
}

(Try this in NetRun now!)

#include <pmmintrin.h>

enum {n=4};
__m128 mat[n]; /* rows */
__m128 vec;
float outvector[n];

int foo(void) {
 for (int row=0;row<n;row++) {
   __m128 mrow=mat[row];
   __m128 v=vec;
   __m128 sum=mrow*v;
   sum=_mm_hadd_ps(sum,sum); /* adds adjacent-two floats */
   _mm_store_ss(&outvector[row],_mm_hadd_ps(sum,sum)); /* adds those floats */
 }
 return 0;
}

(Try this in NetRun now!)

#include <xmmintrin.h>

enum {n=4};
__m128 mat[n]; /* by column */
float vec[n];
__m128 outvector;

int foo(void) {
 float z=0.0;
 __m128 sum=_mm_load1_ps(&z);
 for (int col=0;col<n;col++) {
   __m128 mcol=mat[col];
   float v=vec[col];
   sum+=mcol*_mm_load1_ps(&v);
 }
 outvector=sum;
 return 0;
}

(Try this in NetRun now!)

The bottom line is that, as is *typical*, just making the inner loop parallel requires extensive changes to the data layout and processing throughout the program.

Per-Float Branching in SSE

• Comparison instructions, like cmpeq, that fill the corresponding SSE register's float with 0xffffffff (binary all ones) if the corresponding float component comparison comes out true, and all zeros if the comparison comes out false.  That is, these are per-component comparisons, so you could get an output like {0xffffffff,0x0,0xffffffff,0x0} if the even two components come out true, and the odd two components come out false.
  
• Logical operations, like and, nand, or, and xor, that bitwise-and SSE registers.  Note that bitwising AND'ing two floats is totally useless (the sign, exponent, and mantissa fields get anded!).  But bitwise AND'ing with all-zeros zeros out the corresponding float, and AND'ing with all ones leaves the float untouched.
  

	for (int i=0;i<n;i++) { 
	        if (vec[i]<7) 
			vec[i]=vec[i]*a+b;
		else
			vec[i]=c;
	}

	for (int i=0;i<n;i++) { 
	        unsigned int mask=(vec[i]<7)?0xffFFffFF:0;
		vec[i]=((vec[i]*a+b)&mask) | (c&~mask);
	}

	__m128 A=_mm_load1_ps(&a), B=_mm_load1_ps(&b), C=_mm_load1_ps(&c);
	__m128 Thresh=_mm_load1_ps(&thresh);
	for (int i=0;i<n;i+=4) { 
		__m128 V=_mm_load_ps(&vec[i]);
		__m128 mask=_mm_cmplt_ps(V,Thresh);   // Do all four comparisons
		__m128 V_then=_mm_add_ps(_mm_mul_ps(V,A),B); // "then" half of "if"
		__m128 V_else=C;                             // "else" half of "if"
		V=_mm_or_ps( _mm_and_ps(mask,V_then), _mm_andnot_ps(mask,V_else) );
		_mm_store_ps(&vec[i],V);
	}

(Try this in NetRun now!)

Hiding SSE Nastiness In a Wrapper Class

#include <xmmintrin.h>

class fourfloats; /* forward declaration */

/* Wrapper around four bitmasks: 0 if false, all-ones (NAN) if true. 
  This class is used to implement comparisons on SSE registers.
*/
class fourmasks {
	__m128 mask;
public:
	fourmasks(__m128 m) {mask=m;}
	__m128 if_then_else(fourfloats dthen,fourfloats delse);
};

/* Nice wrapper around __m128:
   it represents four floating point values. */
class fourfloats {
	__m128 v;
public:
	fourfloats(float onevalue) { v=_mm_load1_ps(&onevalue); }
	
	fourfloats(__m128 ssevalue) {v=ssevalue;} // put in an SSE value
	operator __m128 () const {return v;} // take out an SSE value

	fourfloats(const float *fourvalues) { v=_mm_load_ps(fourvalues);}
	void store(float *fourvalues) {_mm_store_ps(fourvalues,v);}

/* arithmetic operations return blocks of floats */
	fourfloats operator+(const fourfloats &right) {
		return _mm_add_ps(v,right.v);
	}	

/* comparison operations return blocks of masks (bools) */
	fourmasks operator<(const fourfloats &right) {
		return _mm_cmplt_ps(v,right.v);
	}
};

inline __m128 fourmasks::if_then_else(fourfloats dthen,fourfloats delse) {
	return _mm_and_ps(mask,dthen)+
		_mm_andnot_ps(mask,delse);
}

float src[4]={1.0,5.0,3.0,4.0};
float dest[4];
int foo(void) {
/*
	// Serial code
	for (int i=0;i<4;i++) {
		if (src[i]<4.0) dest[i]=src[i]*2.0; else dest[i]=17.0;
	}
*/
	// Parallel code
	fourfloats s(src);
	fourfloats d=(s<4.0).if_then_else(s+s,17.0);
	d.store(dest);

	//farray_print(dest,4); // <- for debugging
	return 0;
}

(Try this in NetRun now!)

Compilers are now good enough that there is zero speed penalty due to the nice "fourfloats" class: the nice syntax comes for free!