source: proiecte/NBody/NBody 2.0/Kernels/simpleIntegratorKernel.cl @ 34

/*
 * Each work-item invocation of this kernel, calculates the position for 
 * one particle
 *
 * Work-items use local memory to reduce memory bandwidth and reuse of data
 */

__kernel void simple_integrator(        __global float4* pos ,

                                                                        __global float4* vel,

                                                                        int numBodies,

                                                                    float deltaTime,
                
                                                                    float epsSqr,

                                                                    __local float4* localPos)

{
        //---------------------------------------
    unsigned int tid = get_local_id(0);         

    unsigned int gid = get_global_id(0);                // Index of the particle that it must update;

    unsigned int localSize = get_local_size(0); // This can be thought of as the number of threads
                                                                                                // executing this kernel, for different work-items, 
                                                                                                // concurrently and synchronously;

        //---------------------------------------
    // Number of tiles we need to iterate

    unsigned int numTiles = numBodies / localSize;


        //---------------------------------------
        // Reads the particle position (and mass) and velocity 
        // of particle "i" for which this kernel invocation 
        // is tasked to update. 
        
    float4 oldPos = pos[gid];
    float4 oldVel = vel[gid];

        //---------------------------------------
        // Initializes the float4 we will use to accumulate
        // the acceleration on particle "i".
        
    float4 acc = (float4)(0.0f, 0.0f, 0.0f, 0.0f);


    for(int i = 0; i < numTiles; ++i)

    {

        // load one tile into local memory

        int idx = i * localSize + tid;

        localPos[tid] = pos[idx];



        // Synchronize to make sure data is available for processing

        barrier(CLK_LOCAL_MEM_FENCE);


        // calculate acceleration effect due to each body

        // a[i->j] = m[j] * r[i->j] / (r^2 + epsSqr)^(3/2)

        for(int j = 0; j < localSize; ++j)

        {

            // Calculate acceleartion caused by particle j on particle i

            float4 r = localPos[j] - oldPos;

            float distSqr = r.x * r.x  +  r.y * r.y  +  r.z * r.z;

            float invDist = 1.0f / sqrt(distSqr + epsSqr);

            float invDistCube = invDist * invDist * invDist;

            float s = localPos[j].w * invDistCube;



            // accumulate effect of all particles

            acc += s * r;

        }

        // Synchronize so that next tile can be loaded

        barrier(CLK_LOCAL_MEM_FENCE);

    }



    // updated position and velocity

    float4 newPos = oldPos + oldVel * deltaTime + acc * 0.5f * deltaTime * deltaTime;

    newPos.w = oldPos.w;



    float4 newVel = oldVel + acc * deltaTime;



    // write to global memory

    pos[gid] = newPos;

    vel[gid] = newVel;
}