Publish the source code

            


What I did not blog or online space. Here only released the source code.
Some explanation:
1, I have to write code that may not be very nice.
2, the code for the Intel Core 2 32KB L1 cache tuning, performance on other platforms may not be perfect.
3, the code on the Core 2 test, reaching the required precision. But it has not been tested on other CPU.

Reply:
/ * Compute the potential energy of a collection of * /
/ * Particles interacting via pairwise potential * /

#include & lt; stdio.h & gt;
#include & lt; stdlib.h & gt;
#include & lt; windows.h & gt;
#include & lt; time.h & gt;


#define NPARTS 1000
#define NITER 201
#define DIMS 3


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 mulcc [2] = {- 0.5, -0.5};
__declspec (align (64)) double addcc [2] = {1.5,1.5};

__declspec (align (64)) double potfinal [264];

HANDLE semmiddle;
HANDLE semafter;

Reply:

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;
}
}

inline int myrand1 () {
return (((holdrand1 = holdrand1 * 214013L + 2531011L) & gt; & gt; 16) & amp; 0x7fff);
}

inline int myrand2 () {
return (((holdrand2 = holdrand2 * 214013L + 2531011L) & gt; & gt; 16) & amp; 0x7fff);
}

Reply:

DWORD WINAPI mythread (void * myarg) {
int i;
SetThreadAffinityMask (GetCurrentThread (), 2);
dist2 [1035] = 1;
initPositions2 ();
updatePositions2 ();
for (i = 0; i & lt; = 190; i + = 2) {
pot2 = 0.0;
computePot2 ();
pot2 * = (double) RAND_MAX;
potfinal [i] = pot2;
updatePositions2 ();
updatePositions2 ();
}
ReleaseSemaphore (semmiddle, 1, NULL);
for (i = 192; i & lt; = 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;
semmiddle = CreateSemaphore (NULL, 0, 1, NULL);
semafter = CreateSemaphore (NULL, 0, 1, NULL);
CreateThread (0, 8 * 1024, mythread, (void *) & amp; myarg, 0, NULL);
SetThreadAffinityMask (GetCurrentThread (), 1);
initPositions1 ();
start = clock ();
dist1 [1035] = 1;
for (i = 1; i & lt; NITER; i + = 2) {
pot1 = 0.0;
updatePositions1 ();
updatePositions1 ();
computePot1 ();
pot1 * = (double) RAND_MAX;
potfinal [i] = pot1;
}
WaitForSingleObject (semmiddle, INFINITE);
for (i = 0; i & lt; = 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);

}


Reply:
The landlord, you too, you first look at the top of the assembly that N long learning ~~ ~~~
Reply:
The top one first
Reply:
I was finally able to continue replied, then paste.

void initPositions1 () {
int i, j;
for (i = 0; i & lt; DIMS; i ++)
for (j = 0; j & lt; NPARTS; j ++)
r1 [i] [j] = (double) myrand1 ();
}

void initPositions2 () {
int i, j;
for (i = 0; i & lt; DIMS; i ++)
for (j = 0; j & lt; NPARTS; j ++)
r2 [i] [j] = (double) myrand2 ();
}

void updatePositions1 () {
int i, j;
for (i = 0; i & lt; DIMS; i ++)
for (j = 0; j & lt; NPARTS; j ++)
r1 [i] [j] - = (double) myrand1 ();
}

void updatePositions2 () {
int i, j;
for (i = 0; i & lt; DIMS; i ++)
for (j = 0; j & lt; NPARTS; j ++)
r2 [i] [j] - = (double) myrand2 ();
}

Reply:

int computePot1 () {
int i, j, k;
for (i = 2, k = 0; i & lt; 47; i ++) {
for (j = 0; j & lt; i-1; j ++, k ++) {
dist1 [k] = (r1 [0] [j] - r1 [0] [i]) * (r1 [0] [j] - r1 [0] [i]) + (r1 [1] [j] - r1 [1] [i]) * (r1 [1] [j] - r1 [1] [i]) + (r1 [2] [j] - r1 [2] [i]) * (r1 [2] [j] - r1 [2] [i]);
}
}
myinvsqrt (dist1, dist1 + k);
for (j = 0; j & lt; 1035; j ++) {
pot1 + = dist1 [j];
}
for (i = 47; i & lt; NPARTS; i ++) {
for (j = 0; j & lt; i-1; j ++) {
dist1 [j] = (r1 [0] [j] - r1 [0] [i]) * (r1 [0] [j] - r1 [0] [i]) + (r1 [1] [j] - r1 [1] [i]) * (r1 [1] [j] - r1 [1] [i]) + (r1 [2] [j] - r1 [2] [i]) * (r1 [2] [j] - r1 [2] [i]);
}
myinvsqrt (dist1, dist1 + i-1);
for (j = 0; j & lt; i-1; j ++) {
pot1 + = dist1 [j];
}
}
return 0;
}

int computePot2 () {
int i, j, k;
for (i = 2, k = 0; i & lt; 47; i ++) {
for (j = 0; j & lt; i-1; j ++, k ++) {
dist2 [k] = (r2 [0] [j] - r2 [0] [i]) * (r2 [0] [j] - r2 [0] [i]) + (r2 [1] [j] - r2 [1] [i]) * (r2 [1] [j] - r2 [1] [i]) + (r2 [2] [j] - r2 [2] [i]) * (r2 [2] [j] - r2 [2] [i]);
}
}
myinvsqrt (dist2, dist2 + k);
for (j = 0; j & lt; 1035; j ++) {
pot2 + = dist2 [j];
}
for (i = 47; i & lt; NPARTS; i ++) {
for (j = 0; j & lt; i-1; j ++) {
dist2 [j] = (r2 [0] [j] - r2 [0] [i]) * (r2 [0] [j] - r2 [0] [i]) + (r2 [1] [j] - r2 [1] [i]) * (r2 [1] [j] - r2 [1] [i]) + (r2 [2] [j] - r2 [2] [i]) * (r2 [2] [j] - r2 [2] [i]);
}
myinvsqrt (dist2, dist2 + i-1);
for (j = 0; j & lt; i-1; j ++) {
pot2 + = dist2 [j];
}
}
return 0;
}


All posted good.

Reply:
Download the complete learning! :)
Reply:
Attach the sample period compiled c code. And assembly is still a difference, but on the same algorithm.
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;
}

