Repository navigation
Loop unrolling
salehqt edited this page Sep 4, 2012
·
4 revisions
We examined the file src/swarm/gpu/gravitation_accjerk.hpp and compared the C++ source code and PTX code. The loops are actually getting unrolled. Here are the extracts of the source code and PTX output for reference:
This function is called when calcForces is invoked:
GPUAPI void operator() (int ij,int b,int c,double& pos,double& vel,double& acc,double& jerk) const{
// Write positions to shared (global) memory
if(b < nbod && c < 3)
sys[b][c].pos() = pos , sys[b][c].vel() = vel;
__syncthreads();
if(ij < pair_count)
calc_pair(ij);
__syncthreads();
if(b < nbod && c < 3){
sum(b,c,acc,jerk);
}
}
First few lines are optimized, the next major thing to do is this line:
if(ij < pair_count)
Which is translated to: (PTX Code)
BB0_38:
#line gravitation_accjerk.hpp:101 1
cvt.s64.s32 %rl88, %r4;
mad.lo.s64 %rl32, %rl88, 768, %rl8;
.loc 2 84 1
setp.gt.s32 %p52, %r4, 2;
@%p52 bra BB0_41;
The last line jumps over the statements in the if (bra)
Then this function is called:
GENERIC void calc_pair(int ij)const{
int i = first<nbod>( ij );
int j = second<nbod>( ij );
if(i != j){
// Relative vector from planet i to planet j
double dx[3] = { sys[j][0].pos()- sys[i][0].pos(),sys[j][1].pos()- sys[i][1].pos(), sys[j][2].pos()- sys[i][2].pos() };
// Relative velocity between the planets
double dv[3] = { sys[j][0].vel()- sys[i][0].vel(),sys[j][1].vel()- sys[i][1].vel(), sys[j][2].vel()- sys[i][2].vel() };
// Distance between the planets
double r2 = dx[0]*dx[0]+dx[1]*dx[1]+dx[2]*dx[2];
double jerk_mag = inner_product(dx,dv) * 3.0 / r2;
double acc_mag = rsqrt(r2) / r2;
// this loop actually unrolls in PTX code.
#pragma unroll
for(int c = 0; c < 3; c ++)
{
shared[ij][c].acc() = dx[c]* acc_mag;
shared[ij][c].jerk() = (dv[c] - dx[c] * jerk_mag ) * acc_mag;
}
}
}
Calculation of i and j and the overall if statement is tranlated to:
mov.u32 %r203, 2;
.loc 11 82 1
sub.s32 %r204, %r203, %r4;
.loc 11 87 1
sub.s32 %r205, %r203, %r204;
.loc 11 84 1
setp.gt.s32 %p53, %r204, 0;
selp.b32 %r20, %r204, %r205, %p53;
setp.lt.s32 %p11, %r204, 1;
selp.u32 %r206, 1, 0, %p11;
#line gravitation_accjerk.hpp:85 1
setp.eq.s32 %p54, %r20, %r206;
@%p54 bra BB0_41;
Now in the body of the if statement, we have some calculations to do:
// Relative vector from planet i to planet j
double dx[3] = { sys[j][0].pos()- sys[i][0].pos(),sys[j][1].pos()- sys[i][1].pos(), sys[j][2].pos()- sys[i][2].pos() };
// Relative velocity between the planets
double dv[3] = { sys[j][0].vel()- sys[i][0].vel(),sys[j][1].vel()- sys[i][1].vel(), sys[j][2].vel()- sys[i][2].vel() };
// Distance between the planets
double r2 = dx[0]*dx[0]+dx[1]*dx[1]+dx[2]*dx[2];
double jerk_mag = inner_product(dx,dv) * 3.0 / r2;
double acc_mag = rsqrt(r2) / r2;
These lines are tranlated to:
selp.u64 %rl89, 1, 0, %p11;
shl.b64 %rl90, %rl89, 10;
add.s64 %rl91, %rl11, %rl90;
mul.wide.s32 %rl92, %r20, 1024;
add.s64 %rl93, %rl11, %rl92;
#line gravitation_accjerk.hpp:88 1
ld.f64 %fd150, [%rl93];
ld.f64 %fd151, [%rl91];
sub.f64 %fd152, %fd151, %fd150;
ld.f64 %fd153, [%rl93+256];
ld.f64 %fd154, [%rl91+256];
sub.f64 %fd155, %fd154, %fd153;
ld.f64 %fd156, [%rl93+512];
ld.f64 %fd157, [%rl91+512];
sub.f64 %fd158, %fd157, %fd156;
#line gravitation_accjerk.hpp:90 1
ld.f64 %fd159, [%rl93+128];
ld.f64 %fd160, [%rl91+128];
sub.f64 %fd161, %fd160, %fd159;
ld.f64 %fd162, [%rl93+384];
ld.f64 %fd163, [%rl91+384];
sub.f64 %fd164, %fd163, %fd162;
ld.f64 %fd165, [%rl93+640];
ld.f64 %fd166, [%rl91+640];
sub.f64 %fd167, %fd166, %fd165;
#line gravitation_accjerk.hpp:93 1
mul.f64 %fd168, %fd155, %fd155;
fma.rn.f64 %fd169, %fd152, %fd152, %fd168;
fma.rn.f64 %fd170, %fd158, %fd158, %fd169;
#line gravitation_accjerk.hpp:95 1
mul.f64 %fd171, %fd155, %fd164;
fma.rn.f64 %fd172, %fd152, %fd161, %fd171;
fma.rn.f64 %fd173, %fd158, %fd167, %fd172;
mul.f64 %fd174, %fd173, 0d4008000000000000;
.loc 7 1328 3
div.rn.f64 %fd175, %fd174, %fd170;
.loc 6 470 5
rsqrt.approx.f64 %fd176, %fd170;
.loc 7 1328 3
div.rn.f64 %fd177, %fd176, %fd170;
Once two values acc_mag and jerk_mag are ready (%fd177, %fd175) we can proceed with unrolling the loop that comes on:
#pragma unroll
for(int c = 0; c < 3; c ++)
{
shared[ij][c].acc() = dx[c]* acc_mag;
shared[ij][c].jerk() = (dv[c] - dx[c] * jerk_mag ) * acc_mag;
}
You may recognize that the operations are repeated 3 times in the translation of the above for loop and there is no branching involved:
#line gravitation_accjerk.hpp:101 1
mul.f64 %fd178, %fd152, %fd177;
st.shared.f64 [%rl32], %fd178;
#line gravitation_accjerk.hpp:102 1
neg.f64 %fd179, %fd152;
fma.rn.f64 %fd180, %fd179, %fd175, %fd161;
mul.f64 %fd181, %fd180, %fd177;
st.shared.f64 [%rl32+128], %fd181;
#line gravitation_accjerk.hpp:101 1
mul.f64 %fd182, %fd155, %fd177;
st.shared.f64 [%rl32+256], %fd182;
#line gravitation_accjerk.hpp:102 1
neg.f64 %fd183, %fd155;
fma.rn.f64 %fd184, %fd183, %fd175, %fd164;
mul.f64 %fd185, %fd184, %fd177;
st.shared.f64 [%rl32+384], %fd185;
#line gravitation_accjerk.hpp:101 1
mul.f64 %fd186, %fd158, %fd177;
st.shared.f64 [%rl32+512], %fd186;
#line gravitation_accjerk.hpp:102 1
neg.f64 %fd187, %fd158;
fma.rn.f64 %fd188, %fd187, %fd175, %fd167;
mul.f64 %fd189, %fd188, %fd177;
st.shared.f64 [%rl32+640], %fd189;
Here we exit calc_pair function (there is no return statement since all these functions are inlined)
Going back to operator(), we have these lines:
__syncthreads();
if(b < nbod && c < 3){
that easily translate to:
BB0_41:
.loc 2 84 1
bar.sync 0;
shl.b64 %rl94, %rl10, 8;
add.s64 %rl33, %rl8, %rl94;
shl.b64 %rl95, %rl3, 10;
add.s64 %rl96, %rl5, %rl95;
shl.b64 %rl97, %rl4, 3;
add.s64 %rl34, %rl97, %rl96;
.loc 2 84 1
@%p1 bra BB0_43;
mov.f64 %fd460, 0d0000000000000000;
mov.f64 %fd450, %fd460;
bra.uni BB0_57;
the first bra jumps over the else clause, the else clause is the last 3 lines setting acc and jerk to something (probably NAN).
going into the body of the if statement we enter the sum function:
GENERIC void sum(int b,int c,double& acc, double & jerk)const{
// Total acceleration from other planets
double acc_from_planets = 0.0;
// Total acceleration from body 0 (sun or star)
double acc_from_sun = 0.0;
// Total jerk from other planets
double jerk_from_planets = 0.0;
// Total jerk from body 0 (sun or star)
double jerk_from_sun = 0.0;
/// Find the contribution from/to Sun first
// this loop is actually unrolled for 3 bodies
#pragma unroll
for(int d = 0; d < pair_count; d++){
int x = first<nbod>(d), y= second<nbod>(d);
if(x == b){
if(y == 0) {
acc_from_sun += shared[d][c].acc() * sys[y].mass();
jerk_from_sun += shared[d][c].jerk() * sys[y].mass();
} else {
acc_from_planets += shared[d][c].acc() * sys[y].mass();
jerk_from_planets += shared[d][c].jerk() * sys[y].mass();
}
}else if(y == b){
if(x == 0) {
acc_from_sun -= shared[d][c].acc() * sys[x].mass();
jerk_from_sun -= shared[d][c].jerk() * sys[x].mass();
} else {
acc_from_planets -= shared[d][c].acc() * sys[x].mass();
jerk_from_planets -= shared[d][c].jerk() * sys[x].mass();
}
}
}
acc = acc_from_sun + acc_from_planets;
jerk = jerk_from_sun + jerk_from_planets;
}
The first few lines are just declarations, we enter the loop right away. but as you can see there are no loop variables or branching backwards.
The unrolling is very difficult to verify since all the if statements inside the loop are also transformed and optimized, too. From here to the end of the loop there are several paths that can be taken based on value of b.
BB0_43:
#line gravitation_accjerk.hpp:133 1
setp.eq.s32 %p12, %r5, 2;
@%p12 bra BB0_47;
#line gravitation_accjerk.hpp:141 1
setp.eq.s32 %p55, %r5, 0;
@%p55 bra BB0_46;
mov.f64 %fd432, 0d0000000000000000;
mov.f64 %fd434, %fd432;
mov.f64 %fd431, %fd432;
mov.f64 %fd433, %fd432;
bra.uni BB0_48;
BB0_46:
#line gravitation_accjerk.hpp:143 1
ld.shared.f64 %fd198, [%rl33];
ld.f64 %fd199, [%rl34+2816];
mov.f64 %fd432, 0d0000000000000000;
#line gravitation_accjerk.hpp:146 1
neg.f64 %fd200, %fd198;
fma.rn.f64 %fd433, %fd200, %fd199, %fd432;
#line gravitation_accjerk.hpp:147 1
ld.shared.f64 %fd201, [%rl33+128];
neg.f64 %fd202, %fd201;
fma.rn.f64 %fd434, %fd202, %fd199, %fd432;
mov.f64 %fd431, %fd432;
bra.uni BB0_48;
BB0_47:
#line gravitation_accjerk.hpp:135 1
ld.shared.f64 %fd205, [%rl33];
ld.f64 %fd206, [%rl34+768];
fma.rn.f64 %fd431, %fd205, %fd206, 0d0000000000000000;
mov.f64 %fd434, 0d0000000000000000;
#line gravitation_accjerk.hpp:136 1
ld.shared.f64 %fd207, [%rl33+128];
fma.rn.f64 %fd432, %fd207, %fd206, 0d0000000000000000;
mov.f64 %fd433, %fd434;
BB0_48:
#line gravitation_accjerk.hpp:133 1
setp.eq.s32 %p13, %r5, 1;
@%p13 bra BB0_51;
#line gravitation_accjerk.hpp:141 1
setp.eq.s32 %p56, %r5, 0;
@%p56 bra BB0_50;
bra.uni BB0_52;
BB0_50:
#line gravitation_accjerk.hpp:143 1
ld.f64 %fd208, [%rl34+1792];
ld.shared.f64 %fd209, [%rl33+768];
#line gravitation_accjerk.hpp:146 1
neg.f64 %fd210, %fd209;
fma.rn.f64 %fd433, %fd210, %fd208, %fd433;
#line gravitation_accjerk.hpp:147 1
ld.shared.f64 %fd211, [%rl33+896];
neg.f64 %fd212, %fd211;
fma.rn.f64 %fd434, %fd212, %fd208, %fd434;
bra.uni BB0_52;
BB0_51:
#line gravitation_accjerk.hpp:135 1
ld.f64 %fd213, [%rl34+768];
ld.shared.f64 %fd214, [%rl33+768];
fma.rn.f64 %fd431, %fd214, %fd213, %fd431;
#line gravitation_accjerk.hpp:136 1
ld.shared.f64 %fd215, [%rl33+896];
fma.rn.f64 %fd432, %fd215, %fd213, %fd432;
BB0_52:
#line gravitation_accjerk.hpp:133 1
@%p12 bra BB0_55;
@%p13 bra BB0_54;
bra.uni BB0_56;
BB0_54:
#line gravitation_accjerk.hpp:143 1
ld.f64 %fd216, [%rl34+2816];
ld.shared.f64 %fd217, [%rl33+1536];
#line gravitation_accjerk.hpp:146 1
neg.f64 %fd218, %fd217;
fma.rn.f64 %fd433, %fd218, %fd216, %fd433;
#line gravitation_accjerk.hpp:147 1
ld.shared.f64 %fd219, [%rl33+1664];
neg.f64 %fd220, %fd219;
fma.rn.f64 %fd434, %fd220, %fd216, %fd434;
bra.uni BB0_56;
BB0_55:
#line gravitation_accjerk.hpp:135 1
ld.f64 %fd221, [%rl34+1792];
ld.shared.f64 %fd222, [%rl33+1536];
#line gravitation_accjerk.hpp:138 1
fma.rn.f64 %fd433, %fd222, %fd221, %fd433;
#line gravitation_accjerk.hpp:139 1
ld.shared.f64 %fd223, [%rl33+1664];
fma.rn.f64 %fd434, %fd223, %fd221, %fd434;
Finally, we have to add the contribution from sun and other planets to get the total (4 registers):
BB0_56:
#line gravitation_accjerk.hpp:152 1
add.f64 %fd450, %fd431, %fd433;
#line gravitation_accjerk.hpp:153 1
add.f64 %fd460, %fd432, %fd434;
Now %fd450 and %fd460 hold the value for acceleration and jerk (I don't know which is which).