🔒 Locked

Vectorising a Dot Product

Started by DarkMatter2008 Apr 10, 2009 at 12:41 PM 7 replies 5.2k views

Original Post

DarkMatter2008

Apr 10, 2009 at 12:41 PM

Hello all, I'm trying to vectorise a dot product calculation using SSE intrinsics, I am new to this, what I'm trying to get is 4 dot products:- a . e b . e c . e d . e This is what I have at the moment.


__declspec(align(16)) float x1[4] = {a.x, b.x, c.x, d.x};
__declspec(align(16)) float y1[4] = {a.y, b.y, c.y, d.y};
__declspec(align(16)) float z1[4] = {a.z, b.z, c.z, d.z};

__declspec(align(16)) float vectresult[4];

__m128 m1, m2, m3, m4;

__m128* pSrc1 = (__m128*) x1;
__m128* pSrc2 = (__m128*) y1;
__m128* pSrc3 = (__m128*) z1;
__m128* pOut = (__m128*) vectresult;

while(some condition)
{
updateEvector();

__m128 m0_x = _mm_set_ps1(e.x);
__m128 m0_y = _mm_set_ps1(e.y); 
__m128 m0_z = _mm_set_ps1(e.z);  

m1 = _mm_mul_ps(*pSrc1, m0_x); 
m2 = _mm_mul_ps(*pSrc2, m0_y); 
m3 = _mm_mul_ps(*pSrc3, m0_z);
m4 = _mm_add_ps(m1, m2);
m4 = _mm_add_ps(m4, m3);

*pOut = m4;
}

It is currently getting beaten with the tradional method using the FPU:-


while(some condition)
{
updateEvector();

ae = (a.x*e.x) + (a.y*e.y) + (a.z*e.z);
be = (b.x*e.x) + (b.y*e.y) + (b.z*e.z);
ce = (c.x*e.x) + (c.y*e.y) + (c.z*e.z);
de = (d.x*e.x) + (d.y*e.y) + (d.z*e.z);
}

I did see an SSE implementation for a dot product with asm on this forum, but that was still slower then using the FPU. Is it possible to do what I want to do and beat the FPU with SSE for dot products or am I wasting my time? Thnaks in advance

clb

Apr 10, 2009 at 1:58 PM

What you're doing sounds rational, but (as you've noted by now) for any application of vectorization you need to check the two issues:

1) The performance cost incurred by setup and store. These are things like alignment restrictions (cache-affecting), having to shuffle data between normal float stack and SSE/MMX registers (vectorized accesses interleaved with scalar accesses etc.)
2) The performance gained in vectorizing the actual arithmetic. In an ideal case, like the one you described, you get to fold exactly 4 FPU ops to a single SSE op on a full 4-way data, throughout all the instructions you were doing. If you had branching where you couldn't even apply SSE op masks, the performance is a lot less.

On a basic level that's it - then it comes down to whether 2) outweighs 1).

In your case you have the ideal setup with regards to 2) - You're just doing the same operation (dot product) magnified 4x in a SOA fashion, can't even hope for anything better. If your whole 'meat' of the computation is just in doing dot products on an array of data, then you should be expecting a near 4x increase there. But I suppose that is not the case? If your dot-array-of-vectors-with-e -routine is just a small part of a bigger algorithm, then the overall performance greatly depends on how you are approaching 1). I.e. perhaps the SSE setup doesn't work smoothly with the rest of the processing you're doing before and after that step. For example, it might be that your vector was in FPU stack, needs to be moved into SSE registers just for the dot product, then needs to be put back in FPU stack for further processing.

Now, to approach this dot-product on a local scale (without knowing about the outside code), I'll guess what happens in your code is that __m128* pSrc1 = (__m128*) x1; doesn't result in optimal code. Try using _mm_load_ps and replacing those __m128* with just __m128. But that's just a first guess. Can you post the generated .asm output? (for both FPU and SSE versions)

Also, can you post your profiling results (ticks/ms), how far off are these two versions?

DarkMatter2008

Apr 10, 2009 at 5:54 PM

Thanks for your reply

Here are the times as requested

SSE processing time: 66 ms
Normal processing time: 32 ms

Here's what I have done:-

