SSE and AVX: SIMD for x86

SSE in 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 or newer 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.

Avoiding Branching with the "If Then Else" Trick

if (x) y=a; else y=b;   // conditional

y=b+x*(a-b); // linear version (assumes a-b works)
y=x*a+(1-x)*b; // rearranged version of above
y=(x&a)|(~x)&b); // bitwise version (assumes x is all zero bits, or all one bits)

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!

#include <immintrin.h> /* Intel SIMD intrinsics header (use emmintrin.h if this fails) */

class floats; /* forward declaration */

/* This is a SIMD block of boolean flags. */
class booleans {
	__m128 v;
public:
	enum {n=4}; // how many bools we work on at once
	
	booleans(__m128 n) :v(n) {}

	/* Logical operations on booleans: */
	friend booleans operator&&(const booleans &a,const booleans &b) {return _mm_and_ps(a.v,b.v);}
	friend booleans operator||(const booleans &a,const booleans &b) {return _mm_or_ps(a.v,b.v);}

	/*
	Return "the_then" where we are true, and "the_else" where we are false.
	Basically performs a per-element branch, but without branching.
	Example:
		booleans m=(a<b);
		c=m.then_else(3,7);
	*/
	floats then_else(const floats &the_then,const floats &the_else);
};

/* This is a SIMD block of single-precision floats. */
class floats {
	__m128 v;
public:
	enum {n=4}; // floats::n is how many floats we work on at once

	floats() {}
	floats(__m128 n) :v(n) {}
	floats(const float *src) {v=_mm_loadu_ps(src);}
	floats(float single) {v=_mm_load1_ps(&single);}

	operator __m128 (void) const {return v;}
	
	void store(float *dest) const {_mm_storeu_ps(dest,v);}
	/* Caution: this only works if dest is 16-byte aligned. */
	void store_aligned(float *dest) const {_mm_store_ps(dest,v);}

	/* Printing support */
	friend std::ostream &operator<<(std::ostream &o,const floats &a) {
		float out[floats::n];
		a.store(out);
		for (int i=0;i<floats::n;i++) std::cout<<out[i]<<"  ";
		return o;
	}

	/* Arithmetic operations */
	friend floats operator+(const floats &a,const floats &b) {return _mm_add_ps(a.v,b.v);}
	friend floats operator-(const floats &a,const floats &b) {return _mm_sub_ps(a.v,b.v);}
	friend floats operator*(const floats &a,const floats &b) {return _mm_mul_ps(a.v,b.v);}
	friend floats operator/(const floats &a,const floats &b) {return _mm_div_ps(a.v,b.v);}
	friend floats min(const floats &a,const floats &b) {return _mm_min_ps(a.v,b.v);}
	friend floats max(const floats &a,const floats &b) {return _mm_max_ps(a.v,b.v);}
	friend floats sqrt(const floats &a) {return _mm_sqrt_ps(a.v);}

	/* Comparison operators */
	friend booleans operator<(const floats &a,const floats &b) {return _mm_cmplt_ps(a.v,b.v);}
	friend booleans operator<=(const floats &a,const floats &b) {return _mm_cmple_ps(a.v,b.v);}
	friend booleans operator>=(const floats &a,const floats &b) {return _mm_cmpge_ps(a.v,b.v);}
	friend booleans operator>(const floats &a,const floats &b) {return _mm_cmpgt_ps(a.v,b.v);}
	friend booleans operator==(const floats &a,const floats &b) {return _mm_cmpeq_ps(a.v,b.v);}
	friend booleans operator!=(const floats &a,const floats &b) {return _mm_cmpneq_ps(a.v,b.v);}
	
};

inline floats booleans::then_else(const floats &the_then,const floats &the_else)
{
	return _mm_or_ps(
		_mm_and_ps(v,the_then),
		_mm_andnot_ps(v,the_else)
	);
}

/* A little test code: */
	float a[4]={1.375,1.567,1.876,1.999};
	float b[4]={10.0,0.1,0.2,10.3};
	float c[4]={5.0,5.1,5.2,5.3};
	float d[4]={8.1,8.2,8.3,8.4};
	float out[4]; 

int foo(void) {
	floats A(a), B(b);
	booleans m=(A<B);
	std::cout<<"then_else="<<m.then_else(3,c)<<"\n";
	return 0;
}

