/****   Copyright (c) 2001-2003  
 ****    Seung-Yeop Lee : Original MPI MessagePassingCode Code
 ****    Alexandros Gerbessiotis: RMA adaptation.   
 ****   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

// ** 4-Round Radix Sort using MPI_Isend/Irecv *************
void sprdx(unsigned int *Data, int Num)
{
  unsigned int *gData, *gNewData, *tranData;
  int *gOrder, *GbOrder, *grOrder,*gpOrder; 
  int *GtemOrder, *recTemp, *SenCount, *RevCount;
  int Bucket[BucketSize]; 
  int Global[BucketSize];
  int Order[BucketSize];
  int temOrder[BucketSize]; 
  int TabOrder[BucketSize];
  int ProToGo, OffInPro, GloOrder, ShiftBits;  
  register int Sum, i, j, k, l, m, r, w, Start;
  int Flag, count;
  int recv_count;
  MPI_Status status;
  int myReq;
  MPI_Request request[BucketSize];
  int *Table[16];
  int a1, a2, a3, a4, a5;
  int *GloTab[16];
  int GloUp[BucketSize];
  int loop;
  int myrank, ProcNum;
  myrank=AIPID();
  ProcNum=AINPROCS();


  recTemp = (int *)malloc(Num * sizeof(int));
  gData = (int *)malloc(Num * sizeof(int));
  gNewData = (int *)malloc(Num * sizeof(int));
  gOrder = (int *)malloc(ProcNum * sizeof(int));
  grOrder = (int *)malloc(ProcNum * sizeof(int));
  gpOrder = (int *)malloc(ProcNum * sizeof(int));
  GbOrder = (int *)malloc(ProcNum * sizeof(int));
  GtemOrder = (int *)malloc(ProcNum * sizeof(int));
  SenCount = (int *)malloc(ProcNum * sizeof(int));
  RevCount = (int *)malloc(ProcNum * sizeof(int));
  for (i=0; i< ProcNum; i++){
     Table[i] = (int *)malloc(BucketSize * sizeof(int));
     GloTab[i] = (int *)malloc(BucketSize * sizeof(int));
  }

  memset(gData, 0, Num * sizeof(int));
  memset(gNewData, 0, Num * sizeof(int));

// ******** Sort begins here ***************************************
  for (ShiftBits = 0; ShiftBits < NumSize; ShiftBits += ShiftSize) {
// ********* Initialize **********************  
    memset(Bucket, 0, BucketSize * sizeof(int));
    memset(GloUp, 0, BucketSize * sizeof(int));
    memset(gOrder, 0, ProcNum * sizeof(int));
    memset(grOrder, 0, ProcNum * sizeof(int));
    memset(gpOrder, 0, ProcNum * sizeof(int));
    memset(GbOrder, 0, ProcNum * sizeof(int));
    memset(SenCount, 0, ProcNum * sizeof(int));
    memset(RevCount, 0, ProcNum * sizeof(int));

// ********* Count the data and add its occurrence number to the Bucket
    for (i = 0; i < Num; i++)
      Bucket[(Data[i] >> ShiftBits) & Mask]++;  

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

    MPI_Allreduce(Bucket, Global, BucketSize, MPI_INT, 
                     MPI_SUM, MPI_COMM_WORLD);

    MPI_Scan(Bucket, Order, BucketSize, MPI_INT,
                MPI_SUM, MPI_COMM_WORLD);
    Sum = 0;
    for (i = 0; i < BucketSize; i++) {
      Order[i] = Order[i] + Sum - Bucket[i];
      Sum = Sum + Global[i];
    }
  
   memcpy(temOrder, Order, BucketSize * sizeof(int));
   memcpy(TabOrder, Order, BucketSize * sizeof(int));

// *** Counting numbers of send and receive for each processor 
    for (i = 0; i < Num; i++) {
        j = (Data[i] >> ShiftBits) & Mask;
        GloOrder = temOrder[j]++;
        ProToGo = GloOrder / Num;

        if ((myrank != ProToGo) && ((GloOrder % Num) == 0) && 
                            (Order[j] != (temOrder[j]-1)))
              GloUp[j] = GloOrder; // if one bucket can't go to
                                   // one processor

        gOrder[ProToGo]++;
        if (ProToGo != myrank) 
            SenCount[ProToGo] = 1; // count the number of send for 
                                   // each processor
    }  
    MPI_Allreduce(SenCount, RevCount, ProcNum, MPI_INT,
                MPI_SUM, MPI_COMM_WORLD);
    myReq = RevCount[myrank];     // count the number of receive
                                  // for each processor
    MPI_Scan(gOrder, gpOrder, ProcNum, MPI_INT,
                MPI_SUM, MPI_COMM_WORLD);
    AIOREGISTER((char *)grOrder,ProcNum*sizeof(int));
    AIOINIT((char *)grOrder);
    AIOCOMMIT();
    for(i=0;i<ProcNum;i++) {
        AIOHPPUT(i,(char *)&gOrder[i],(char *)&grOrder[0],myrank*sizeof(int),sizeof(int));
    }
   AIOSYNC((char *)grOrder);
   AIBARRIER();
   AIODEREGISTER((char *)grOrder);


// *** Broadcasr Global Order Table ***********
    GbOrder[0] = gOrder[0];
    for (i = 1; i < ProcNum; i++) 
        GbOrder[i] = GbOrder[i-1] + gOrder[i];

    memcpy(GtemOrder, GbOrder, ProcNum * sizeof(int));

  MPI_Barrier(MPI_COMM_WORLD);
    for (i = 0; i < ProcNum; i++) {
         if (myrank == i) {  
            Table[i] =  TabOrder;
            GloTab[i] = GloUp;    
         }
         MPI_Bcast(Table[i], BucketSize, MPI_INT,
                     i, MPI_COMM_WORLD);
         MPI_Bcast(GloTab[i], BucketSize, MPI_INT,
                     i, MPI_COMM_WORLD);  // every processor knows
    }                                     // global order for all 
                                          // other processors
    
// ******* Rearrange data by processor order *******************
    memcpy(temOrder, Order, BucketSize * sizeof(int));

    Start = 0;
    for (i = 0; i < Num; i++) {
        GloOrder = temOrder[(Data[i] >> ShiftBits) & Mask]++;
        ProToGo = GloOrder / Num;
        OffInPro = GloOrder % Num;

        if (ProToGo == myrank) // ***********Local Data
           gNewData[OffInPro] = Data[i];

        else {
           if (ProToGo == 0) 
               gData[Start++] = Data[i];
           else 
               gData[GtemOrder[ProToGo-1]++] = Data[i];
        }
    }

//****** Send data **********************************************
    memset(Data, 0, Num * sizeof(int));
    memcpy(temOrder, Order, BucketSize * sizeof(int));
    memcpy(GtemOrder, GbOrder, ProcNum * sizeof(int));
    w = 0;
    recv_count = Num;
    i = myrank;
    loop = 0;
    AIOREGISTER((char *)recTemp,Num*sizeof(int));
    AIOINIT((char *)recTemp);
    AIOCOMMIT();

    while (loop < ProcNum) { 
        if (gOrder[i] != 0){
          if (i!=myrank){
                if (i == 0) 
		    AIOHPPUT(i,(char *)&gData[i],(char *)recTemp,(gpOrder[i]-gOrder[i])*sizeof(int),gOrder[i]*sizeof(int));
                else 
		    AIOHPPUT(i,(char *)&gData[GtemOrder[i-1]],(char *)recTemp,(gpOrder[i]-gOrder[i])*sizeof(int),gOrder[i]*sizeof(int));
        }
       }
        i = (i+1)%ProcNum;
        loop++;
   }
   AIOSYNC((char *)recTemp);
   AIBARRIER();
   AIODEREGISTER((char *)recTemp);
//****** Receive Data by Global Order Table ********************** 
     gpOrder[0]=0;
     for(i=1;i<ProcNum;i++){
      gpOrder[i]=gpOrder[i-1]+grOrder[i-1];
     }
     for(i=0;i<ProcNum;i++){
      if (i!= myrank){
        a2 = gpOrder[i];
        for (j=0; j < grOrder[i]; j++){
               a1 = (recTemp[a2+j] >> ShiftBits) & Mask;
               if (GloTab[i][a1] != 0){ 
                   if ((GloTab[i ][a1] / Num) == myrank) {
                          Table[i ][a1] = GloTab[i ][a1];
                          GloTab[i ][a1] = 0;
                   }  
               }
                 GloOrder = Table[i ][a1]++;
                        ProToGo = GloOrder / Num;
                        OffInPro = GloOrder % Num;
                 gNewData[OffInPro] = recTemp[a2+j];
          }
      }
     }
    memcpy(Data, gNewData, Num * sizeof(int));
 }// end of sorting 
}