#include "stdafx.h"#include <stdlib.h>#include <iostream>#include "Windows.h"#include <xmmintrin.h>LARGE_INTEGER m_nFreq;LARGE_INTEGER m_nBeginTime;typedef struct VECTOR{	float x;	float y;	float z;}Vector;void tic(){	QueryPerformanceFrequency(&m_nFreq);	QueryPerformanceCounter(&m_nBeginTime);}__int64 toc(){	LARGE_INTEGER nEndTime;	__int64 nCalcTime;	QueryPerformanceCounter(&nEndTime);	nCalcTime = (nEndTime.QuadPart - m_nBeginTime.QuadPart) * 1000 / m_nFreq.QuadPart;	wprintf(L" processing time: %d ms\n", nCalcTime);	return nCalcTime;}Vector v[1000000];int _tmain(int argc, _TCHAR* argv[]){	for(int i =0; i<1000000; i++){		v.x =  float(rand() % 64 - 32);		v.y =  float(rand() % 64 - 32);		v.z =  float(rand() % 64 - 32);	}	int count = 0;	Vector a, b, c, d;	a.x = float(rand() % 64 - 32);	a.y = float(rand() % 64 - 32);	a.z = float(rand() % 64 - 32);	b.x = float(rand() % 64 - 32);	b.y = float(rand() % 64 - 32);	b.z = float(rand() % 64 - 32);	c.x = float(rand() % 64 - 32);	c.y = float(rand() % 64 - 32);	c.z = float(rand() % 64 - 32);	d.x = float(rand() % 64 - 32);	d.y = float(rand() % 64 - 32);	d.z = float(rand() % 64 - 32);	__declspec(align(16)) float x1[4] = {a.x, b.x, c.x, d.x};	__declspec(align(16)) float y1[4] = {a.y, b.y, c.y, d.y};	__declspec(align(16)) float z1[4] = {a.z, b.z, c.z, d.z};	__declspec(align(16)) float vectresult[4];	__m128 m1, m2, m3, m4, m0_x, m0_y, m0_z;	__m128 pSrc1 =  _mm_load_ps(&x1[0]);	__m128 pSrc2 =  _mm_load_ps(&y1[0]);	__m128 pSrc3 =  _mm_load_ps(&z1[0]);	__m128* pOut = (__m128*) vectresult;	float ae, be, ce, de;	Vector e;	printf("SSE\n");	tic();	while(count < 1000000){		e = v[count];		m0_x = _mm_set_ps1(e.x);		m0_y = _mm_set_ps1(e.y); 		m0_z = _mm_set_ps1(e.z);  		m1 = _mm_mul_ps(pSrc1, m0_x); 		m2 = _mm_mul_ps(pSrc2, m0_y); 		m3 = _mm_mul_ps(pSrc3, m0_z);		m4 = _mm_add_ps(m1, m2);		m4 = _mm_add_ps(m4, m3);		*pOut = m4;		count++;	}	toc();	count = 0;	printf("Normal\n");	tic();	while(count < 1000000){		e = v[count];		ae = (a.x*e.x) + (a.y*e.y) + (a.z*e.z);		be = (b.x*e.x) + (b.y*e.y) + (b.z*e.z);		ce = (c.x*e.x) + (c.y*e.y) + (c.z*e.z);		de = (d.x*e.x) + (d.y*e.y) + (d.z*e.z);		count++;	}	toc();	std::cin.get(); //wait for keypress	return 0;}

Ravyne

Apr 10, 2009 at 6:33 PM

If you haven't already, you have to read This: http://www.cortstratton.org/articles/OptimizingForSSE.php

Its a case study in optimizing Matrix-Vector multiplies, which is essentially the same as what you're doing. It should give you some insight into how to get further towards your goal, but it may be that the setup overhead doesn't overcome the smallish batch size of 4 -- but if you can batch more up at once, you should get massive performance increases. The author takes a naive C++ implementation (100 cycles) to an optimized SSE implementation over several steps, ultimately getting it down to around 17 cycles.

DarkMatter2008

Apr 11, 2009 at 8:02 AM

I had a look at what they did, I copied it over across to see it in action, but it under performs by an increadable ammount. After reports of success every where else on using SSE I feel that I missing something fundamental in my knowledge.

Here's what I now have.
C++

