clfRadixSort.cpp

// C++ class for sorting integer list in OpenCL
// copyright Philippe Helluy, University de Strasbourg, France, 2011, helluy@math.unistra.fr
// licensed under the GNU Lesser General Public License see http://www.gnu.org/copyleft/lesser.html
// if you find this software usefull you can cite the following work in your reports or articles:
// Philippe HELLUY, A portable implementation of the radix sort algorithm in OpenCL, HAL 2011.

// members of the class clfRadixSort
// see a description in the hpp...

#include <fstream>
#include <Base/Common/BIASpragma.hh>
#include <OpenCLFramework/Algorithm/clfRadixSort.hh>

#define TRANSPOSE 0  // transpose the initial vector (faster memory access)
////////////////////////////////////////////////////////


// the following parameters are computed from the previous
#define _RADIX (1 << 5) //  radix  = 2^5
#define _PASS (30/5) // number of needed passes to sort the list
#define _HISTOSIZE (16 * 16 * _RADIX ) // size of the histogram
// maximal value of integers for the sort to be correct
#define _MAXINT (1 << (30-1))


using namespace std; 
using namespace BIAS;

clfRadixSort::clfRadixSort(clfContext *ctx, bool sharedGL, unsigned int device)
: clfAlgorithm(ctx, sharedGL, device), nkeys((1<<23))
{
  nkeys_rounded=nkeys;
  // check some conditions
  assert(30 % 5 == 0);
  assert((1<<23) % (16 * 16) == 0);
  assert( (16 * 16 * _RADIX) % 512 == 0);
  assert(pow(float(2),(int) log2(float(16))) == 16);
  assert(pow(float(2),(int) log2(float(16))) == 16);

  // init the timers
  histo_time=0;
  scan_time=0;
  reorder_time=0;
  transpose_time=0;

  try {
    programCL_->AddSource("/OpenCLFramework/Algorithm/cl/radixSort.cl");
    programCL_->Build();
    
    programCL_->AddKernel("histogram");
    programCL_->AddKernel("scanhistograms");
    programCL_->AddKernel("pastehistograms");
    programCL_->AddKernel("reorder");
    programCL_->AddKernel("transpose");
    programCL_->AddKernel("applypermutationHighToLow");
    programCL_->AddKernel("applypermutation");
  
    d_inKeys = context_->CreateBuffer();
    d_inKeys->Allocate(sizeof(unsigned int)* (1<<23), false,false);
  
    d_outKeys = context_->CreateBuffer();
    d_outKeys->Allocate(sizeof(unsigned int)* (1<<23), false,false);
  
    d_inPermut = context_->CreateBuffer();
    d_inPermut->Allocate(sizeof(unsigned int)* (1<<23), false,false);
  
    d_outPermut = context_->CreateBuffer();
    d_outPermut->Allocate(sizeof(unsigned int)* (1<<23), false,false);
  
    // allocate the histogram on the GPU
    d_Histograms = context_->CreateBuffer();
    d_Histograms->Allocate(sizeof(unsigned int)* _RADIX * 16 * 16, false, false);
  
    d_globsum = context_->CreateBuffer();
    d_globsum->Allocate(sizeof(unsigned int)* 512, false, false);
    // allocate the auxiliary histogram on GPU
  
    d_temp = context_->CreateBuffer();
    d_temp->Allocate(sizeof(unsigned int)* 512, false, false);
    // temporary vector when the sum is not needed
  
    // we set here the fixed arguments of the OpenCL kernels
    // the changing arguments are modified elsewhere in the class
    
    programCL_->KernelSetArgument("histogram", 1, *d_Histograms);
    programCL_->KernelSetLocalArgument("histogram", 3, sizeof(unsigned int)*_RADIX*16);
  
    programCL_->KernelSetArgument("pastehistograms", 0, *d_Histograms);
    programCL_->KernelSetArgument("pastehistograms", 1, *d_globsum);
    
    programCL_->KernelSetArgument("reorder", 2, *d_Histograms);
    programCL_->KernelSetLocalArgument("reorder", 6, sizeof(unsigned int)* _RADIX * 16);

    int maxmemcache=max(512,16 * 16 * _RADIX / 512);
    programCL_->KernelSetLocalArgument("scanhistograms", 1, sizeof(unsigned int)* maxmemcache);

  } catch (clfException &e) {
    cout << e.GetDetailedString() << endl;
  }

  // construction of the initial permutation
  for(int i = 0; i < (1<<23); i++){
    h_initialPermut[i] = i;
  }
}

void clfRadixSort::SetData(int num, float *data, int scale) {
  Resize(num);
  for (int i=0;i<num;i++) {
    h_Keys[i] = (int)(data[i]*scale);
  }
  // construction of the initial permutation
  for(int i = 0; i < num; i++){
    h_Permut[i] = i;
  }
  Host2GPU();
}

