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

__constant float VERY_LARGE_NUMBER = 99999999; // A big number used fo scalling time steps;


__kernel void leapfrog_integrator(      __global float4* pos ,

                                                                        __global float4* vel,

                                                                        int numBodies,

                                                                    float deltaTime,
                
                                                                    float epsSqr,

                                                                    __local float4* localPos,
                                                                    
                                                                    __local float4* localVel,
                                                                    
                                                                    __global float4* acc ,
                                                                    
                                                                    int initialAccComputation,
                                                                   
                                                                    __global float* collTime)

{
        //---------------------------------------
    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 oldAcc = (float4)(0.0f, 0.0f, 0.0f, 0.0f);

        float time_scale_sq = VERY_LARGE_NUMBER;        // Init the collision time [read theory
                                                                                                // for more information;
                                                                                                
    // Compute the partial v_(n+1) and x_(n+1):
        //oldVel += acc[gid]*deltaTime/2;
        //oldPos += oldVel*deltaTime;
        
    for(int i = 0; i < numTiles; ++i)

    {
        // load one tile into local memory

        int idx = i * localSize + tid;

        localPos[tid] = pos[idx];
        localVel[tid] = vel[idx];

        // Synchronize to make sure data is available for processing

        barrier(CLK_LOCAL_MEM_FENCE);

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

        {
                        float4 rji, vji;
            float r2= 0.0f, r3= 0.0f, v2=0.0f;  // r^2, r^3, v^2;
                
                        rji = localPos[j] - oldPos;
                        vji = localVel[j] - oldVel;
           
                        r2 = pow(rji.x, 2) + pow(rji.y, 2) + pow(rji.z, 2);
                        v2 = pow(vji.x, 2) + pow(vji.y, 2) + pow(vji.z, 2);

                        if (r2 == 0)
                                r2 = 0.000000000001;
           
                        r3= r2*sqrt(r2);
                
                //----------------------
                        // The computation of acc
                        
                        oldAcc += oldPos.w*rji/r3;
                        
                        
                        //----------------------
                        // The computation of time_scale_sq:
                        
                        float estimate_sq= r2/v2;                       // [distance]^2/[velocity]^2 = [time]^2
                if (time_scale_sq > estimate_sq){
                        time_scale_sq = estimate_sq;
                        }
                        
                        float a = (oldPos.w + localPos[j].w)/r2;
                    estimate_sq = sqrt(r2)/a;           // [distance]/[acceleration] = [time]^2
                        
                        if (time_scale_sq > estimate_sq){
                        time_scale_sq = estimate_sq;
                        }
            }

        // Synchronize so that next tile can be loaded

        barrier(CLK_LOCAL_MEM_FENCE);
                
    }
    
    collTime[gid] = sqrt(time_scale_sq);
    
    // Now we compute the rest of the v_(n+1):
    //oldVel += oldAcc*deltaTime/2;
    
//      pos[gid] = oldPos;
//      vel[gid] = oldVel;
        acc[gid] = oldAcc;

}