#include "stdafx.h"#include <stdlib.h>#include <iostream>#include "Windows.h"#include <xmmintrin.h>LARGE_INTEGER m_nFreq;LARGE_INTEGER m_nBeginTime;typedef float Matrix[4][3];typedef __declspec(align(16)) float Matrix4f[4][4];typedef struct VECTOR{	float x;	float y;	float z;	float w;}Vector;void tic(){	QueryPerformanceFrequency(&m_nFreq);	QueryPerformanceCounter(&m_nBeginTime);}__int64 toc(){	LARGE_INTEGER nEndTime;	__int64 nCalcTime;	QueryPerformanceCounter(&nEndTime);	nCalcTime = (nEndTime.QuadPart - m_nBeginTime.QuadPart) * 1000 / m_nFreq.QuadPart;	wprintf(L" processing time: %d ms\n", nCalcTime);	return nCalcTime;}Vector v[1000000];void MatrixMultiply3(Matrix4f &m, float *vin, float *vout){   // Get a pointer to the elements of m   float *row0 = reinterpret_cast<float*>(m);   __asm {      mov         esi, vin      mov         edi, vout      // load columns of matrix into xmm4-7      mov         edx, row0      movaps   xmm4, [edx]      movaps   xmm5, [edx+0x10]      movaps   xmm6, [edx+0x20]      movaps   xmm7, [edx+0x30]      // load v into xmm0.      movaps   xmm0, [esi]      // we'll store the final result in xmm2; initialize it      // to zero      xorps      xmm2, xmm2      // broadcast x into xmm1, multiply it by the first      // column of the matrix (xmm4), and add it to the total      movaps   xmm1, xmm0      shufps   xmm1, xmm1, 0x00      mulps      xmm1, xmm4      addps      xmm2, xmm1      // repeat the process for y, z and w      movaps   xmm1, xmm0      shufps   xmm1, xmm1, 0x55      mulps      xmm1, xmm5      addps      xmm2, xmm1      movaps   xmm1, xmm0      shufps   xmm1, xmm1, 0xAA      mulps      xmm1, xmm6      addps      xmm2, xmm1      movaps   xmm1, xmm0      shufps   xmm1, xmm1, 0xFF      mulps      xmm1, xmm7      addps      xmm2, xmm1      // write the results to vout      movups   [edi], xmm2   }}void BatchTransform1(Matrix4f &m, float *vin, float *vout, int len){   // initializations in C++ land  // Matrix4f mt(m, Matrix4f::TRANSPOSE); // work from a    float *row0 = reinterpret_cast<float*>(m);   static const int vecSize = 2 * sizeof(float[4]);      // if there are an odd number of vectors, process the first one   // separately and advance the pointers   if (len & 0x1) {      MatrixMultiply3(m, vin, vout);      ++vin;      ++vout;   }   len >>= 1; // we process two vectors at a time   __asm {      mov      esi, vin      mov      edi, vout      mov      ecx, len      // load columns of matrix into xmm4-7      mov      edx, row0      movaps   xmm4, [edx]      movaps   xmm5, [edx+0x10]      movaps   xmm6, [edx+0x20]      movaps   xmm7, [edx+0x30]BT2_START:      // process x (hiding the prefetches in the delays)      movss      xmm1, [esi+0x00]      movss      xmm3, [esi+0x10]      shufps   xmm1, xmm1, 0x00      prefetchnta [edi+0x30]      shufps   xmm3, xmm3, 0x00      mulps      xmm1, xmm4      prefetchnta   [esi+0x30]      mulps      xmm3, xmm4      // process y      movss      xmm0, [esi+0x04]      movss      xmm2, [esi+0x14]      shufps   xmm0, xmm0, 0x00      shufps   xmm2, xmm2, 0x00      mulps      xmm0, xmm5      mulps      xmm2, xmm5      addps      xmm1, xmm0      addps      xmm3, xmm2      // process z (hiding some pointer arithmetic between      // the multiplies)      movss      xmm0, [esi+0x08]      movss      xmm2, [esi+0x18]      shufps   xmm0, xmm0, 0x00      shufps   xmm2, xmm2, 0x00      mulps      xmm0, xmm6      add         esi, vecSize      mulps      xmm2, xmm6      add         edi, vecSize      addps      xmm1, xmm0      addps      xmm3, xmm2      // process w      addps      xmm1, xmm7      addps      xmm3, xmm7      // write output vectors to memory and loop      movaps   [edi-0x20], xmm1      movaps   [edi-0x10], xmm3      dec         ecx      jnz         BT2_START   }}inline void transformVector1(Matrix matrix, Vector *in, Vector *out, const float zAssumption){	out->x  = in->x * matrix[0][0] +		in->y * matrix[0][1] +		zAssumption * matrix[0][2];	out->y = in->x * matrix[1][0] +		in->y * matrix[1][1] +		zAssumption * matrix[1][2];	out->z = in->x * matrix[2][0] +		in->y * matrix[2][1] +		zAssumption * matrix[2][2];	out->w = in->x * matrix[3][0] +		in->y * matrix[3][1] +		zAssumption * matrix[3][2];}inline void transformVector2(Matrix matrix, Vector *in, const float zAssumption, const int len){	int i = 0;	float x;	float y;	do 	{		x = in->x;		y = in->y;		in->x  = x * matrix[0][0] +			y * matrix[0][1] +			zAssumption * matrix[0][2];		in->x = x * matrix[1][0] +			y * matrix[1][1] +			zAssumption * matrix[1][2];		in->x = x * matrix[2][0] +			y * matrix[2][1] +			zAssumption * matrix[2][2];		in->x = x * matrix[3][0] +			y * matrix[3][1] +			zAssumption * matrix[3][2];		in++;		i++;	}while(i < len);}__declspec(align(16)) float in[1000000][4];__declspec(align(16)) float out[1000000][4];int _tmain(int argc, _TCHAR* argv[]){	float zAssumption = float(rand() % 64 - 32);	for(int i =0; i<1000000; i++){		v.x =  float(rand() % 64 - 32);		v.y =  float(rand() % 64 - 32);		v.z =  zAssumption;	}	int count = 0;	Vector a, b, c, d;	a.x = float(rand() % 64 - 32);	a.y = float(rand() % 64 - 32);	a.z = float(rand() % 64 - 32);	b.x = float(rand() % 64 - 32);	b.y = float(rand() % 64 - 32);	b.z = float(rand() % 64 - 32);	c.x = float(rand() % 64 - 32);	c.y = float(rand() % 64 - 32);	c.z = float(rand() % 64 - 32);	d.x = float(rand() % 64 - 32);	d.y = float(rand() % 64 - 32);	d.z = float(rand() % 64 - 32);	Matrix mat;	mat[0][0] = a.x;	mat[0][1] = a.y;	mat[0][2] = a.z;	mat[1][0] = b.x;	mat[1][1] = b.y;	mat[1][2] = b.z;	mat[2][0] = c.x;	mat[2][1] = c.y;	mat[2][2] = c.z;	mat[3][0] = d.x;	mat[3][1] = d.y;	mat[3][2] = d.z;	Matrix4f matCol;	mat[0][0] = a.x;	mat[1][0] = a.y;	mat[2][0] = a.z;	mat[3][0] = 1;	mat[0][1] = b.x;	mat[1][1] = b.y;	mat[2][1] = b.z;	mat[3][1] = 1;	mat[0][2] = c.x;	mat[1][2] = c.y;	mat[2][2] = c.z;	mat[3][2] = 1;	mat[0][3] = d.x;	mat[1][3] = d.y;	mat[2][3] = d.z;	mat[3][3] = 1;	__declspec(align(16)) float x1[4] = {a.x, b.x, c.x, d.x};	__declspec(align(16)) float y1[4] = {a.y, b.y, c.y, d.y};	__declspec(align(16)) float z1[4] = {a.z, b.z, c.z, d.z};	__declspec(align(16)) float vectresult[4];	__m128 m1, m2, m3, m4, m0_x, m0_y, m0_z;	__m128 pSrc1 =  _mm_load_ps(&x1[0]);	__m128 pSrc2 =  _mm_load_ps(&y1[0]);	__m128 pSrc3 =  _mm_load_ps(&z1[0]);	__m128* pOut = (__m128*) vectresult;	float ae, be, ce, de;	Vector e;	printf("SSE\n");	tic();	while(count < 1000000){		e = v[count];		m0_x = _mm_set_ps1(e.x);		m0_y = _mm_set_ps1(e.y); 		m0_z = _mm_set_ps1(e.z);  		m1 = _mm_mul_ps(pSrc1, m0_x); 		m2 = _mm_mul_ps(pSrc2, m0_y); 		m3 = _mm_mul_ps(pSrc3, m0_z);		m4 = _mm_add_ps(m1, m2);		m4 = _mm_add_ps(m4, m3);		*pOut = m4;		count++;	}	toc();	count = 0;	printf("Normal\n");	tic();	while(count < 1000000){		e = v[count];		ae = (a.x*e.x) + (a.y*e.y) + (a.z*e.z);		be = (b.x*e.x) + (b.y*e.y) + (b.z*e.z);		ce = (c.x*e.x) + (c.y*e.y) + (c.z*e.z);		de = (d.x*e.x) + (d.y*e.y) + (d.z*e.z);		count++;	}	toc();	count = 0;	printf("Matrix\n");	Vector result;	tic();	while(count < 1000000){		transformVector1(mat, &v[count], &result, zAssumption);		count++;	}	toc();	//MatrixMultiply3	count = 0;	printf("SSE Matrix1\n");	tic();BatchTransform1(matCol, ∈[0][0], &out[0][0],1000000);	/*while(count < 1000000){		MatrixMultiply3(matCol, ∈[0], &out[0]);		count++;	}*/	toc();	std::cin.get(); //wait for keypress	return 0;}

