|
|
if (i%10 == 0) printf("%5d: Potential: %20.7f\n", i, pot);
|
distx = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]);
disty = (r[1][j] - r[1][i])*(r[1][j] - r[1][i]);
distz = (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);
#pragma unroll(n)
for(i=0;i<10000;i++){……}
#include <stdio.h>
int main()
{
__int64 start, stop, elaspe;
int i;
int arr[10000];
start= _rdtsc();
for(i=0; i<10000; i++)
{
arr[i]=i;
}
stop= _rdtsc();
elaspe = stop -start;
printf("Processor cycles\n %I64u\n", elaspe);
return 0;
}
for( i=2; i<NPARTS; i++ ) {
for( j=0; j<i-1; j++ ) {
dist = sqrt( (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]) );
pot += 1.0 / dist;
}
}
for( i=2; i<NPARTS; i++ ) {
for( j=0; j<i-1; j++ ) {
dist[j] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);
}
for( j=0; j<i-1; j++ ) {
dist[j] = 1.0 / sqrt(dist[j]);
}
for( j=0; j<i-1; j++ ) {
pot += dist[j];
}
}
for( i=2; i<NPARTS; i++ ) {
for( j=0; j<i-1; j++ ) {
dist[j] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);
}
vdInvSqrt(i-1,dist,dist);
for( j=0; j<i-1; j++ ) {
pot += dist[j];
}
}
There are two extreme cases: so-called "short" and "long" vectors (logarithmic scale is used to show both cases). For short vectors there are cycle organization and initialization overheads. The cost of such overheads is amortized with increasing vector length, and for vectors longer than a few dozens of elements the performance remains quite flat until the L2 cache size is exceeded with the length of the vector.
for( i=2,k=0; i<47; i++ ) {
for( j=0; j<i-1; j++,k++ ) {
dist[k] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);
}
}
vdInvSqrt(k,dist,dist);
for( j=0; j<k; j++ ) {
pot += dist[j];
}
for( i=47; i<NPARTS; i++ ) {
for( j=0; j<i-1; j++ ) {
dist[j] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);
}
vdInvSqrt(i-1,dist,dist);
for( j=0; j<i-1; j++ ) {
pot += dist[j];
}
}
void initPositions() {
int i, j;
for( i=0; i<DIMS; i++ )
for( j=0; j<NPARTS; j++ )
r[i][j] = (double) rand();
}
void updatePositions() {
int i, j;
for( i=0; i<DIMS; i++ )
for( j=0; j<NPARTS; j++ )
r[i][j] -= (double) rand();
}
在main函數(shù)中：
pot = 0.0;
computePot();
pot*=(double)RAND_MAX;
if (i%10 == 0) printf("%5d: Potential: %20.7f\n", i, pot);
int holdrand=1;
inline int myrand (){
return( ((holdrand = holdrand * 214013L+ 2531011L) >> 16) & 0x7fff );
}
__declspec(align(64)) int holdrand=1;
__declspec(align(64)) double r[DIMS][NPARTS];
__declspec(align(64)) double pot;
__declspec(align(64)) double dist[1048];
for( i=2,k=0; i<47; i++ ) {
for( j=0; j<i-1; j++,k++ ) {
dist[k] = (r[0][j] - r[0][i])*(r[0][j] - r[0][i]) + (r[1][j] - r[1][i])*(r[1][j] - r[1][i]) + (r[2][j] - r[2][i])*(r[2][j] - r[2][i]);
}
}
vdInvSqrt(k,dist,dist);
for( j=0; j<1035; j++ ) {
pot += dist[j];
}
|
Full Name
|
CPI
|
Clockticks events
|
Clockticks %
|
|
void updatePositions(void)
|
3.214080375
|
8274621
|
0.287907869
|
|
int computePot(void)
|
1.294881302
|
926757552
|
32.24568138
|
|
mkl_vml_core_t7_vml_dInvSqrt_50
|
0.91981472
|
1925228486
|
66.9865643
|
int computePot1(void);
void initPositions1(void);
void updatePositions1(void);
int computePot2(void);
void initPositions2(void);
void updatePositions2(void);
__declspec(align(64)) int holdrand1=1;
__declspec(align(64)) double r1[DIMS][NPARTS];
__declspec(align(64)) double pot1;
__declspec(align(64)) double dist1[1048];
__declspec(align(64)) int holdrand2=1;
__declspec(align(64)) double r2[DIMS][NPARTS];
__declspec(align(64)) double pot2;
__declspec(align(64)) double dist2[1048];
__declspec(align(64)) double potfinal[264];
DWORD WINAPI mythread( void *myarg ){
int i;
SetThreadAffinityMask(GetCurrentThread(), 2);
initPositions2();
updatePositions2();
for(i=0;i<=190;i+=2){
pot2 = 0.0;
computePot2();
pot2*=(double)RAND_MAX;
potfinal[i]=pot2;
updatePositions2();
updatePositions2();
}
ReleaseSemaphore(semmiddle, 1, NULL);
for(i=192;i<=NITER;i+=2){
pot2 = 0.0;
computePot2();
pot2*=(double)RAND_MAX;
potfinal[i]=pot2;
updatePositions2();
updatePositions2();
}
ReleaseSemaphore(semafter, 1, NULL);
return 0;
}//從線程
int main() {
int i;
int myarg=0;
clock_t start, stop;
omp_set_num_threads(1);
vmlSetMode(VML_LA);
vmlSetMode(VML_DOUBLE_CONSISTENT);
vmlSetMode(VML_ERRMODE_IGNORE);
vmlSetMode(VML_NUM_THREADS_OMP_FIXED);
semmiddle = CreateSemaphore(NULL, 0, 1, NULL);
semafter = CreateSemaphore(NULL, 0, 1, NULL);
CreateThread(0, 8*1024, mythread, (void *)&myarg, 0, NULL);
SetThreadAffinityMask(GetCurrentThread(), 1);
initPositions1();
start=clock();
for(i=1;i<NITER;i+=2){
pot1 = 0.0;
updatePositions1();
updatePositions1();
computePot1();
pot1*=(double)RAND_MAX;
potfinal[i]=pot1;
}
WaitForSingleObject(semmiddle, INFINITE);
for(i=0;i<=190;i+=10)
printf("%5d: Potential: %20.7f\n", i, potfinal[i]);
WaitForSingleObject(semafter , INFINITE);
i=200;
printf("%5d: Potential: %20.7f\n", i, potfinal[i]);
stop=clock();
printf ("Seconds = %10.9f\n",(double)(stop-start)/ CLOCKS_PER_SEC);
}//主線程
float InvSqrt(float x){
float xhalf = 0.5f*x;
int i = *(int*)&x; // get bits for floating value
i = 0x5f3759df - (i>>1); // gives initial guess y0
x = *(float*)&i; // convert bits back to float
x = x*(1.5f-xhalf*x*x); // Newton step, repeating increases accuracy
return x;
}
double myinvsqrt (double x)
{
double xhalf = 0.5*x;
__int64 i = *(__int64*)&x;
i = 0x5fe6ec85e7de30da - (i>>1);
x = *(double*)&i;
x = x*(1.5-xhalf*x*x);
x = x*(1.5-xhalf*x*x);
x = x*(1.5-xhalf*x*x);
x = x*(1.5-xhalf*x*x);
return x;
}
double myinvsqrt (double x)
{
double xhalf = 0.5*x;
__int64 i = *(__int64*)&x;
i = newmagicnum - (i>>1);
x = *(double*)&i;
x = x*(1.5-xhalf*x*x);
x = x*(1.5-xhalf*x*x);
x = x*(1.5-xhalf*x*x);
x = x*offset
return x;
}
double myinvsqrt (double x)
{
double xhalf = 0.5*x;
float xf=(float)x;
__asm{
movss xmm1,xf;
rsqrtss xmm1,xmm1;
movss xf,xmm1;
}
x=(double)xf;
x = x*(1.5-xhalf*x*x);
x = x*(1.5-xhalf*x*x);
return x;
}
double myinvsqrt (double x)
{
double xhalf = -0.5*x;
float xf=(float)x;
__asm{
movss xmm1,xf;
rsqrtss xmm1,xmm1;
movss xf,xmm1;
}
x=(double)xf;
x = x*(1.5+xhalf*x*x);
x = x*(1.5+xhalf*x*x);
return x;
}
|
指令名
|
延遲
|
吞吐率
|
計(jì)算部件
|
|
movapd
|
1
|
0.33
|
FP_MOVE
|
|
cvtpd2ps
|
4
|
1
|
FP_ADD,MMX_SHFT
|
|
cvtps2pd
|
2
|
2
|
FP_ADD,MMX_SHFT,MMX_ALU
|
|
shufps
|
2
|
1
|
MMX_SHFT
|
|
rsqrtps
|
3
|
2
|
MMX_MISC
|
|
mulpd
|
5
|
1
|
FP_MUL
|
|
addpd
|
3
|
1
|
FP_ADD
|
void myinvsqrt (double *start,double *end)
{
__asm{
mov esi,start;
mov edi,end;
test edi,0x0000001f;
jz myalign;
and edi,0xffffffe0;
add edi,32;
myalign:
myagain:
movapd xmm0,[esi];
movapd xmm3,[esi+16];
cvtpd2ps xmm6,xmm0;
cvtpd2ps xmm7,xmm3;
shufps xmm6,xmm7,01000100b;
rsqrtps xmm6,xmm6;
cvtps2pd xmm1,xmm6;
shufps xmm6,xmm6,01001110b;
cvtps2pd xmm4,xmm6;
mulpd xmm0,mulcc;
mulpd xmm3,mulcc;
movapd xmm2,xmm1;
movapd xmm5,xmm4;
mulpd xmm1,xmm1;
mulpd xmm4,xmm4;
mulpd xmm1,xmm0;
mulpd xmm4,xmm3;
addpd xmm1,addcc;
addpd xmm4,addcc;
mulpd xmm1,xmm2;
mulpd xmm4,xmm5;//前半段
movapd xmm2,xmm1;
movapd xmm5,xmm4;
mulpd xmm1,xmm1;
mulpd xmm4,xmm4;
mulpd xmm1,xmm0;
mulpd xmm4,xmm3;
addpd xmm1,addcc;
addpd xmm4,addcc;
mulpd xmm1,xmm2;
mulpd xmm4,xmm5;
movapd [esi],xmm1;
movapd [esi+16],xmm4;
add esi,32;
cmp esi,edi;
jne myagain;
}
}
//后半段
myinvsqrt(dist1,dist1+k); //調(diào)用方法
發(fā)表于 @ 2007年01月21日 13:58:00|評(píng)論(2 )|編輯
|
|
英特爾多核平臺(tái)編程優(yōu)化大賽報(bào)告
