英特尔多核平台编码优化大赛的优化过程--补充
HouSisong@GMail.com 2007.01.20
tag: 多核编程,sse2,牛顿迭代,代码优化,优化大赛,invsqrt,开方
主要文章请参看我的《英特尔多核平台编码优化大赛的优化过程》： http://blog.csdn.net/housisong/archive/2007/01/20/1488465.aspx
本文章是其补充；提供一个完整的float实现版本、double到float的手工转换、手工得到invSqrt的粗略起始迭代值 等其它几个不算成功的实现；
我测试和优化过程中用的 CPU：AMD64x2 3600+ (双核CPU)
操作系统：Windows XP 32bit
编译器：Visual Studio 2005
大赛公布的原始代码执行时间 3.97秒
一个float完整实现版本(牺牲了计算精度)，源代码如下：
(如果用汇编实现应该还可以把速度提高一些，或者使用ICC编译器:)
/* compute the potential energy of a collection of */
/* particles interacting via pairwise potential */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>
#include <xmmintrin.h>//使用SSE1
#include <process.h>
#include <vector>
////
// 作者：侯思松 HouSisong@263.net
//计算结果的精度为float单精度浮点 版本
#define _IS_FAST
//以牺牲精度的方式加快速度，否则就运算达到float单精度
////
//#define _NEW_TIME_CLOCK
#ifdef _NEW_TIME_CLOCK
#define clock_t double
double CLOCKS_PER_SEC=0.0 ;
inline double clock() {
__int64 result;
if (CLOCKS_PER_SEC==0 )
{
QueryPerformanceFrequency((LARGE_INTEGER *)& result);
CLOCKS_PER_SEC=(double )result;
}
QueryPerformanceCounter((LARGE_INTEGER *)& result);
return (double )result;
}
#else
#include <time.h>
#endif
#define _IS_USES_MY_RAND
//单线程执行rand函数，所以使用自定义rand是安全的
const long DefthreadCount=2; //1,2,4,8,16,.. 把计算任务分成多个任务并行执行
float& m128_value(__m128& x,const long index) { return ((float*)(& x))[index]; }
#define NPARTS 1000
#define NITER 201
#define DIMS 3
#ifdef _IS_USES_MY_RAND
class CMyRand
{
private :
unsigned long _my_holdrand;
public :
CMyRand():_my_holdrand(1 ){}
inline int _my_rand (void )
{
unsigned long result=_my_holdrand * 214013 + 2531011 ;
_my_holdrand = result;
return ( (result>>16) & RAND_MAX );
}
};
CMyRand _MyRand;
inline int _my_rand (void){ return _MyRand._my_rand(); }
#else
#define _my_rand rand
#endif
int rand( void );
int computePot();
void initPositions(void );
void updatePositions(void );
__declspec(align(16)) float r[DIMS][(NPARTS+3)/4*4]; //16byte对齐
double pot;
int main() {
int i;
clock_t start, stop;
//char ctmp; std::cin>>ctmp;
initPositions();
updatePositions();
start = clock();
for( i=0; i<NITER; i++ ) {
pot = 0.0 ;
computePot();
if (i%10 == 0 )
{
printf("%5d: Potential: %20.7f " , i, pot);}
updatePositions();
}
pot = 0.0 ;
stop= clock();
printf ("Seconds = %10.9f ",(double)(stop-start)/ CLOCKS_PER_SEC);
return 0 ;
}
void initPositions() {
int i, j;
for( i=0; i<DIMS; i++ ){
for( j=0; j<NPARTS; j++ )
r[i][j]= (float)( 0.5 + _my_rand() *(1.0/ RAND_MAX) );
}
}
void updatePositions() {
int i,j;
for (i=0;i<DIMS;++ i)
{
for( j=0; j<NPARTS; ++ j )
{
r[i][j] -=(float)( 0.5 + _my_rand() *(1.0/ RAND_MAX) );
}
}
}
struct TWorkData
{
long iBegin;
long iEnd;
double fResult;
};
const __m128 xmms1_5={ (1.5),(1.5),(1.5),(1.5 ) };
const __m128 xmms_0_5={ (-0.5),(-0.5),(-0.5),(-0.5 ) };
void computePotPart_forj(int i,__m128* pResult)
{
__m128 lcpot= _mm_setzero_ps();
__m128 _mmi0=_mm_set1_ps(-r[0 ][i]);
__m128 _mmi1=_mm_set1_ps(-r[1 ][i]);
__m128 _mmi2=_mm_set1_ps(-r[2 ][i]);
int j=0 ;
//for( j=0; j<i-1; j++ ) { // "j<i-1"比较奇怪，疑为"j<i" ！
//* 把这个循环做4次展开
for(;j+4<i;j+=4 )
{
__m128 sm0=_mm_add_ps(*(__m128*)&r[0 ][j],_mmi0);
sm0= _mm_mul_ps(sm0,sm0);
__m128 sm1=_mm_add_ps(*(__m128*)&r[1 ][j],_mmi1);
sm1= _mm_mul_ps(sm1,sm1);
__m128 sm2=_mm_add_ps(*(__m128*)&r[2 ][j],_mmi2);
sm2= _mm_mul_ps(sm2,sm2);
sm0= _mm_add_ps(sm0,sm1);
sm0= _mm_add_ps(sm0,sm2);
sm1=_mm_rsqrt_ps(sm0); // 1/sqrt(,,,)
#ifndef _IS_FAST
// 牛顿迭代，提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
sm0=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0= _mm_mul_ps(sm0,sm1);
sm0= _mm_mul_ps(sm0,sm1);
sm0= _mm_add_ps(sm0,xmms1_5);
sm0= _mm_mul_ps(sm0,sm1);
#else
sm0 = sm1;
#endif
lcpot = _mm_add_ps(lcpot,sm0);
}//*/
for(;j<i-1;++ j)
{
__m128 sm0=_mm_set_ss(r[0 ][j]);
sm0= _mm_add_ss(sm0,_mmi0);
sm0= _mm_mul_ss(sm0,sm0);
__m128 sm1=_mm_set_ss(r[1 ][j]);
sm1= _mm_add_ss(sm1,_mmi1);
sm1= _mm_mul_ss(sm1,sm1);
__m128 sm2=_mm_set_ss(r[2 ][j]);
sm2= _mm_add_ss(sm2,_mmi2);
sm2= _mm_mul_ss(sm2,sm2);
sm0= _mm_add_ss(sm0,sm1);
sm0= _mm_add_ss(sm0,sm2);
sm1= _mm_rsqrt_ss(sm0);
#ifndef _IS_FAST
// 牛顿迭代，提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
sm0=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0= _mm_mul_ps(sm0,sm1);
sm0= _mm_mul_ps(sm0,sm1);
sm0= _mm_add_ps(sm0,xmms1_5);
sm0= _mm_mul_ps(sm0,sm1);
#else
sm0 = sm1;
#endif
lcpot = _mm_add_ss(lcpot,sm0);
}
*pResult=_mm_add_ps(* pResult,lcpot);
}
void computePotPart(TWorkData* work_data) {
int i;
__m128 lcpot= _mm_setzero_ps();
// #pragma omp parallel for schedule(static)
for( i=work_data->iBegin; i<work_data->iEnd; i++ ) {
computePotPart_forj(i,& lcpot);
}
__m128 dt0;
dt0= _mm_movehl_ps(lcpot,lcpot);
lcpot= _mm_add_ps(lcpot,dt0);
dt0=_mm_shuffle_ps(lcpot,lcpot,1 );
lcpot= _mm_add_ss(lcpot,dt0);
work_data->fResult=m128_value(lcpot,0 );
}
/
// 工作线程池 TWorkThreadPool
// 用于把一个任务拆分成多个线程任务
// 要求每个小任务任务量相近
//todo:改成任务领取模式
class TWorkThreadPool;
typedef void (*TThreadCallBack)(void * pData);
enum TThreadState{ thrStartup=0 , thrReady, thrBusy, thrTerminate, thrDeath };
class TWorkThread
{
public :
volatile HANDLE thread_handle;
volatile enum TThreadState state;
volatile TThreadCallBack func;
volatile void * pdata; //work data
volatile HANDLE waitfor_event;
TWorkThreadPool* pool;
volatile DWORD thread_ThreadAffinityMask;
TWorkThread() { memset(this,0,sizeof (TWorkThread)); }
};
void do_work_end(TWorkThread* thread_data);
void __cdecl thread_dowork(TWorkThread* thread_data) //void __stdcall thread_dowork(TWorkThread* thread_data)
{
volatile TThreadState& state=thread_data-> state;
SetThreadAffinityMask(GetCurrentThread(),thread_data-> thread_ThreadAffinityMask);
state = thrStartup;
while(true )
{
WaitForSingleObject(thread_data->waitfor_event, -1 );
if(state == thrTerminate)
break ;
state = thrBusy;
volatile TThreadCallBack& func=thread_data-> func;
if (func!=0 )
func((void *)thread_data-> pdata);
do_work_end(thread_data);
}
state = thrDeath;
_endthread();//ExitThread(0);
}
class TWorkThreadPool
{
private :
volatile HANDLE thread_event;
volatile HANDLE new_thread_event;
std::vector<TWorkThread> work_threads;
inline int passel_count() const { return (int)work_threads.size()+1 ; }
void inti_threads(long work_count) {
SYSTEM_INFO SystemInfo;
GetSystemInfo(& SystemInfo);
long cpu_count = SystemInfo.dwNumberOfProcessors;
long best_count = cpu_count;
if (cpu_count>work_count) best_count= work_count;
long newthrcount=best_count - 1 ;
work_threads.resize(newthrcount);
thread_event = CreateSemaphore(0, 0,newthrcount , 0 );
new_thread_event = CreateSemaphore(0, 0,newthrcount , 0 );
for(long i = 0; i < newthrcount; ++ i)
{
work_threads[i].waitfor_event= thread_event;
work_threads[i].state = thrTerminate;
work_threads[i].pool=this ;
work_threads[i].thread_ThreadAffinityMask=1<<(i+1 );
//DWORD thr_id;
work_threads[i].thread_handle =(HANDLE)_beginthread((void (__cdecl *)(void *))thread_dowork, 0, (void*)& work_threads[i]);
//CreateThread(0, 0, (LPTHREAD_START_ROUTINE)thread_dowork,(void*) &work_threads[i], 0, &thr_id);
}
SetThreadAffinityMask(GetCurrentThread(),0x01 );
for(long i = 0; i < newthrcount; ++ i)
{
while(true ) {
if (work_threads[i].state == thrStartup) break ;
else Sleep(0 );
}
work_threads[i].state = thrReady;
}
}
void free_threads(void )
{
long thr_count=(long )work_threads.size();
long i;
for(i = 0; i <thr_count; ++ i)
{
while(true ) {
if (work_threads[i].state == thrReady) break ;
else Sleep(0 );
}
work_threads[i].state= thrTerminate;
}
if (thr_count>0 )
ReleaseSemaphore(thread_event,thr_count, 0 );
for(i = 0; i <thr_count; ++ i)
{
while(true ) {
if (work_threads[i].state == thrDeath) break ;
else Sleep(0 );
}
}
CloseHandle(thread_event);
CloseHandle(new_thread_event);
work_threads.clear();
}
void passel_work(TThreadCallBack work_proc,void** word_data_list,int work_count) {
// assert(work_count>=1);
//assert(work_count<=passel_count());
if (work_count==1 )
{
work_proc(word_data_list[0 ]);
}
else
{
long i;
long thr_count=(long )work_threads.size();
for(i = 0; i < work_count-1; ++ i)
{
work_threads[i].func = work_proc;
work_threads[i].pdata = word_data_list[i];
work_threads[i].state = thrBusy;
}
for(i = work_count-1; i < thr_count; ++ i)
{
work_threads[i].func = 0 ;
work_threads[i].pdata =0 ;
work_threads[i].state = thrBusy;
}
if (thr_count>0 )
ReleaseSemaphore(thread_event,thr_count, 0 );
//current thread do a work
work_proc(word_data_list[work_count-1 ]);
//wait for work finish
for(i = 0; i <thr_count; ++ i)
{
while(true ) {
if (work_threads[i].state == thrReady) break ;
else Sleep(0 );
}
}
std::swap(thread_event,new_thread_event);
}
}
public :
explicit TWorkThreadPool(unsigned long work_count):thread_event(0 ),work_threads() {
//assert(work_count>=1);
inti_threads(work_count); }
~ TWorkThreadPool() { free_threads(); }
long best_work_count() const { return passel_count(); }
void work_execute(TThreadCallBack work_proc,void** word_data_list,int work_count) {
while (work_count>0 )
{
long passel_work_count;
if (work_count>= passel_count())
passel_work_count= passel_count();
else
passel_work_count = work_count;
passel_work(work_proc,word_data_list,passel_work_count);
work_count-= passel_work_count;
word_data_list=& word_data_list[passel_work_count];
}
}
inline void DoWorkEnd(TWorkThread* thread_data){
thread_data->waitfor_event= new_thread_event;
thread_data->func=0 ;
thread_data->state = thrReady;
}
};
void do_work_end(TWorkThread* thread_data)
{
thread_data->pool-> DoWorkEnd(thread_data);
}
//TWorkThreadPool end;
////
static TWorkThreadPool g_work_thread_pool(DefthreadCount);//线程池
int computePot() {
static bool is_inti_work=false ;
static TWorkData work_list[DefthreadCount];
static TWorkData* pwork_list[DefthreadCount];
int i;
if (! is_inti_work)
{
long fi=0 ;
for (int i=0;i<DefthreadCount;++ i)
{
if (0== i)
work_list[i].iBegin=0 ;
else
work_list[i].iBegin =work_list[i-1 ].iEnd;
if (i==DefthreadCount-1 )
work_list[i].iEnd=(long )NPARTS;
else
work_list[i].iEnd =(long)( (double)(NPARTS-1)*sqrt((double)(i+1)/DefthreadCount)+1+0.5 );
pwork_list[i]=& work_list[i];
}
is_inti_work=true ;
}
g_work_thread_pool.work_execute((TThreadCallBack)computePotPart,(void **)(&pwork_list[0 ]),DefthreadCount);
for (i=0;i<DefthreadCount;++ i)
pot+= work_list[i].fResult;
return 0 ;
}
/* particles interacting via pairwise potential */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <windows.h>
#include <xmmintrin.h>//使用SSE1
#include <process.h>
#include <vector>
////
// 作者：侯思松 HouSisong@263.net
//计算结果的精度为float单精度浮点 版本
#define _IS_FAST
//以牺牲精度的方式加快速度，否则就运算达到float单精度
////
//#define _NEW_TIME_CLOCK
#ifdef _NEW_TIME_CLOCK
#define clock_t double
double CLOCKS_PER_SEC=0.0 ;
inline double clock() {
__int64 result;
if (CLOCKS_PER_SEC==0 )
{
QueryPerformanceFrequency((LARGE_INTEGER *)& result);
CLOCKS_PER_SEC=(double )result;
}
QueryPerformanceCounter((LARGE_INTEGER *)& result);
return (double )result;
}
#else
#include <time.h>
#endif
#define _IS_USES_MY_RAND
//单线程执行rand函数，所以使用自定义rand是安全的
const long DefthreadCount=2; //1,2,4,8,16,.. 把计算任务分成多个任务并行执行
float& m128_value(__m128& x,const long index) { return ((float*)(& x))[index]; }
#define NPARTS 1000
#define NITER 201
#define DIMS 3
#ifdef _IS_USES_MY_RAND
class CMyRand
{
private :
unsigned long _my_holdrand;
public :
CMyRand():_my_holdrand(1 ){}
inline int _my_rand (void )
{
unsigned long result=_my_holdrand * 214013 + 2531011 ;
_my_holdrand = result;
return ( (result>>16) & RAND_MAX );
}
};
CMyRand _MyRand;
inline int _my_rand (void){ return _MyRand._my_rand(); }
#else
#define _my_rand rand
#endif
int rand( void );
int computePot();
void initPositions(void );
void updatePositions(void );
__declspec(align(16)) float r[DIMS][(NPARTS+3)/4*4]; //16byte对齐
double pot;
int main() {
int i;
clock_t start, stop;
//char ctmp; std::cin>>ctmp;
initPositions();
updatePositions();
start = clock();
for( i=0; i<NITER; i++ ) {
pot = 0.0 ;
computePot();
if (i%10 == 0 )
{
printf("%5d: Potential: %20.7f " , i, pot);}
updatePositions();
}
pot = 0.0 ;
stop= clock();
printf ("Seconds = %10.9f ",(double)(stop-start)/ CLOCKS_PER_SEC);
return 0 ;
}
void initPositions() {
int i, j;
for( i=0; i<DIMS; i++ ){
for( j=0; j<NPARTS; j++ )
r[i][j]= (float)( 0.5 + _my_rand() *(1.0/ RAND_MAX) );
}
}
void updatePositions() {
int i,j;
for (i=0;i<DIMS;++ i)
{
for( j=0; j<NPARTS; ++ j )
{
r[i][j] -=(float)( 0.5 + _my_rand() *(1.0/ RAND_MAX) );
}
}
}
struct TWorkData
{
long iBegin;
long iEnd;
double fResult;
};
const __m128 xmms1_5={ (1.5),(1.5),(1.5),(1.5 ) };
const __m128 xmms_0_5={ (-0.5),(-0.5),(-0.5),(-0.5 ) };
void computePotPart_forj(int i,__m128* pResult)
{
__m128 lcpot= _mm_setzero_ps();
__m128 _mmi0=_mm_set1_ps(-r[0 ][i]);
__m128 _mmi1=_mm_set1_ps(-r[1 ][i]);
__m128 _mmi2=_mm_set1_ps(-r[2 ][i]);
int j=0 ;
//for( j=0; j<i-1; j++ ) { // "j<i-1"比较奇怪，疑为"j<i" ！
//* 把这个循环做4次展开
for(;j+4<i;j+=4 )
{
__m128 sm0=_mm_add_ps(*(__m128*)&r[0 ][j],_mmi0);
sm0= _mm_mul_ps(sm0,sm0);
__m128 sm1=_mm_add_ps(*(__m128*)&r[1 ][j],_mmi1);
sm1= _mm_mul_ps(sm1,sm1);
__m128 sm2=_mm_add_ps(*(__m128*)&r[2 ][j],_mmi2);
sm2= _mm_mul_ps(sm2,sm2);
sm0= _mm_add_ps(sm0,sm1);
sm0= _mm_add_ps(sm0,sm2);
sm1=_mm_rsqrt_ps(sm0); // 1/sqrt(,,,)
#ifndef _IS_FAST
// 牛顿迭代，提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
sm0=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0= _mm_mul_ps(sm0,sm1);
sm0= _mm_mul_ps(sm0,sm1);
sm0= _mm_add_ps(sm0,xmms1_5);
sm0= _mm_mul_ps(sm0,sm1);
#else
sm0 = sm1;
#endif
lcpot = _mm_add_ps(lcpot,sm0);
}//*/
for(;j<i-1;++ j)
{
__m128 sm0=_mm_set_ss(r[0 ][j]);
sm0= _mm_add_ss(sm0,_mmi0);
sm0= _mm_mul_ss(sm0,sm0);
__m128 sm1=_mm_set_ss(r[1 ][j]);
sm1= _mm_add_ss(sm1,_mmi1);
sm1= _mm_mul_ss(sm1,sm1);
__m128 sm2=_mm_set_ss(r[2 ][j]);
sm2= _mm_add_ss(sm2,_mmi2);
sm2= _mm_mul_ss(sm2,sm2);
sm0= _mm_add_ss(sm0,sm1);
sm0= _mm_add_ss(sm0,sm2);
sm1= _mm_rsqrt_ss(sm0);
#ifndef _IS_FAST
// 牛顿迭代，提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
sm0=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0= _mm_mul_ps(sm0,sm1);
sm0= _mm_mul_ps(sm0,sm1);
sm0= _mm_add_ps(sm0,xmms1_5);
sm0= _mm_mul_ps(sm0,sm1);
#else
sm0 = sm1;
#endif
lcpot = _mm_add_ss(lcpot,sm0);
}
*pResult=_mm_add_ps(* pResult,lcpot);
}
void computePotPart(TWorkData* work_data) {
int i;
__m128 lcpot= _mm_setzero_ps();
// #pragma omp parallel for schedule(static)
for( i=work_data->iBegin; i<work_data->iEnd; i++ ) {
computePotPart_forj(i,& lcpot);
}
__m128 dt0;
dt0= _mm_movehl_ps(lcpot,lcpot);
lcpot= _mm_add_ps(lcpot,dt0);
dt0=_mm_shuffle_ps(lcpot,lcpot,1 );
lcpot= _mm_add_ss(lcpot,dt0);
work_data->fResult=m128_value(lcpot,0 );
}
/
// 工作线程池 TWorkThreadPool
// 用于把一个任务拆分成多个线程任务
// 要求每个小任务任务量相近
//todo:改成任务领取模式
class TWorkThreadPool;
typedef void (*TThreadCallBack)(void * pData);
enum TThreadState{ thrStartup=0 , thrReady, thrBusy, thrTerminate, thrDeath };
class TWorkThread
{
public :
volatile HANDLE thread_handle;
volatile enum TThreadState state;
volatile TThreadCallBack func;
volatile void * pdata; //work data
volatile HANDLE waitfor_event;
TWorkThreadPool* pool;
volatile DWORD thread_ThreadAffinityMask;
TWorkThread() { memset(this,0,sizeof (TWorkThread)); }
};
void do_work_end(TWorkThread* thread_data);
void __cdecl thread_dowork(TWorkThread* thread_data) //void __stdcall thread_dowork(TWorkThread* thread_data)
{
volatile TThreadState& state=thread_data-> state;
SetThreadAffinityMask(GetCurrentThread(),thread_data-> thread_ThreadAffinityMask);
state = thrStartup;
while(true )
{
WaitForSingleObject(thread_data->waitfor_event, -1 );
if(state == thrTerminate)
break ;
state = thrBusy;
volatile TThreadCallBack& func=thread_data-> func;
if (func!=0 )
func((void *)thread_data-> pdata);
do_work_end(thread_data);
}
state = thrDeath;
_endthread();//ExitThread(0);
}
class TWorkThreadPool
{
private :
volatile HANDLE thread_event;
volatile HANDLE new_thread_event;
std::vector<TWorkThread> work_threads;
inline int passel_count() const { return (int)work_threads.size()+1 ; }
void inti_threads(long work_count) {
SYSTEM_INFO SystemInfo;
GetSystemInfo(& SystemInfo);
long cpu_count = SystemInfo.dwNumberOfProcessors;
long best_count = cpu_count;
if (cpu_count>work_count) best_count= work_count;
long newthrcount=best_count - 1 ;
work_threads.resize(newthrcount);
thread_event = CreateSemaphore(0, 0,newthrcount , 0 );
new_thread_event = CreateSemaphore(0, 0,newthrcount , 0 );
for(long i = 0; i < newthrcount; ++ i)
{
work_threads[i].waitfor_event= thread_event;
work_threads[i].state = thrTerminate;
work_threads[i].pool=this ;
work_threads[i].thread_ThreadAffinityMask=1<<(i+1 );
//DWORD thr_id;
work_threads[i].thread_handle =(HANDLE)_beginthread((void (__cdecl *)(void *))thread_dowork, 0, (void*)& work_threads[i]);
//CreateThread(0, 0, (LPTHREAD_START_ROUTINE)thread_dowork,(void*) &work_threads[i], 0, &thr_id);
}
SetThreadAffinityMask(GetCurrentThread(),0x01 );
for(long i = 0; i < newthrcount; ++ i)
{
while(true ) {
if (work_threads[i].state == thrStartup) break ;
else Sleep(0 );
}
work_threads[i].state = thrReady;
}
}
void free_threads(void )
{
long thr_count=(long )work_threads.size();
long i;
for(i = 0; i <thr_count; ++ i)
{
while(true ) {
if (work_threads[i].state == thrReady) break ;
else Sleep(0 );
}
work_threads[i].state= thrTerminate;
}
if (thr_count>0 )
ReleaseSemaphore(thread_event,thr_count, 0 );
for(i = 0; i <thr_count; ++ i)
{
while(true ) {
if (work_threads[i].state == thrDeath) break ;
else Sleep(0 );
}
}
CloseHandle(thread_event);
CloseHandle(new_thread_event);
work_threads.clear();
}
void passel_work(TThreadCallBack work_proc,void** word_data_list,int work_count) {
// assert(work_count>=1);
//assert(work_count<=passel_count());
if (work_count==1 )
{
work_proc(word_data_list[0 ]);
}
else
{
long i;
long thr_count=(long )work_threads.size();
for(i = 0; i < work_count-1; ++ i)
{
work_threads[i].func = work_proc;
work_threads[i].pdata = word_data_list[i];
work_threads[i].state = thrBusy;
}
for(i = work_count-1; i < thr_count; ++ i)
{
work_threads[i].func = 0 ;
work_threads[i].pdata =0 ;
work_threads[i].state = thrBusy;
}
if (thr_count>0 )
ReleaseSemaphore(thread_event,thr_count, 0 );
//current thread do a work
work_proc(word_data_list[work_count-1 ]);
//wait for work finish
for(i = 0; i <thr_count; ++ i)
{
while(true ) {
if (work_threads[i].state == thrReady) break ;
else Sleep(0 );
}
}
std::swap(thread_event,new_thread_event);
}
}
public :
explicit TWorkThreadPool(unsigned long work_count):thread_event(0 ),work_threads() {
//assert(work_count>=1);
inti_threads(work_count); }
~ TWorkThreadPool() { free_threads(); }
long best_work_count() const { return passel_count(); }
void work_execute(TThreadCallBack work_proc,void** word_data_list,int work_count) {
while (work_count>0 )
{
long passel_work_count;
if (work_count>= passel_count())
passel_work_count= passel_count();
else
passel_work_count = work_count;
passel_work(work_proc,word_data_list,passel_work_count);
work_count-= passel_work_count;
word_data_list=& word_data_list[passel_work_count];
}
}
inline void DoWorkEnd(TWorkThread* thread_data){
thread_data->waitfor_event= new_thread_event;
thread_data->func=0 ;
thread_data->state = thrReady;
}
};
void do_work_end(TWorkThread* thread_data)
{
thread_data->pool-> DoWorkEnd(thread_data);
}
//TWorkThreadPool end;
////
static TWorkThreadPool g_work_thread_pool(DefthreadCount);//线程池
int computePot() {
static bool is_inti_work=false ;
static TWorkData work_list[DefthreadCount];
static TWorkData* pwork_list[DefthreadCount];
int i;
if (! is_inti_work)
{
long fi=0 ;
for (int i=0;i<DefthreadCount;++ i)
{
if (0== i)
work_list[i].iBegin=0 ;
else
work_list[i].iBegin =work_list[i-1 ].iEnd;
if (i==DefthreadCount-1 )
work_list[i].iEnd=(long )NPARTS;
else
work_list[i].iEnd =(long)( (double)(NPARTS-1)*sqrt((double)(i+1)/DefthreadCount)+1+0.5 );
pwork_list[i]=& work_list[i];
}
is_inti_work=true ;
}
g_work_thread_pool.work_execute((TThreadCallBack)computePotPart,(void **)(&pwork_list[0 ]),DefthreadCount);
for (i=0;i<DefthreadCount;++ i)
pot+= work_list[i].fResult;
return 0 ;
}
代码执行时间 0.125秒 相对于原始代码加速比：3176.0%
注意到一个事实，float比double版快出了很多，why?
原来，double版中为了使用SSE而在float和double之间的转换花费了很多的时间！
我不知道这个问题是AMD64x2 CPU的问题还是在酷睿2上也一样;
double版中为了优化这个转换，我预先保存一份r数组的float转化值计算，这样就能节约double到float转换；
(double版完整源代码参见《英特尔多核平台编码优化大赛的优化过程》)
定义一个临时数组rf:
__declspec(align(16)) float rf[DIMS][(NPARTS+3)/4*4];
修改updatePositions函数:
void
updatePositions() {
int i,j;
for (i=0;i<DIMS;++ i)
{
for( j=0; j<NPARTS; ++ j )
{
r[i][j] -=( 0.5 + _my_rand() *(1.0/ RAND_MAX) );
rf[i][j]=(float )r[i][j];
}
}
}
int i,j;
for (i=0;i<DIMS;++ i)
{
for( j=0; j<NPARTS; ++ j )
{
r[i][j] -=( 0.5 + _my_rand() *(1.0/ RAND_MAX) );
rf[i][j]=(float )r[i][j];
}
}
}
然后重新实现computePotPart_forj：
void computePotPart_forj_double_float(int i,__m128d*
pResult)
{
__m128d lcpot= _mm_setzero_pd();
__m128d _mmi0=_mm_set1_pd(-r[0 ][i]);
__m128d _mmi1=_mm_set1_pd(-r[1 ][i]);
__m128d _mmi2=_mm_set1_pd(-r[2 ][i]);
__m128 _mmi0f=_mm_set1_ps(-rf[0 ][i]);
__m128 _mmi1f=_mm_set1_ps(-rf[1 ][i]);
__m128 _mmi2f=_mm_set1_ps(-rf[2 ][i]);
int j=0 ;
//*
for(;j+4<i;j+=4 )
{
__m128d dm0=_mm_add_pd(*(__m128d*)&r[0 ][j],_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_add_pd(*(__m128d*)&r[1 ][j],_mmi1);
dm1= _mm_mul_pd(dm1,dm1);
__m128d dm2=_mm_add_pd(*(__m128d*)&r[2 ][j],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm0= _mm_add_pd(dm0,dm1);
dm0= _mm_add_pd(dm0,dm2);
__m128d dm5=_mm_add_pd(*(__m128d*)&r[0][j+2 ],_mmi0);
dm5= _mm_mul_pd(dm5,dm5);
__m128d dm6=_mm_add_pd(*(__m128d*)&r[1][j+2 ],_mmi1);
dm6= _mm_mul_pd(dm6,dm6);
dm2=_mm_add_pd(*(__m128d*)&r[2][j+2 ],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm5= _mm_add_pd(dm5,dm6);
dm5= _mm_add_pd(dm5,dm2);
// 用SSE的rsqrt近似计算1/sqrt(a)； 然后使用牛顿迭代来提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
{
/* __m128 sm0=_mm_cvtpd_ps(dm0);
__m128 sm1=_mm_cvtpd_ps(dm5);
sm0=_mm_movelh_ps(sm0,sm1);*/
__m128 sm0 =_mm_add_ps(*(__m128*)&rf[0 ][j],_mmi0f);
sm0= _mm_mul_ps(sm0,sm0);
__m128 sm1=_mm_add_ps(*(__m128*)&rf[1 ][j],_mmi1f);
sm1= _mm_mul_ps(sm1,sm1);
__m128 sm2=_mm_add_ps(*(__m128*)&rf[2 ][j],_mmi2f);
sm2= _mm_mul_ps(sm2,sm2);
sm0= _mm_add_ps(sm0,sm1);
sm0= _mm_add_ps(sm0,sm2);
__m128 sma=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ps(sm0); // 计算1/sqrt(a)
//牛顿迭代，提高开方精度
sma= _mm_mul_ps(sma,sm0);
sma= _mm_mul_ps(sma,sm0);
sma= _mm_add_ps(sma,xmms1_5);
sm0= _mm_mul_ps(sm0,sma);
__m128d dma=_mm_mul_pd(dm0,xmmd_0_5); //a*(-0.5)
__m128d dmb= _mm_mul_pd(dm5,xmmd_0_5);
sm1= _mm_movehl_ps(sm1,sm0);
dm0=_mm_cvtps_pd(sm0); //
dm5=_mm_cvtps_pd(sm1); //
//再次迭代，加倍精度
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_add_pd(dma,xmmd1_5);
dmb= _mm_add_pd(dmb,xmmd1_5);
dm0= _mm_mul_pd(dm0,dma);
dm5= _mm_mul_pd(dm5,dmb);
}
lcpot= _mm_add_pd(lcpot,dm0);
lcpot= _mm_add_pd(lcpot,dm5);
}
for (;j+1<i;++ j)
{
__m128d dm0=_mm_set_sd(r[0 ][j]);
dm0= _mm_add_pd(dm0,_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_set_sd(r[1 ][j]);
dm1= _mm_add_sd(dm1,_mmi1);
dm1= _mm_mul_sd(dm1,dm1);
__m128d dm2=_mm_set_sd(r[2 ][j]);
dm2= _mm_add_sd(dm2,_mmi2);
dm2= _mm_mul_sd(dm2,dm2);
dm0= _mm_add_sd(dm0,dm1);
dm0= _mm_add_sd(dm0,dm2);
{
__m128 sm0= _mm_cvtpd_ps(dm0);
__m128d dma=_mm_mul_sd(dm0,xmmd_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ss(sm0); //计算1/sqrt(a)
dm0=_mm_cvtps_pd(sm0); //
//牛顿迭代，提高开方精度
dm1= _mm_mul_sd(dm0,dm0);
dm1= _mm_mul_sd(dm1,dma);
dm1= _mm_add_sd(dm1,xmmd1_5);
dm0= _mm_mul_sd(dm0,dm1);
//再次迭代
dma= _mm_mul_sd(dma,dm0);
dma= _mm_mul_sd(dma,dm0);
dma= _mm_add_sd(dma,xmmd1_5);
dm0= _mm_mul_sd(dm0,dma);
}
lcpot= _mm_add_sd(lcpot,dm0);
}
*pResult=_mm_add_pd(* pResult,lcpot);
}
{
__m128d lcpot= _mm_setzero_pd();
__m128d _mmi0=_mm_set1_pd(-r[0 ][i]);
__m128d _mmi1=_mm_set1_pd(-r[1 ][i]);
__m128d _mmi2=_mm_set1_pd(-r[2 ][i]);
__m128 _mmi0f=_mm_set1_ps(-rf[0 ][i]);
__m128 _mmi1f=_mm_set1_ps(-rf[1 ][i]);
__m128 _mmi2f=_mm_set1_ps(-rf[2 ][i]);
int j=0 ;
//*
for(;j+4<i;j+=4 )
{
__m128d dm0=_mm_add_pd(*(__m128d*)&r[0 ][j],_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_add_pd(*(__m128d*)&r[1 ][j],_mmi1);
dm1= _mm_mul_pd(dm1,dm1);
__m128d dm2=_mm_add_pd(*(__m128d*)&r[2 ][j],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm0= _mm_add_pd(dm0,dm1);
dm0= _mm_add_pd(dm0,dm2);
__m128d dm5=_mm_add_pd(*(__m128d*)&r[0][j+2 ],_mmi0);
dm5= _mm_mul_pd(dm5,dm5);
__m128d dm6=_mm_add_pd(*(__m128d*)&r[1][j+2 ],_mmi1);
dm6= _mm_mul_pd(dm6,dm6);
dm2=_mm_add_pd(*(__m128d*)&r[2][j+2 ],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm5= _mm_add_pd(dm5,dm6);
dm5= _mm_add_pd(dm5,dm2);
// 用SSE的rsqrt近似计算1/sqrt(a)； 然后使用牛顿迭代来提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
{
/* __m128 sm0=_mm_cvtpd_ps(dm0);
__m128 sm1=_mm_cvtpd_ps(dm5);
sm0=_mm_movelh_ps(sm0,sm1);*/
__m128 sm0 =_mm_add_ps(*(__m128*)&rf[0 ][j],_mmi0f);
sm0= _mm_mul_ps(sm0,sm0);
__m128 sm1=_mm_add_ps(*(__m128*)&rf[1 ][j],_mmi1f);
sm1= _mm_mul_ps(sm1,sm1);
__m128 sm2=_mm_add_ps(*(__m128*)&rf[2 ][j],_mmi2f);
sm2= _mm_mul_ps(sm2,sm2);
sm0= _mm_add_ps(sm0,sm1);
sm0= _mm_add_ps(sm0,sm2);
__m128 sma=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ps(sm0); // 计算1/sqrt(a)
//牛顿迭代，提高开方精度
sma= _mm_mul_ps(sma,sm0);
sma= _mm_mul_ps(sma,sm0);
sma= _mm_add_ps(sma,xmms1_5);
sm0= _mm_mul_ps(sm0,sma);
__m128d dma=_mm_mul_pd(dm0,xmmd_0_5); //a*(-0.5)
__m128d dmb= _mm_mul_pd(dm5,xmmd_0_5);
sm1= _mm_movehl_ps(sm1,sm0);
dm0=_mm_cvtps_pd(sm0); //
dm5=_mm_cvtps_pd(sm1); //
//再次迭代，加倍精度
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_add_pd(dma,xmmd1_5);
dmb= _mm_add_pd(dmb,xmmd1_5);
dm0= _mm_mul_pd(dm0,dma);
dm5= _mm_mul_pd(dm5,dmb);
}
lcpot= _mm_add_pd(lcpot,dm0);
lcpot= _mm_add_pd(lcpot,dm5);
}
for (;j+1<i;++ j)
{
__m128d dm0=_mm_set_sd(r[0 ][j]);
dm0= _mm_add_pd(dm0,_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_set_sd(r[1 ][j]);
dm1= _mm_add_sd(dm1,_mmi1);
dm1= _mm_mul_sd(dm1,dm1);
__m128d dm2=_mm_set_sd(r[2 ][j]);
dm2= _mm_add_sd(dm2,_mmi2);
dm2= _mm_mul_sd(dm2,dm2);
dm0= _mm_add_sd(dm0,dm1);
dm0= _mm_add_sd(dm0,dm2);
{
__m128 sm0= _mm_cvtpd_ps(dm0);
__m128d dma=_mm_mul_sd(dm0,xmmd_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ss(sm0); //计算1/sqrt(a)
dm0=_mm_cvtps_pd(sm0); //
//牛顿迭代，提高开方精度
dm1= _mm_mul_sd(dm0,dm0);
dm1= _mm_mul_sd(dm1,dma);
dm1= _mm_add_sd(dm1,xmmd1_5);
dm0= _mm_mul_sd(dm0,dm1);
//再次迭代
dma= _mm_mul_sd(dma,dm0);
dma= _mm_mul_sd(dma,dm0);
dma= _mm_add_sd(dma,xmmd1_5);
dm0= _mm_mul_sd(dm0,dma);
}
lcpot= _mm_add_sd(lcpot,dm0);
}
*pResult=_mm_add_pd(* pResult,lcpot);
}
该代码对速度有很小的改进，但精度确只到小数点后5/6位(过早的降低了运算精度，后面的计算使误差放大)，
再增加一次牛顿叠代就能弥补精度的缺失，但速度上就没有优势了；只能放弃该方案；
既然硬件的float和double之间的转换慢，那我手工来写一个会不会更好一些呢？ (比较奇怪的尝试:)
见代码:
(该函数利用了double/float的IEE浮点编码结构)
const __m64 _mmd2f_esub =_mm_set_pi32((1023-127) << (52-32),0
);
const __m128i _xmmd2f_esub= _mm_set_epi64(_mmd2f_esub,_mmd2f_esub);
void computePotPart_forj_d2f(int i,__m128d* pResult)
{
__m128d lcpot= _mm_setzero_pd();
__m128d _mmi0=_mm_set1_pd(-r[0 ][i]);
__m128d _mmi1=_mm_set1_pd(-r[1 ][i]);
__m128d _mmi2=_mm_set1_pd(-r[2 ][i]);
int j=0 ;
//*
for(;j+4<i;j+=4 )
{
__m128d dm0=_mm_add_pd(*(__m128d*)&r[0 ][j],_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_add_pd(*(__m128d*)&r[1 ][j],_mmi1);
dm1= _mm_mul_pd(dm1,dm1);
__m128d dm2=_mm_add_pd(*(__m128d*)&r[2 ][j],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm0= _mm_add_pd(dm0,dm1);
dm0= _mm_add_pd(dm0,dm2);
__m128d dm5=_mm_add_pd(*(__m128d*)&r[0][j+2 ],_mmi0);
dm5= _mm_mul_pd(dm5,dm5);
__m128d dm6=_mm_add_pd(*(__m128d*)&r[1][j+2 ],_mmi1);
dm6= _mm_mul_pd(dm6,dm6);
dm2=_mm_add_pd(*(__m128d*)&r[2][j+2 ],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm5= _mm_add_pd(dm5,dm6);
dm5= _mm_add_pd(dm5,dm2);
// 用SSE的rsqrt近似计算1/sqrt(a)； 然后使用牛顿迭代来提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
{
__m128d dma=_mm_mul_pd(dm0,xmmd_0_5); //a*(-0.5)
__m128d dmb= _mm_mul_pd(dm5,xmmd_0_5);
(*(__m128i*)&dm0)=_mm_sub_epi64((*(__m128i*)& dm0),_xmmd2f_esub);
(*(__m128i*)&dm5)=_mm_sub_epi64((*(__m128i*)& dm5),_xmmd2f_esub);
(*(__m128i*)&dm0)=_mm_srli_epi64(*(__m128i*)&dm0,32-3 );
(*(__m128i*)&dm5)=_mm_srli_epi64(*(__m128i*)&dm5,32-3 );
(*(__m128i*)&dm0)=_mm_slli_epi64(*(__m128i*)&dm0,32 );
__m128 sm0;
(*(__m128i*)&sm0)=_mm_xor_si128(*(__m128i*)&dm0,*(__m128i*)& dm5);
__m128 sma=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ps(sm0); // 计算1/sqrt(a)
//牛顿迭代，提高开方精度
sma= _mm_mul_ps(sma,sm0);
sma= _mm_mul_ps(sma,sm0);
sma= _mm_add_ps(sma,xmms1_5);
sm0= _mm_mul_ps(sm0,sma);
(*(__m128i*)&dm0)=_mm_srli_epi64(*(__m128i*)&sm0,3); //
(*(__m128i*)&dm5)=_mm_slli_epi64(*(__m128i*)&sm0,32); //
(*(__m128i*)&dm5)=_mm_srli_epi64(*(__m128i*)&dm5,3); //
(*(__m128i*)&dm0)=_mm_add_epi64((*(__m128i*)& dm0),_xmmd2f_esub);
(*(__m128i*)&dm5)=_mm_add_epi64((*(__m128i*)& dm5),_xmmd2f_esub);
//再次迭代，加倍精度
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_add_pd(dma,xmmd1_5);
dmb= _mm_add_pd(dmb,xmmd1_5);
dm0= _mm_mul_pd(dm0,dma);
dm5= _mm_mul_pd(dm5,dmb);
}
lcpot= _mm_add_pd(lcpot,dm0);
lcpot= _mm_add_pd(lcpot,dm5);
}
for (;j+1<i;++ j)
{
__m128d dm0=_mm_set_sd(r[0 ][j]);
dm0= _mm_add_pd(dm0,_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_set_sd(r[1 ][j]);
dm1= _mm_add_sd(dm1,_mmi1);
dm1= _mm_mul_sd(dm1,dm1);
__m128d dm2=_mm_set_sd(r[2 ][j]);
dm2= _mm_add_sd(dm2,_mmi2);
dm2= _mm_mul_sd(dm2,dm2);
dm0= _mm_add_sd(dm0,dm1);
dm0= _mm_add_sd(dm0,dm2);
{
__m128 sm0= _mm_cvtpd_ps(dm0);
__m128d dma=_mm_mul_sd(dm0,xmmd_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ss(sm0); //计算1/sqrt(a)
dm0=_mm_cvtps_pd(sm0); //
//牛顿迭代，提高开方精度
dm1= _mm_mul_sd(dm0,dm0);
dm1= _mm_mul_sd(dm1,dma);
dm1= _mm_add_sd(dm1,xmmd1_5);
dm0= _mm_mul_sd(dm0,dm1);
//再次迭代
dma= _mm_mul_sd(dma,dm0);
dma= _mm_mul_sd(dma,dm0);
dma= _mm_add_sd(dma,xmmd1_5);
dm0= _mm_mul_sd(dm0,dma);
}
lcpot= _mm_add_sd(lcpot,dm0);
}
*pResult=_mm_add_pd(* pResult,lcpot);
}
const __m128i _xmmd2f_esub= _mm_set_epi64(_mmd2f_esub,_mmd2f_esub);
void computePotPart_forj_d2f(int i,__m128d* pResult)
{
__m128d lcpot= _mm_setzero_pd();
__m128d _mmi0=_mm_set1_pd(-r[0 ][i]);
__m128d _mmi1=_mm_set1_pd(-r[1 ][i]);
__m128d _mmi2=_mm_set1_pd(-r[2 ][i]);
int j=0 ;
//*
for(;j+4<i;j+=4 )
{
__m128d dm0=_mm_add_pd(*(__m128d*)&r[0 ][j],_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_add_pd(*(__m128d*)&r[1 ][j],_mmi1);
dm1= _mm_mul_pd(dm1,dm1);
__m128d dm2=_mm_add_pd(*(__m128d*)&r[2 ][j],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm0= _mm_add_pd(dm0,dm1);
dm0= _mm_add_pd(dm0,dm2);
__m128d dm5=_mm_add_pd(*(__m128d*)&r[0][j+2 ],_mmi0);
dm5= _mm_mul_pd(dm5,dm5);
__m128d dm6=_mm_add_pd(*(__m128d*)&r[1][j+2 ],_mmi1);
dm6= _mm_mul_pd(dm6,dm6);
dm2=_mm_add_pd(*(__m128d*)&r[2][j+2 ],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm5= _mm_add_pd(dm5,dm6);
dm5= _mm_add_pd(dm5,dm2);
// 用SSE的rsqrt近似计算1/sqrt(a)； 然后使用牛顿迭代来提高开方精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
{
__m128d dma=_mm_mul_pd(dm0,xmmd_0_5); //a*(-0.5)
__m128d dmb= _mm_mul_pd(dm5,xmmd_0_5);
(*(__m128i*)&dm0)=_mm_sub_epi64((*(__m128i*)& dm0),_xmmd2f_esub);
(*(__m128i*)&dm5)=_mm_sub_epi64((*(__m128i*)& dm5),_xmmd2f_esub);
(*(__m128i*)&dm0)=_mm_srli_epi64(*(__m128i*)&dm0,32-3 );
(*(__m128i*)&dm5)=_mm_srli_epi64(*(__m128i*)&dm5,32-3 );
(*(__m128i*)&dm0)=_mm_slli_epi64(*(__m128i*)&dm0,32 );
__m128 sm0;
(*(__m128i*)&sm0)=_mm_xor_si128(*(__m128i*)&dm0,*(__m128i*)& dm5);
__m128 sma=_mm_mul_ps(sm0,xmms_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ps(sm0); // 计算1/sqrt(a)
//牛顿迭代，提高开方精度
sma= _mm_mul_ps(sma,sm0);
sma= _mm_mul_ps(sma,sm0);
sma= _mm_add_ps(sma,xmms1_5);
sm0= _mm_mul_ps(sm0,sma);
(*(__m128i*)&dm0)=_mm_srli_epi64(*(__m128i*)&sm0,3); //
(*(__m128i*)&dm5)=_mm_slli_epi64(*(__m128i*)&sm0,32); //
(*(__m128i*)&dm5)=_mm_srli_epi64(*(__m128i*)&dm5,3); //
(*(__m128i*)&dm0)=_mm_add_epi64((*(__m128i*)& dm0),_xmmd2f_esub);
(*(__m128i*)&dm5)=_mm_add_epi64((*(__m128i*)& dm5),_xmmd2f_esub);
//再次迭代，加倍精度
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_add_pd(dma,xmmd1_5);
dmb= _mm_add_pd(dmb,xmmd1_5);
dm0= _mm_mul_pd(dm0,dma);
dm5= _mm_mul_pd(dm5,dmb);
}
lcpot= _mm_add_pd(lcpot,dm0);
lcpot= _mm_add_pd(lcpot,dm5);
}
for (;j+1<i;++ j)
{
__m128d dm0=_mm_set_sd(r[0 ][j]);
dm0= _mm_add_pd(dm0,_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_set_sd(r[1 ][j]);
dm1= _mm_add_sd(dm1,_mmi1);
dm1= _mm_mul_sd(dm1,dm1);
__m128d dm2=_mm_set_sd(r[2 ][j]);
dm2= _mm_add_sd(dm2,_mmi2);
dm2= _mm_mul_sd(dm2,dm2);
dm0= _mm_add_sd(dm0,dm1);
dm0= _mm_add_sd(dm0,dm2);
{
__m128 sm0= _mm_cvtpd_ps(dm0);
__m128d dma=_mm_mul_sd(dm0,xmmd_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ss(sm0); //计算1/sqrt(a)
dm0=_mm_cvtps_pd(sm0); //
//牛顿迭代，提高开方精度
dm1= _mm_mul_sd(dm0,dm0);
dm1= _mm_mul_sd(dm1,dma);
dm1= _mm_add_sd(dm1,xmmd1_5);
dm0= _mm_mul_sd(dm0,dm1);
//再次迭代
dma= _mm_mul_sd(dma,dm0);
dma= _mm_mul_sd(dma,dm0);
dma= _mm_add_sd(dma,xmmd1_5);
dm0= _mm_mul_sd(dm0,dma);
}
lcpot= _mm_add_sd(lcpot,dm0);
}
*pResult=_mm_add_pd(* pResult,lcpot);
}
该函数的速度比原来的代码稍慢:) 放弃之
(我还尝试过这样的代码（代码没有保存:( ）
代码原理和上面的很接近，利用IEE的浮点格式，强制把double转成float后(通过指数平衡和移位)，
分两路使用_mm_rsqrt_ps、牛顿叠代等；但这个慢慢，也放弃了)
既然硬件的float和double之间的转换慢，那我就不转换来看看怎样实现；
不使用SSE的_mm_rsqrt_ps指令，而利用IEE浮点格式生成一个粗略的近似解，然后迭代(迭代了3次)；
(也可以用查表的方式来得到初始解，但在SSE体系中，这种实现很可能得不偿失，所以没有去实现)；(该函数利用了double的IEE浮点编码结构)这样以后，速度有了一些提高，但精度还有点不够(大概有6位小数位精度)，再次迭代的话就失去了速度优势；
不知道有没有比魔法数0x5fe6ec85,0xe7de30da更好的魔法数：）
放弃;
const __m64 _mmi_mn=_mm_set_pi32(0x5fe6ec85,0xe7de30da
);
const __m128i xmmi64_mn= _mm_set1_epi64(_mmi_mn);
void computePotPart_forj_int(int i,__m128d* pResult)
{
__m128d lcpot= _mm_setzero_pd();
__m128d _mmi0=_mm_set1_pd(-r[0 ][i]);
__m128d _mmi1=_mm_set1_pd(-r[1 ][i]);
__m128d _mmi2=_mm_set1_pd(-r[2 ][i]);
int j=0 ;
//*
for(;j+4<i;j+=4 )
{
__m128d dm0=_mm_add_pd(*(__m128d*)&r[0 ][j],_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_add_pd(*(__m128d*)&r[1 ][j],_mmi1);
dm1= _mm_mul_pd(dm1,dm1);
__m128d dm2=_mm_add_pd(*(__m128d*)&r[2 ][j],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm0= _mm_add_pd(dm0,dm1);
dm0= _mm_add_pd(dm0,dm2);
__m128d dm5=_mm_add_pd(*(__m128d*)&r[0][j+2 ],_mmi0);
dm5= _mm_mul_pd(dm5,dm5);
__m128d dm6=_mm_add_pd(*(__m128d*)&r[1][j+2 ],_mmi1);
dm6= _mm_mul_pd(dm6,dm6);
dm2=_mm_add_pd(*(__m128d*)&r[2][j+2 ],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm5= _mm_add_pd(dm5,dm6);
dm5= _mm_add_pd(dm5,dm2);
// 利用IEE double 浮点格式的编码生成1/sqrt(a)的一个近似值； 然后使用牛顿迭代来提高精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
{
__m128i xmmi0= xmmi64_mn;
__m128d dma=_mm_mul_pd(dm0,xmmd_0_5); //a*(-0.5)
__m128d dmb= _mm_mul_pd(dm5,xmmd_0_5);
*(__m128i*)&dm0=_mm_srli_epi64(*(__m128i*)&dm0,1 );
*(__m128i*)&dm5=_mm_srli_epi64(*(__m128i*)&dm5,1 );
*(__m128i*)&dm0=_mm_sub_epi64(xmmi0,*(__m128i*)& dm0);
*(__m128i*)&dm5=_mm_sub_epi64(xmmi0,*(__m128i*)& dm5);
//迭代，加倍精度
dm1 = _mm_mul_pd(dma,dm0);
dm2= _mm_mul_pd(dmb,dm5);
dm1= _mm_mul_pd(dm1,dm0);
dm2= _mm_mul_pd(dm2,dm5);
dm1= _mm_add_pd(dm1,xmmd1_5);
dm2= _mm_add_pd(dm2,xmmd1_5);
dm0= _mm_mul_pd(dm0,dm1);
dm5= _mm_mul_pd(dm5,dm2);
dm1= _mm_mul_pd(dma,dm0);
dm2= _mm_mul_pd(dmb,dm5);
dm1= _mm_mul_pd(dm1,dm0);
dm2= _mm_mul_pd(dm2,dm5);
dm1= _mm_add_pd(dm1,xmmd1_5);
dm2= _mm_add_pd(dm2,xmmd1_5);
dm0= _mm_mul_pd(dm0,dm1);
dm5= _mm_mul_pd(dm5,dm2);
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_add_pd(dma,xmmd1_5);
dmb= _mm_add_pd(dmb,xmmd1_5);
dm0= _mm_mul_pd(dm0,dma);
dm5= _mm_mul_pd(dm5,dmb);
//精度还有点不够 但再次迭代的话就失去了速度优势
}
lcpot= _mm_add_pd(lcpot,dm0);
lcpot= _mm_add_pd(lcpot,dm5);
}
for (;j+1<i;++j)
{
__m128d dm0=_mm_set_sd(r[0][j]);
dm0=_mm_add_pd(dm0,_mmi0);
dm0=_mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_set_sd(r[1][j]);
dm1=_mm_add_sd(dm1,_mmi1);
dm1=_mm_mul_sd(dm1,dm1);
__m128d dm2=_mm_set_sd(r[2][j]);
dm2=_mm_add_sd(dm2,_mmi2);
dm2=_mm_mul_sd(dm2,dm2);
dm0=_mm_add_sd(dm0,dm1);
dm0=_mm_add_sd(dm0,dm2);
const __m128i xmmi64_mn= _mm_set1_epi64(_mmi_mn);
void computePotPart_forj_int(int i,__m128d* pResult)
{
__m128d lcpot= _mm_setzero_pd();
__m128d _mmi0=_mm_set1_pd(-r[0 ][i]);
__m128d _mmi1=_mm_set1_pd(-r[1 ][i]);
__m128d _mmi2=_mm_set1_pd(-r[2 ][i]);
int j=0 ;
//*
for(;j+4<i;j+=4 )
{
__m128d dm0=_mm_add_pd(*(__m128d*)&r[0 ][j],_mmi0);
dm0= _mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_add_pd(*(__m128d*)&r[1 ][j],_mmi1);
dm1= _mm_mul_pd(dm1,dm1);
__m128d dm2=_mm_add_pd(*(__m128d*)&r[2 ][j],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm0= _mm_add_pd(dm0,dm1);
dm0= _mm_add_pd(dm0,dm2);
__m128d dm5=_mm_add_pd(*(__m128d*)&r[0][j+2 ],_mmi0);
dm5= _mm_mul_pd(dm5,dm5);
__m128d dm6=_mm_add_pd(*(__m128d*)&r[1][j+2 ],_mmi1);
dm6= _mm_mul_pd(dm6,dm6);
dm2=_mm_add_pd(*(__m128d*)&r[2][j+2 ],_mmi2);
dm2= _mm_mul_pd(dm2,dm2);
dm5= _mm_add_pd(dm5,dm6);
dm5= _mm_add_pd(dm5,dm2);
// 利用IEE double 浮点格式的编码生成1/sqrt(a)的一个近似值； 然后使用牛顿迭代来提高精度
// 1/sqrt(a)的牛顿迭代公式x_next=(3-a*x*x)*x*0.5 =( 1.5 + (a*(-0.5)) * x*x) ) * x
{
__m128i xmmi0= xmmi64_mn;
__m128d dma=_mm_mul_pd(dm0,xmmd_0_5); //a*(-0.5)
__m128d dmb= _mm_mul_pd(dm5,xmmd_0_5);
*(__m128i*)&dm0=_mm_srli_epi64(*(__m128i*)&dm0,1 );
*(__m128i*)&dm5=_mm_srli_epi64(*(__m128i*)&dm5,1 );
*(__m128i*)&dm0=_mm_sub_epi64(xmmi0,*(__m128i*)& dm0);
*(__m128i*)&dm5=_mm_sub_epi64(xmmi0,*(__m128i*)& dm5);
//迭代，加倍精度
dm1 = _mm_mul_pd(dma,dm0);
dm2= _mm_mul_pd(dmb,dm5);
dm1= _mm_mul_pd(dm1,dm0);
dm2= _mm_mul_pd(dm2,dm5);
dm1= _mm_add_pd(dm1,xmmd1_5);
dm2= _mm_add_pd(dm2,xmmd1_5);
dm0= _mm_mul_pd(dm0,dm1);
dm5= _mm_mul_pd(dm5,dm2);
dm1= _mm_mul_pd(dma,dm0);
dm2= _mm_mul_pd(dmb,dm5);
dm1= _mm_mul_pd(dm1,dm0);
dm2= _mm_mul_pd(dm2,dm5);
dm1= _mm_add_pd(dm1,xmmd1_5);
dm2= _mm_add_pd(dm2,xmmd1_5);
dm0= _mm_mul_pd(dm0,dm1);
dm5= _mm_mul_pd(dm5,dm2);
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_mul_pd(dma,dm0);
dmb= _mm_mul_pd(dmb,dm5);
dma= _mm_add_pd(dma,xmmd1_5);
dmb= _mm_add_pd(dmb,xmmd1_5);
dm0= _mm_mul_pd(dm0,dma);
dm5= _mm_mul_pd(dm5,dmb);
//精度还有点不够 但再次迭代的话就失去了速度优势
}
lcpot= _mm_add_pd(lcpot,dm0);
lcpot= _mm_add_pd(lcpot,dm5);
}
for (;j+1<i;++j)
{
__m128d dm0=_mm_set_sd(r[0][j]);
dm0=_mm_add_pd(dm0,_mmi0);
dm0=_mm_mul_pd(dm0,dm0);
__m128d dm1=_mm_set_sd(r[1][j]);
dm1=_mm_add_sd(dm1,_mmi1);
dm1=_mm_mul_sd(dm1,dm1);
__m128d dm2=_mm_set_sd(r[2][j]);
dm2=_mm_add_sd(dm2,_mmi2);
dm2=_mm_mul_sd(dm2,dm2);
dm0=_mm_add_sd(dm0,dm1);
dm0=_mm_add_sd(dm0,dm2);
{
__m128 sm0=_mm_cvtpd_ps(dm0);
__m128d dma=_mm_mul_sd(dm0,xmmd_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ss(sm0); //计算1/sqrt(a)
dm0=_mm_cvtps_pd(sm0); //
//牛顿迭代，提高开方精度
dm1=_mm_mul_sd(dm0,dm0);
dm1=_mm_mul_sd(dm1,dma);
dm1=_mm_add_sd(dm1,xmmd1_5);
dm0=_mm_mul_sd(dm0,dm1);
__m128 sm0=_mm_cvtpd_ps(dm0);
__m128d dma=_mm_mul_sd(dm0,xmmd_0_5); //a*(-0.5)
sm0=_mm_rsqrt_ss(sm0); //计算1/sqrt(a)
dm0=_mm_cvtps_pd(sm0); //
//牛顿迭代，提高开方精度
dm1=_mm_mul_sd(dm0,dm0);
dm1=_mm_mul_sd(dm1,dma);
dm1=_mm_add_sd(dm1,xmmd1_5);
dm0=_mm_mul_sd(dm0,dm1);
//再次迭代
dma=_mm_mul_sd(dma,dm0);
dma=_mm_mul_sd(dma,dm0);
dma=_mm_add_sd(dma,xmmd1_5);
dm0=_mm_mul_sd(dm0,dma);
}
lcpot= _mm_add_sd(lcpot,dm0);
}
*pResult=_mm_add_pd(* pResult,lcpot);
}
dma=_mm_mul_sd(dma,dm0);
dma=_mm_mul_sd(dma,dm0);
dma=_mm_add_sd(dma,xmmd1_5);
dm0=_mm_mul_sd(dm0,dma);
}
lcpot= _mm_add_sd(lcpot,dm0);
}
*pResult=_mm_add_pd(* pResult,lcpot);
}
还想到一个没有时间去实现的方案：由于SSE和x87是独立的两个硬件，那么可以使它们并行执行；
SSE部件代码不变，但循环做更大的展开，然后让x87承担一路或两路运算;
完
扫一扫
3764
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