Erik Rufelt

Apr 11, 2009 at 8:05 AM

Do you check the generated assembly?
It might be possible that the compiler doesn't inline the asm-function or something like that..

DarkMatter2008

Apr 11, 2009 at 8:47 AM

I've put inline declrations on the SSE functions and told the compiler to inline fucntions using that declarations, there's no calls to these functions in the assmebly now. Has reduced the result but not enough, these are the times now:-

SSE
processing time: 49 ms
Normal
processing time: 22 ms
Matrix
processing time: 33 ms
SSE Matrix1
processing time: 1804 ms

Obviousely I'm probably missing something in my knowledge here because that last time is way off when compaired to the traditional way. Can anybody spot what I am doing wrong?

HellRaiZer

Apr 11, 2009 at 10:18 AM

I think the comparison is unfair. The "normal" way of doing things isn't the same as the SSE one. Two reasons:
1) You keep the result of all the transformations to the same variable. This can be optimized away by the compiler if it's clever enough. And if it's not, it will certainly help the cache.
2) The 'result' variable is kept on the stack and it's not used by the code following the loop, so again it can be optimized away.

Why don't you write the transformVector1() function in the same way BatchTransform1() works?
a) Get the z coord from the input vector.
b) Use the out[] array in order to keep the results.

In order to be sure than the actual work isn't omitted by the compiler, print all the values to a file after you are done (outside of the tic()/toc() blocks, of course).