void clfRadixSort::SetData(int num, clfBuffer *data) {
  Resize(num);
  data->CopyBuffer(*d_inKeys, 0,0, num*sizeof(unsigned int));
  d_inPermut->WriteToBuffer( h_initialPermut, 0, sizeof(unsigned int)  * num);
  context_->Finish();
}

// resize the sorted vector
void clfRadixSort::Resize(int nn){

  assert(nn <= (1<<23));

  nkeys=nn;

  // length of the vector has to be divisible by (16 * 16)
  int reste=nkeys % (16 * 16);
  nkeys_rounded=nkeys;
  unsigned int pad[16 * 16];
  for(int ii=0;ii<16 * 16;ii++){
    pad[ii]=_MAXINT-(unsigned int)1;
  }
  if (reste !=0) {
    nkeys_rounded=nkeys-reste+(16 * 16);
    // pad the vector with big values
    assert(nkeys_rounded <= (1<<23));
    d_inKeys->WriteToBuffer( pad, sizeof(unsigned int)*nkeys, sizeof(unsigned int) *(16 * 16 - reste));
  }

}

// transpose the list for faster memory access
void clfRadixSort::Transpose(int nbrow,int nbcol){

  int tilesize=32;

  // if the matrix is too small, avoid using local memory
  if (nbrow%tilesize != 0) tilesize=1;
  if (nbcol%tilesize != 0) tilesize=1;

  if (tilesize == 1) {
    cout << "Warning, small list, avoiding cache..."<<endl;
  }

  programCL_->KernelSetArgument("transpose", 0, *d_inKeys);
  programCL_->KernelSetArgument("transpose", 1, *d_outKeys);
  programCL_->KernelSetArgument("transpose", 2, nbcol);
  programCL_->KernelSetArgument("transpose", 3, nbrow);
  programCL_->KernelSetArgument("transpose", 4, *d_inPermut);
  programCL_->KernelSetArgument("transpose", 5, *d_outPermut);
  programCL_->KernelSetLocalArgument("transpose", 6, sizeof(unsigned int)*tilesize*tilesize);
  programCL_->KernelSetLocalArgument("transpose", 7, sizeof(unsigned int)*tilesize*tilesize);
  programCL_->KernelSetArgument("transpose", 8, tilesize);

  size_t global_work_size[2];
  size_t local_work_size[2];

  assert(nbrow%tilesize == 0);
  assert(nbcol%tilesize == 0);

  global_work_size[0]=nbrow/tilesize;
  global_work_size[1]=nbcol;

  local_work_size[0]=1;
  local_work_size[1]=tilesize;

  context_->RunOn2DRange(*programCL_, "transpose", global_work_size[0], global_work_size[1], local_work_size[0], local_work_size[1]);
  context_->Finish();
  //exchange the pointers

  // swap the old and new vectors of keys
  // swap the old and new vectors of keys
  clfBuffer *d_temp = d_inKeys;
  d_temp=d_inKeys;
  d_inKeys=d_outKeys;
  d_outKeys=d_temp;

  // swap the old and new permutations
  d_temp=d_inPermut;
  d_inPermut=d_outPermut;
  d_outPermut=d_temp;

}

// global sorting algorithm

void clfRadixSort::Sort(){

  assert(nkeys_rounded <= (1<<23));
  assert(nkeys <= nkeys_rounded);
  int nbcol=nkeys_rounded/(16 * 16);
  int nbrow= 16 * 16;

  if (TRANSPOSE){
    Transpose(nbrow,nbcol);
  }

  for(unsigned int pass=0;pass<_PASS;pass++){
    //for(unsigned int pass=0;pass<1;pass++){
    Histogram(pass);
    ScanHistogram();
    Reorder(pass);
  }

  if (TRANSPOSE){
    Transpose(nbcol,nbrow);
  }

  sort_time=histo_time+scan_time+reorder_time+transpose_time;
}

void clfRadixSort::ApplyPermutation(clfBuffer *in, clfBuffer *out, int elemsize, int numElems, bool hightolow) {
  string kname = "applypermutation";
  if (hightolow) kname = "applypermutationHighToLow";
  programCL_->KernelSetArgument(kname, 0, *in);
  programCL_->KernelSetArgument(kname, 1, *out);
  programCL_->KernelSetArgument(kname, 2, *d_inPermut);
  programCL_->KernelSetArgument(kname, 3, elemsize);
  programCL_->KernelSetArgument(kname, 4, numElems);
  context_->RunOn1DRange(*programCL_, kname, context_->DivUp(32, numElems), 32);
}

// check the computation at the end
void clfRadixSort::Check(){
  
  cout << "Get the data from the GPU"<<endl;

  RecupGPU();

  cout << "Test order"<<endl;

  // first see if the final list is ordered
  for(unsigned int i=0;i<nkeys-1;i++){
    if (!(h_Keys[i] <= h_Keys[i+1])) {
      cout <<"erreur tri "<< i<<" "<<h_Keys[i]<<" ,"<<i+1<<" "<<h_Keys[i+1]<<endl;
    }
    assert(h_Keys[i] <= h_Keys[i+1]);
  }

  cout << "test OK !"<<endl;

}

