/****   Copyright (c) 2001-2003  
 ****    Seung-Yeop Lee : Original MPI MessagePassingCode Code
 ****    Alexandros Gerbessiotis:  Converted RemoteMemoryAccess extension
 ****   Permission to use, copy, modify, and distribute this software,
 ****   and to incorporate it, in whole or in part, into other software,
 ****   is hereby granted without fee, provided that
 ****     (1) the above copyright notice and this permission notice appear in
 ****         all copies of the source code, and the above copyright notice
 ****         appear in clearly visible form on all supporting documentation
 ****         and distribution media, and
 ****     (2) any redistribution of the software, in original
 ****         or modified form, be without fee and subject to
 ****         these same conditions.
 ****   No guarantee is offered that the code works as
 ****   advertised or is absent of any, even damaging,
 ****   side-effects; use this code on your own personal risk
 ****/
#include	"rdxai.h"	
#include "ai.h"
#define  BucketSize   256
#define  Mask         BucketSize - 1
#define  ShiftSize    8
#define  NumSize      32
#define LOOPOPT
#undef LOOPOPT
#ifdef  LOOPOPT
#define NNM  100000
#define NN   16
#endif


void sbrdx(unsigned int *Data, int Num)
{
  unsigned int *locData;    // local data storage
  int Bucket[BucketSize];   // Bucket
  int Global[BucketSize];
  int Order[BucketSize];    // Global Order
  int LocOrder[BucketSize];
  int BucOrder[BucketSize];
  int Rcount; // Number of data to receive
  register int ProToGo, OffInPro, GloOrder, ShiftBits;
  register int Sum, i, j, k, m,  Flag, Start;
  register int Buc, Size;
  int myrank, ProcNum;
#ifdef LOOPOPT
  int mpistat1;
#endif
  myrank=AIPID();
  ProcNum=AINPROCS();

  locData = (int *)malloc(Num * sizeof(int));
  memset(locData, 0, Num* sizeof(unsigned int));

 AIOREGISTER((char *)Order,BucketSize*sizeof(int));
 AIOREGISTER((char *)Global,BucketSize*sizeof(int));
 AIOCOMMIT();

// ** Sort begins here ***************************************
  for (ShiftBits = 0; ShiftBits < NumSize; ShiftBits += ShiftSize) {

    memset(Bucket, 0, BucketSize * sizeof(int));
    for (i = 0; i < Num; i++)
       Bucket[(Data[i] >> ShiftBits) & Mask]++; // count

// ** Local sort within each processor (count sort) **********
    LocOrder[0] = Bucket[0];
    for (i = 1; i < BucketSize; i++)
       LocOrder[i] = LocOrder[i-1] + Bucket[i];

    memcpy(BucOrder, LocOrder, BucketSize * sizeof(int));

    Start = 0;
    for (i = 0; i < Num; i++) {
      j = (Data[i] >> ShiftBits) & Mask;
      if (j == 0)
          locData[Start++] = Data[i];
      else
          locData[LocOrder[j-1]++] = Data[i];
    }

// ** Compute the global order *****************************

/*
*/
#ifdef MPILIB
    MPI_Allreduce(Bucket, Global, BucketSize, MPI_LONG,
                     MPI_SUM, MPI_COMM_WORLD);

    MPI_Scan(Bucket, Order, BucketSize, MPI_LONG,
                MPI_SUM, MPI_COMM_WORLD);
#else
    AIOINIT((char *)Order);
    AIOINIT((char *)Global);
    ai2scan(oper_add,BucketSize,(char *)Bucket,(char *)Order,(char *)Global,sizeof(int));
#endif

    Sum = 0;
    for (i = 0; i < BucketSize; i++) {
      Order[i] = Order[i] + Sum - Bucket[i];
      Sum = Sum + Global[i];
    }
    Rcount = Num;

    AIOREGISTER((char *)Data,Num*sizeof(unsigned int));
    AIOINIT((char *)Data);
    AIOCOMMIT();
// ** Sort Data globally by Global Order **************
#ifdef LOOPOPT
    mpistat1=(BucketSize/ProcNum)*myrank;
    for (j=0; j < BucketSize; j++) {
        i=(mpistat1+j)%BucketSize;
#else
    for (i=0; i < BucketSize; i++) {
#endif
        Size = Bucket[i];
        if (Size > 0) { // ** Bucket is not empty *******
           if (i == 0)
               Buc = 0;
           else
               Buc = BucOrder[i-1];
           GloOrder = Order[(locData[Buc] >> ShiftBits) & Mask];
           ProToGo = GloOrder / Num;
           OffInPro = GloOrder % Num;

           k = Num - OffInPro; //How much space is left?
           // *** send One Bucket at once ********
           if (Size <= k) {
               if (ProToGo == myrank) { // ** all Local data **********
                   memcpy(&Data[OffInPro], &locData[Buc],
                              Size * sizeof(unsigned int));
                   Rcount -= Size;
               }
               else { // ** all Non-local data ***********************
/*
                   m = locData[Buc-1];
                   locData[Buc-1] = OffInPro;
                   MPI_Isend(&locData[Buc-1], Size+1, MPI_INT, ProToGo,
                                    1, MPI_COMM_WORLD, &request);
                   locData[Buc-1] = m;
                   check(Data, Num, &Rcount);
*/
                   AIOHPPUT(ProToGo,(char *)&locData[Buc],(char *)Data,
                           OffInPro*sizeof(unsigned int),Size*sizeof(unsigned int));
               }
           }
           // ** One Bucket can be sent twice *******************
           else { //**** Size > k
               if (ProToGo == myrank) { // ** Local data ***************
                   memcpy(&Data[OffInPro], &locData[Buc],
                              k * sizeof(unsigned int));
                   Rcount -= k;
                   // ********************* Non-local data *************
                   GloOrder += k;
                   ProToGo = GloOrder / Num;
                   OffInPro = GloOrder % Num;

/*
                   m = locData[Buc+k-1];
                   locData[Buc+k-1] = OffInPro;
                   MPI_Isend(&locData[Buc+k-1], Size-k+1, MPI_INT, ProToGo,
                                    1, MPI_COMM_WORLD, &request);
                   locData[Buc+k-1] = m;
                   check(Data, Num, &Rcount);
*/
                   AIOHPPUT(ProToGo,(char *)&locData[Buc+k],(char *)Data,
                           OffInPro*sizeof(unsigned int),(Size-k)*sizeof(unsigned int));
               }
               else { // ** All Non-local data to Processor A  ********
/*
                   m = locData[Buc-1];
                   locData[Buc-1] = OffInPro;
                   MPI_Isend(&locData[Buc-1], k+1, MPI_INT, ProToGo,
                                    1, MPI_COMM_WORLD, &request);
                   locData[Buc-1] = m;
*/
                   AIOHPPUT(ProToGo,(char *)&locData[Buc],(char *)Data,
                           OffInPro*sizeof(unsigned int),(k)*sizeof(unsigned int));
                   // ************* send to next to A processor
                   GloOrder += k;
                   ProToGo = GloOrder / Num;
                   OffInPro = GloOrder % Num;

/*
                   m = locData[Buc+k-1];
                   locData[Buc+k-1] = OffInPro;
                   MPI_Isend(&locData[Buc+k-1], Size-k+1, MPI_INT, ProToGo,
                                    1, MPI_COMM_WORLD, &request);
                   locData[Bucket[i]+k-1] = m;
                   check(Data, Num, &Rcount);
*/
                   AIOHPPUT(ProToGo,(char *)&locData[Buc+k],(char *)Data,
                           OffInPro*sizeof(unsigned int),(Size-k)*sizeof(unsigned int));
               }
           } // else (Size > k)
        } // if (Size>0)
    } // for
    AIOSYNC((char *)Data);
    AIBARRIER();
    AIODEREGISTER((char *)Data);
  } // for
  free((void *)locData);
} // radix_sort

