Radix Sort Parallel Algorithm (n.elem > blockDim * blockDim)

#include "utils.h"

#define MAXTHREADS 1024
typedef unsigned int uint;
#define sizeuint sizeof( uint )


__global__ void krn_new_scan_1(
    const uint *d_inputVals,
    uint *d_tmpScan0, uint *d_tmpScan1,
    uint *d_tmpInterm0, uint *d_tmpInterm1,
    int numElems, int nbit)
{
    int tgrd = threadIdx.x + blockDim.x * blockIdx.x;

    if( tgrd >= numElems )
        return;

    int tblk = threadIdx.x;

    uint tmp0, tmp1;

    __shared__ uint shr_scan0[ MAXTHREADS ];
    __shared__ uint shr_scan1[ MAXTHREADS ];

    int bit = (d_inputVals[ tgrd ] >> nbit) & 1;

    shr_scan0[ tblk ] = !bit;
    shr_scan1[ tblk ] = bit;
    __syncthreads();

    // inclusive scan
    for( int i = 1; i <= MAXTHREADS / 2; i <<= 1 )
    {
        if( tblk >= i )
        {
            tmp0 = shr_scan0[ tblk - i ];
            tmp1 = shr_scan1[ tblk - i ];
        }
        __syncthreads();

        if( tblk >= i )
        {
            shr_scan0[ tblk ] += tmp0;
            shr_scan1[ tblk ] += tmp1;
        }
        __syncthreads();
    }

    // change to exclusive prefix sum
    d_tmpScan0[ tgrd ] = tblk ? shr_scan0[ tblk - 1 ] : 0;
    d_tmpScan1[ tgrd ] = tblk ? shr_scan1[ tblk - 1 ] : 0;

    if( tblk == 0 ) // choosing one thread (0) to do this processing
    {
        int ult_elem;

        if( blockIdx.x == gridDim.x - 1 )
            ult_elem = (numElems - 1) % blockDim.x;
        else
            ult_elem = blockDim.x - 1;

        // store the offset at each block; preparing to do their own prefix sum
        d_tmpInterm0[ blockIdx.x ] = shr_scan0[ ult_elem ];
        d_tmpInterm1[ blockIdx.x ] = shr_scan1[ ult_elem ];
    }
}


__global__ void krn_new_scan_2_1(
    uint *d_tmpInterm0, uint *d_tmpInterm1,
    uint *d_tmpIntermAux0, uint *d_tmpIntermAux1,
    int tmpIntermSize )
{
    int tgrd = threadIdx.x + blockDim.x * blockIdx.x;

    if( tgrd >= tmpIntermSize )
        return;

    int tblk = threadIdx.x;

    uint tmp0, tmp1;

    __shared__ uint shr_interm0[ MAXTHREADS ];
    __shared__ uint shr_interm1[ MAXTHREADS ];

    shr_interm0[ tblk ] = d_tmpInterm0[ tgrd ];
    shr_interm1[ tblk ] = d_tmpInterm1[ tgrd ];
    __syncthreads();

    // inclusive scan
    for( int i = 1; i <= MAXTHREADS / 2; i <<= 1 )
    {
        if( tblk >= i )
        {
            tmp0 = shr_interm0[ tblk - i ];
            tmp1 = shr_interm1[ tblk - i ];
        }
        __syncthreads();

        if( tblk >= i )
        {
            shr_interm0[ tblk ] += tmp0;
            shr_interm1[ tblk ] += tmp1;
        }
        __syncthreads();
    }

    // change to exclusive prefix sum
    d_tmpInterm0[ tgrd ] = tblk ? shr_interm0[ tblk - 1 ] : 0;
    d_tmpInterm1[ tgrd ] = tblk ? shr_interm1[ tblk - 1 ] : 0;


    if (tblk == 0) // randomly choosing just one thread (0) to do the processing
    {
        int ult_elem;

        if( blockIdx.x == gridDim.x - 1 )
            ult_elem = (tmpIntermSize - 1) % blockDim.x;
        else
            ult_elem = blockDim.x - 1;

        // store the offset at each block; preparing to do their own prefix sum
        d_tmpIntermAux0[ blockIdx.x ] = shr_interm0[ ult_elem ];
        d_tmpIntermAux1[ blockIdx.x ] = shr_interm1[ ult_elem ];
    }
}


__global__ void krn_new_scan_2_2(
    uint *d_tmpIntermAux0, uint *d_tmpIntermAux1)
{
    int tblk = threadIdx.x; // assuming grid dimension = 1 (only 1 block)

    if (tblk >= blockDim.x)
        return;

    uint tmp0, tmp1;

    __shared__ uint shr_interm0[ MAXTHREADS ];
    __shared__ uint shr_interm1[ MAXTHREADS ];

    shr_interm0[ tblk ] = d_tmpIntermAux0[ tblk ];
    shr_interm1[ tblk ] = d_tmpIntermAux1[ tblk ];
    __syncthreads();

    // inclusive scan
    for( int i = 1; i <= MAXTHREADS / 2; i <<= 1 )
    {
        if( tblk >= i )
        {
            tmp0 = shr_interm0[ tblk - i ];
            tmp1 = shr_interm1[ tblk - i ];
        }
        __syncthreads();

        if( tblk >= i )
        {
            shr_interm0[ tblk ] += tmp0;
            shr_interm1[ tblk ] += tmp1;
        }
        __syncthreads();
    }

    // change to exclusive prefix sum
    d_tmpIntermAux0[ tblk ] = tblk ? shr_interm0[ tblk - 1 ] : 0;
    d_tmpIntermAux1[ tblk ] = tblk ? shr_interm1[ tblk - 1 ] : 0;
    d_tmpIntermAux1[ tblk ] += shr_interm0[ blockDim.x - 1 ]; // add offset (# of bits 0) to bits 1
}


