lavaMD

The code calculates particle potential and relocation due to mutual forces between particles within a large 3D space. This space is divided into cubes, or large boxes, that are allocated to individual cluster nodes. The large box at each node is further divided into cubes, called boxes. 26 neighbor boxes surround each box (the home box). Home boxes at the boundaries of the particle space have fewer neighbors. Particles only interact with those other particles that are within a cutoff radius since ones at larger distances exert negligible forces. Thus the box size is chosen so that cutoff radius does not span beyond any neighbor box for any particle in a home box, thus limiting the reference space to a finite number of boxes. (from the Rodinia README)

Code snippet

#deine NUMBER_PAR_PER_BOX 100
#define DOT(A,B) ((A.x)*(B.x)+(A.y)*(B.y)+(A.z)*(B.z))

dim.number_boxes = 512;

for(l=0; l<dim.number_boxes; l=l+1){
   //home box - box parameters
   first_i = box[l].offset;      // offset to common arrays

   //home box - distance, force, charge and type parameters from common arrays
   rA = &rv[first_i];
   fA = &fv[first_i];

   //	Do for the # of (home+neighbor) boxes
   for (k=0; k<(1+box[l].nn); k++){      //executes for 22 iterations on average
      numIter++;
      avgNum += (1+box[l].nn);

      //neighbor box - get pointer to the right box
      if(k==0){
         pointer = l;      // set first box to be processed to home box
      }else{
         pointer = box[l].nei[k-1].number;      // remaining boxes are neighbor boxes
      }

      //neighbor box - box parameters
      first_j = box[pointer].offset; 
      //	neighbor box - distance, force, charge and type parameters
      rB = &rv[first_j];
      qB = &qv[first_j];

      //Do for the # of particles in home box
      for (i=0; i<NUMBER_PAR_PER_BOX; i=i+1){

         //do for the # of particles in current (home or neighbor) box
         for (j=0; j<NUMBER_PAR_PER_BOX; j=j+1){
            //coefficients
            r2 = rA[i].v+rB[j].v - (rA[i].x*rB[j].x + rA[i].y*rB[j].y + rA[i].z*rB[j].z); 
            u2 = a2*r2;
            vij= expf(-u2);
            fs = 2.*vij;
            d.x = rA[i].x  - rB[j].x; 
            d.y = rA[i].y  - rB[j].y; 
            d.z = rA[i].z  - rB[j].z; 
            fxij=fs*d.x;
            fyij=fs*d.y;
            fzij=fs*d.z;
            // forces
            fA[i].v +=  qB[j]*vij;
            fA[i].x +=  qB[j]*fxij;
            fA[i].y +=  qB[j]*fyij;
            fA[i].z +=  qB[j]*fzij;
         }   // for j
      }   // for i
   }   // for k
}   // for l