Finally, using intrinsics instead of inline asm, can let the compiler optimize the function when it's inlined.

HellRaiZer

Krypt0n

Apr 11, 2009 at 11:42 PM

SSE for just one dot was really useless, it was slower.
the new dot product instruction (_mm_dp_ps) really boost it by eons.

but it's still not worth if you have to prepare for a longer time than the fpu dots would take. e.g. _mm_set_ps1 is horribly slow. to get a benefit keep all your vectors as __m128

typedef struct VECTOR{  __m128 V;}Vector;

to keep some compatibility with the non-sse version, make it a class, implement accessor functions.

if you do this, you can get your matrix44*vector4 done in 4muls+3adds. compared to 16mul+12adds with FPU that will give you a boost for sure.

samoth

Apr 12, 2009 at 6:33 AM

Quote:
Original post by Krypt0n
SSE for just one dot was really useless, it was slower.
the new dot product instruction (_mm_dp_ps) really boost it by eons.

Not my experience.
In my opinion, SSE is certainly better. It sure does not run 4 times as fast as one might hope, but a careful implementation of dot product is still about 15-20% faster than plain C (not that it matters, because the few thousand dot products you usually do probably aren't a bottleneck anyway).

The real reason why using SSE when you can is better, however, is that it may pair very nicely with non-SSE instructions. For example, I've written an OFB coder in SSE to encode network packets. Later, I added a checksum calculation in C because I thought, well, one could always add a MAC if it turns out that it isn't computionally prohibitive (after all, you don't want to burn too many cycles on sorting out your network packets). So I ran both versions and got a big surprise.
It turned out that the compiler coalesces the instructions for calculating the current block's key into the instructions to update the checksum on the last one, and the CPU is able to perform them mostly in parallel.
With "encode" and "checksum" being about equally fast alone, "encode+checksum" ran about 20% slower than "encode" alone on Intel. On AMD, it took exactly the same time with or without checksum, which is quite cool.

The same may happen if you code your dot product in SSE. With any luck, it might happen to be entirely free thanks to a clever compiler.

As of DPSS, I've measured it roughly 2 cycles faster than a SSE implementation without, so the gain is not that terribly amazing. Seeing that it requires checking for SSE 4.1 and alternate code paths (which honestly is a pain), I didn't deem it worthwile.

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.