__global__ void krn_new_scan_2_3(
    uint *d_tmpInterm0, uint *d_tmpInterm1,
    uint *d_tmpIntermAux0, uint *d_tmpIntermAux1,
    int tmpIntermSize )
{
    int tgrd = threadIdx.x + blockDim.x * blockIdx.x;

    if (tgrd >= tmpIntermSize)
        return;

    d_tmpInterm0[ tgrd ] += d_tmpIntermAux0[ blockIdx.x ];
    d_tmpInterm1[ tgrd ] += d_tmpIntermAux1[ blockIdx.x ];
}


__global__ void krn_new_scan_3(
    uint *d_inputVals, uint *d_outputVals,
    uint *d_tmpScan0, uint *d_tmpScan1,
    uint *d_tmpInterm0, uint *d_tmpInterm1,
    int numElems, int nbit)
{

    int tgrd = threadIdx.x + blockDim.x * blockIdx.x;

    if( tgrd >= numElems )
        return;

    int pos;

    int bit = (d_inputVals[ tgrd ] >> nbit) & 1;

    if( !bit )
        pos = d_tmpScan0[ tgrd ] + d_tmpInterm0[ blockIdx.x ];
    else
        pos = d_tmpScan1[ tgrd ] + d_tmpInterm1[ blockIdx.x ];

    d_outputVals[ pos ] = d_inputVals[ tgrd ];
}



//=================================================================

void computeHistogram( unsigned int*  d_vals, unsigned int numElems )
{

    int blockDim = numBins; // 1024 (can't be > MAXTHREADS)
    int gridDim = (numElems + blockDim - 1) / blockDim;
    int subGridDim = (gridDim + blockDim - 1) / blockDim;


    uint *d_vals_out;
    uint *h_vals = new uint[ numElems ];
        checkCudaErrors( cudaMalloc( &d_vals_out, numElems * sizeuint ) );

    int size_tmpScan = numElems;
    uint *d_tmpScan0; uint *d_tmpScan1; uint *h_tmpScan0 = new uint[ size_tmpScan ]; uint *h_tmpScan1 = new uint[ size_tmpScan ];
        checkCudaErrors( cudaMalloc( &d_tmpScan0, size_tmpScan * sizeuint ) ); checkCudaErrors( cudaMalloc( &d_tmpScan1, size_tmpScan * sizeuint ) );

    int size_tmpInterm = gridDim;
    uint *d_tmpInterm0; uint *d_tmpInterm1; uint *h_tmpInterm0 = new uint[ size_tmpInterm ]; uint *h_tmpInterm1 = new uint[ size_tmpInterm ];
        checkCudaErrors( cudaMalloc( &d_tmpInterm0, size_tmpInterm * sizeuint ) ); checkCudaErrors( cudaMalloc( &d_tmpInterm1, size_tmpInterm * sizeuint ) );

    int size_tmpIntermAux = subGridDim;
    uint *d_tmpIntermAux0; uint *d_tmpIntermAux1; uint *h_tmpIntermAux0 = new uint[ size_tmpIntermAux ]; uint *h_tmpIntermAux1 = new uint[ size_tmpIntermAux ];
        checkCudaErrors( cudaMalloc( &d_tmpIntermAux0, size_tmpIntermAux * sizeuint ) ); checkCudaErrors( cudaMalloc( &d_tmpIntermAux1, size_tmpIntermAux * sizeuint ) );


    printf( "\n Starting. "
         "numElems= %d  numBins= %d  gridDim= %d  blockDim= %d\n",
         numElems,     numBins,     gridDim,     blockDim );

    const int bitsInUint = sizeuint * 8;

    uint *d_ptr_in = (uint*) d_vals;
    uint *d_ptr_out = d_vals_out;

    for( int nbit = 0; nbit < bitsInUint; nbit++ )
    {
        krn_new_scan_1 << < gridDim, blockDim >> > (
            d_ptr_in,
            d_tmpScan0, d_tmpScan1,
            d_tmpInterm0, d_tmpInterm1,
            numElems, nbit);
        cudaDeviceSynchronize(); checkCudaErrors( cudaGetLastError() );

        krn_new_scan_2_1 << < subGridDim, blockDim >> > (
            d_tmpInterm0, d_tmpInterm1,
            d_tmpIntermAux0, d_tmpIntermAux1,
            gridDim );
        cudaDeviceSynchronize(); checkCudaErrors( cudaGetLastError() );

        krn_new_scan_2_2 << < 1, subGridDim >> > (
            d_tmpIntermAux0, d_tmpIntermAux1);
        cudaDeviceSynchronize(); checkCudaErrors( cudaGetLastError() );

        krn_new_scan_2_3 << < subGridDim, blockDim >> > (
            d_tmpInterm0, d_tmpInterm1,
            d_tmpIntermAux0, d_tmpIntermAux1,
            gridDim);
        cudaDeviceSynchronize(); checkCudaErrors( cudaGetLastError() );

        krn_new_scan_3 << < gridDim, blockDim >> > (
            d_ptr_in, d_ptr_out,
            d_tmpScan0, d_tmpScan1,
            d_tmpInterm0, d_tmpInterm1,
            numElems, nbit);
        cudaDeviceSynchronize(); checkCudaErrors( cudaGetLastError() );

        if (nbit < bitsInUint - 1)
            std::swap( d_ptr_in, d_ptr_out );
    }

    // Here: sorted list is in d_ptr_out;
}