Reply:
Learn about, thanks
Reply:
to: flyingdog (flyingdog)
Oh, we are very close to the core code ah; do not know who's speed will be faster;
You try to improve myinvsqrt function, see there will be any increase:
Such a change: After rsqrtps, first conduct a Newton iteration; and then converted to a double iteration
Reply:
to: housisong (HouSisong)
You strong! Even the whole computePot are a compilation.

But I have a problem.
After rsqrtps, the float is already accurate result. Why in the first iteration float?
I double iteration 2 in place, because a lack of precision iteration.
Is it true that the first iteration float, double precision iterative once enough?

Reply:
Are optimized to the assembler level, ah

Humbled, I can only optimized to 0.609 seconds; however, the test platform is Intel's first dual-core platform Pentium D, I do not know will not change in a hurry Core2duo platform
Reply:
to: flyingdog (flyingdog)
rsqrtps calculate approximations only, so you need to iterate to get more accurate results
float a double iteration iteration time, the accuracy is sufficient
9 You can try, I do so on AMD64x2 and p4, the accuracy can be decimal
Reply:
Intel's mkl test it in invsqrt only precision to eight decimal places, the first nine and the original a little deviation!
----------------------------------
housisong
9 You can try, I do so on AMD64x2 and p4, the accuracy can be decimal
----------------------------------
Learn from you!
Reply:
Just then at the document. The original rsqrtps are approximate, I always thought that directly correspond to the CPU instruction set is precisely calculated. This time loss. Fortunately, if the pipeline, then to do a Newton iteration float and double do a Newton iteration almost time.

In fact, it should be said that the accuracy is not accurate to several decimal places, but several significant figures. For example, this program last number, the integer part is relatively small, accurate to seven decimal places is easy. Digital integer part of the first few relatively large. But hard to one decimal seven.

Reply:
To learn from you!

Reply:
defined rsqrtps approximation is calculated, and allowing different CPU computing have different implementations (a little nervous)
rcpps (a reciprocal approximation) is such instruction;

I think this design is ready for the game, the game speed is more important in a lot of time; people who want high-precision only need a simple Newton iteration just once (CPU own internal iteration this time it would not be able to be faster );
This design is more common;
Reply:
I also think this design is more universal, the better.

Reply:
I have just opened a blog, my reports put up optimization.
http://blog.csdn.net/flyingdog/archive/2007/01/21/1489183.aspx
Also please correct me if wrong.

Since I was copied directly from the word go, the less the picture, format is also somewhat problematic.
If you want to complete reporting and optimization program, probably a little more than 1M. Please leave e-mail, I sent you.

Reply:
I read your report; one of the "data parallelism" is my greatest achievement, thank

And it seems that you have done a lot of research data and qualitative performance analysis; and I do not have much tools are available;
(AMD's CPU is not running intel tools, have to rely on their own)
blog support map, you try;
Give me a optimization report, thank you; my mailbox HouSisong@263.net
Reply:
`` Learning
Reply:
to housisong
Has been sent, please check.

Reply:
Strongman !!!!

Reply:
House Party really powerful ah.
Although they are not my people, the poor people are eager to come up with this myself, I have been Cougerenao, my report:
http://blog.csdn.net/aero_boy/archive/2007/01/21/1489357.aspx
Reply:
It seems that my program might be the least of the changes to the program, and how well you learn.
Reply:
The landlord too, thanks to a another landlord seeking Optimization Report
Email: m2213231@yahoo.com.cn
Thank you
Reply:
to m2213231
Has been sent, please check.

Reply:
Mark
Reply:
The landlord is too strong,
Email: wai631@gmail.com
Reply:
The report looked at the landlord, very powerful! I learned a lot ah!
Give me a full version of the report it. Thank you.
email: lanphaday@126.com
Reply:
Receive, thank landlord

Reply:
Learning
I have to give a copy of it, thank you!
shishi76@126.com
Reply:
Learning
luxu01207@163.com
Reply:
Has been sent
Reply:
Learn about, look also gave me a copy, thank you!
sg47@163.com
Reply:
Strongman ah! Assembler and C mixed ah!
Reply:
Sent
luxu01207@163.com Send failed.
Reply:
LZ good and strong, give me a handle, to learn from you
CharlesLuan@hotmail.com
Reply:
Super strong!
Reply:
Has been sent
Reply:
wyvernthinker@gmail.com
Thank you!
Reply:
Has been sent
Reply:
like you ,,,,,, lz .....

mark
Reply:
To see such a good optimization, really very happy ~~~
flyingdog trouble also gave me a copy of it, I read about ~~~
Reply:
playboy66@21cn.com
Reply:
Can be, then I also want to report, thank floor living
liuyanches@163.com
Reply:
mark
Reply:
The landlord wanted to learn, can be, then I also want to report, thank floor living
syf7901@163.com

Reply:
Severe learning