clfRadixSort::~clfRadixSort()
{
  delete d_inKeys;
  delete d_outKeys;
  delete d_Histograms;
  delete d_globsum;
  delete d_inPermut;
  delete d_outPermut;
}


// get the data from the GPU
void clfRadixSort::RecupGPU(void){
  d_inKeys->ReadFromBuffer( h_Keys, 0, sizeof(unsigned int)*(1<<23));
  d_inPermut->ReadFromBuffer( h_Permut, 0, sizeof(unsigned int)*(1<<23));
  d_Histograms->ReadFromBuffer( h_Histograms, 0, sizeof(unsigned int)  * _RADIX * 16 * 16);
  d_globsum->ReadFromBuffer( h_globsum, 0, sizeof(unsigned int)  * 512);
  context_->Finish();
}

// put the data to the GPU
void clfRadixSort::Host2GPU(void){
  d_inKeys->WriteToBuffer( h_Keys, 0, sizeof(unsigned int)  * (1<<23));
  d_inPermut->WriteToBuffer( h_Permut, 0, sizeof(unsigned int)  * (1<<23));
  context_->Finish();
}

unsigned int *clfRadixSort::GetPermutation(unsigned int num) {
  if (num == 0)
    d_inPermut->ReadFromBuffer( h_Permut, 0, sizeof(unsigned int)*(1<<23));
  else {
    d_inPermut->ReadFromBuffer( h_Permut, 0, sizeof(unsigned int)*num);
  }
  return h_Permut;
}


// compute the histograms
void clfRadixSort::Histogram(unsigned int pass){

  size_t nblocitems=16;
  size_t nbitems=16*16;

  assert(_RADIX == pow(float(2),int(5)));
  
  programCL_->KernelSetArgument("histogram", 0, *d_inKeys);
  programCL_->KernelSetArgument("histogram", 2, (int)pass);

  assert( nkeys_rounded%(16 * 16) == 0);
  assert( nkeys_rounded <= (1<<23));

  programCL_->KernelSetArgument("histogram", 4, (int)nkeys_rounded);

  context_->RunOn1DRange(*programCL_, "histogram", nbitems, nblocitems);
  context_->Finish();
}

// scan the histograms
void clfRadixSort::ScanHistogram(void){

  // numbers of processors for the local scan
  // half the size of the local histograms
  size_t nbitems=_RADIX* 16*16 / 2;


  size_t nblocitems= nbitems/512 ;


  // scan locally the histogram (the histogram is split into several
  // parts that fit into the local memory)

  programCL_->KernelSetArgument("scanhistograms", 0, *d_Histograms);
  programCL_->KernelSetArgument("scanhistograms", 2, *d_globsum);

  context_->RunOn1DRange(*programCL_, "scanhistograms", nbitems, nblocitems);
  context_->Finish();

  programCL_->KernelSetArgument("scanhistograms", 0, *d_globsum);
  programCL_->KernelSetArgument("scanhistograms", 2, *d_temp);

  nbitems= 512 / 2;
  nblocitems=nbitems;
  //nblocitems=1;

  context_->RunOn1DRange(*programCL_, "scanhistograms", nbitems, nblocitems);
  context_->Finish();

  // loops again in order to paste together the local histograms
  nbitems = _RADIX* 16*16/2;
  nblocitems=nbitems/512;

  context_->RunOn1DRange(*programCL_, "pastehistograms", nbitems, nblocitems);
  context_->Finish();
}

// reorder the data from the scanned histogram
void clfRadixSort::Reorder(unsigned int pass){

  size_t nblocitems=16;
  size_t nbitems=16*16;

  programCL_->KernelSetArgument("reorder", 0, *d_inKeys);
  programCL_->KernelSetArgument("reorder", 1, *d_outKeys);
  programCL_->KernelSetArgument("reorder", 3, (int)pass);
  programCL_->KernelSetArgument("reorder", 4, *d_inPermut);
  programCL_->KernelSetArgument("reorder", 5, *d_outPermut);
  programCL_->KernelSetArgument("reorder", 7, (int)nkeys_rounded);
    
  assert(_RADIX == pow(float(2),int(5)));

  context_->RunOn1DRange(*programCL_, "reorder", nbitems, nblocitems);
  context_->Finish();

  // swap the old and new vectors of keys
  clfBuffer *d_temp = d_inKeys;
  d_temp=d_inKeys;
  d_inKeys=d_outKeys;
  d_outKeys=d_temp;

  // swap the old and new permutations
  d_temp=d_inPermut;
  d_inPermut=d_outPermut;
  d_outPermut=d_temp;

}