(Try this in NetRun now!)

Intel's Advanced Vector Extensions (AVX)

#include <immintrin.h> /* AVX + SSE4 intrinsics header */

// Internal class: do not use...
class not_vec8 {
	__m256 v; // bitwise inverse of our value (!!)
public:
	not_vec8(__m256 val) {v=val;}
	__m256 get(void) const {return v;} // returns INVERSE of our value (!!)
};

// This is the class to use!
class vec8 {
	__m256 v;
public:
	vec8(__m256 val) {v=val;}
	vec8(const float *src) {v=_mm256_loadu_ps(src);}
	vec8(float x) {v=_mm256_broadcast_ss(&x);}
	
	vec8 operator+(const vec8 &rhs) const {return _mm256_add_ps(v,rhs.v);}
	vec8 operator-(const vec8 &rhs) const {return _mm256_sub_ps(v,rhs.v);}
	vec8 operator*(const vec8 &rhs) const {return _mm256_mul_ps(v,rhs.v);}
	vec8 operator/(const vec8 &rhs) const {return _mm256_div_ps(v,rhs.v);}
	vec8 operator&(const vec8 &rhs) const {return _mm256_and_ps(v,rhs.v);}
	vec8 operator|(const vec8 &rhs) const {return _mm256_or_ps(v,rhs.v);}
	vec8 operator^(const vec8 &rhs) const {return _mm256_xor_ps(v,rhs.v);}
	vec8 operator==(const vec8 &rhs) const {return _mm256_cmp_ps(v,rhs.v,_CMP_EQ_OQ);}
	vec8 operator!=(const vec8 &rhs) const {return _mm256_cmp_ps(v,rhs.v,_CMP_NEQ_OQ);}
	vec8 operator<(const vec8 &rhs) const {return _mm256_cmp_ps(v,rhs.v,_CMP_LT_OQ);}
	vec8 operator<=(const vec8 &rhs) const {return _mm256_cmp_ps(v,rhs.v,_CMP_LE_OQ);}
	vec8 operator>(const vec8 &rhs) const {return _mm256_cmp_ps(v,rhs.v,_CMP_GT_OQ);}
	vec8 operator>=(const vec8 &rhs) const {return _mm256_cmp_ps(v,rhs.v,_CMP_GT_OQ);}

	not_vec8 operator~(void) const {return not_vec8(v);}

	__m256 get(void) const {return v;}

	float *store(float *ptr) {
		_mm256_store_ps(ptr,v);
		return ptr;
	}

	float &operator[](int index) { return ((float *)&v)[index]; }
	float operator[](int index) const { return ((const float *)&v)[index]; }

	friend ostream &operator<<(ostream &o,const vec8 &y) {
		o<<y[0]<<" "<<y[1]<<" "<<y[2]<<" "<<y[3];
		return o;
	}
	friend vec8 operator&(const vec8 &lhs,const not_vec8 &rhs)
		{return _mm256_andnot_ps(rhs.get(),lhs.get());}
	friend vec8 operator&(const not_vec8 &lhs, const vec8 &rhs)
		{return _mm256_andnot_ps(lhs.get(),rhs.get());}
	vec8 if_then_else(const vec8 &then,const vec8 &else_part) const {
		return _mm256_or_ps( _mm256_and_ps(v,then.v),
			 _mm256_andnot_ps(v, else_part.v)
		);
	}
};

vec8 sqrt(const vec8 &v) {
	return _mm256_sqrt_ps(v.get());
}
vec8 rsqrt(const vec8 &v) {
	return _mm256_rsqrt_ps(v.get());
}
/* Return value = dot product of a & b, replicated 4 times */
inline vec8 dot(const vec8 &a,const vec8 &b) {
	vec8 t=a*b;
	__m256 vt=_mm256_hadd_ps(t.get(),t.get());
	vt=_mm256_hadd_ps(vt,vt);
	return vt;
}

float data[8]={1.2,2.3,3.4,1.5,10.0,100.0,1000.0,10000.0};
vec8 a(data);
vec8 b(10.0);
vec8 c(0.0);

int foo(void) {
	c=(a<b).if_then_else(3.7,c);
	return 0;
}

(Try this in NetRun now!)