TreeBcast.hpp

#ifndef _PEXSI_TREE_HPP_
#define _PEXSI_TREE_HPP_

#include "pexsi/environment.hpp"
#include "pexsi/timer.h"

#include <vector>
#include <map>
#include <algorithm>
#include <string>
//#include <random>

// options to switch from a flat bcast/reduce tree to a binary tree

#ifndef FTREE_LIMIT
#define FTREE_LIMIT 16
#endif



namespace PEXSI{


extern std::map< MPI_Comm , std::vector<int> > commGlobRanks;

#ifdef NEW_BCAST

template< typename T>
  class TreeBcast2{
    protected:

      T * myData_;
      MPI_Request recvRequest_;
      NumVec<char> myRecvBuffer_;

      NumVec<MPI_Request> myRequests_;
      NumVec<MPI_Status> myStatuses_;

      bool done_;
      bool fwded_;
      //    bool isAllocated_;

      Int myRoot_;
      MPI_Comm comm_;
      vector<Int> myDests_;
      Int myRank_;
      Int msgSize_;
      bool isReady_;
      Int mainRoot_;
      Int tag_;
      Int numRecv_;

#ifdef COMM_PROFILE_BCAST
    protected:
      Int myGRoot_;
      Int myGRank_;
      //vector<int> Granks_;
    public:
      void SetGlobalComm(const MPI_Comm & pGComm){
        if(commGlobRanks.count(comm_)==0){
          MPI_Group group2 = MPI_GROUP_NULL;
          MPI_Comm_group(pGComm, &group2);
          MPI_Group group1 = MPI_GROUP_NULL;
          MPI_Comm_group(comm_, &group1);

          Int size;
          MPI_Comm_size(comm_,&size);
          vector<int> globRanks(size);
          vector<int> Lranks(size);
          for(int i = 0; i<size;++i){Lranks[i]=i;}
          MPI_Group_translate_ranks(group1, size, &Lranks[0],group2, &globRanks[0]);
          commGlobRanks[comm_] = globRanks;
        }
        myGRoot_ = commGlobRanks[comm_][myRoot_];
        myGRank_ = commGlobRanks[comm_][myRank_];
        //Granks_.resize(myDests_.size());
        //for(int i = 0; i<myDests_.size();++i){
        //  Granks_[i] = globRanks[myDests_[i]];
        //}

        //statusOFS<<myDests_<<std::endl;
        //statusOFS<<Granks_<<std::endl;
      }
#endif



    protected:
      virtual void buildTree(Int * ranks, Int rank_cnt)=0;





    public:

      TreeBcast2(){
        comm_ = MPI_COMM_WORLD;
        myRank_=-1;
        myRoot_ = -1; 
        msgSize_ = -1;
        numRecv_ = -1;
        tag_=-1;
        mainRoot_=-1;
        isReady_ = false;
        myData_ = NULL;
        recvRequest_ = MPI_REQUEST_NULL;
        fwded_=false;
        //      isAllocated_=false;
        done_ = false;
      }


      TreeBcast2(const MPI_Comm & pComm, Int * ranks, Int rank_cnt,Int msgSize):TreeBcast2(){
        comm_ = pComm;
        MPI_Comm_rank(comm_,&myRank_);
        myRoot_ = -1; 
        msgSize_ = msgSize;
        numRecv_ = 0;
        tag_=-1;
        mainRoot_=ranks[0];
        isReady_ = false;
      }


      virtual TreeBcast2 * clone() const = 0; 

      TreeBcast2(const TreeBcast2 & Tree){
        this->Copy(Tree);
      }

      virtual void Copy(const TreeBcast2 & Tree){
        this->comm_ = Tree.comm_;
        this->myRank_ = Tree.myRank_;
        this->myRoot_ = Tree.myRoot_; 
        this->msgSize_ = Tree.msgSize_;

        this->numRecv_ = Tree.numRecv_;
        this->tag_= Tree.tag_;
        this->mainRoot_= Tree.mainRoot_;
        this->isReady_ = Tree.isReady_;
        this->myDests_ = Tree.myDests_;


        this->recvRequest_ = Tree.recvRequest_;
        this->myRecvBuffer_ = Tree.myRecvBuffer_;
        this->myRequests_ = Tree.myRequests_;
        this->myStatuses_ = Tree.myStatuses_;
        this->myData_ = Tree.myData_;
        if(Tree.myData_==(T*)&Tree.myRecvBuffer_[0]){
          this->myData_=(T*)&this->myRecvBuffer_[0];
        }




        this->fwded_= Tree.fwded_;
        //      this->isAllocated_= Tree.isAllocated_;
        this->done_= Tree.done_;
      }

      void Reset(){
        assert(done_);
        CleanupBuffers();
        done_=false;
        myData_ = NULL;
        recvRequest_ = MPI_REQUEST_NULL;
        fwded_=false;
        //      isAllocated_=false;
        isReady_=false;
        numRecv_ = 0;
      }

      //bool IsAllocated(){return isAllocated_;}

      virtual ~TreeBcast2(){
        CleanupBuffers();
      }


      static TreeBcast2<T> * Create(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize,double rseed);

      virtual inline Int GetNumRecvMsg(){return numRecv_;}
      virtual inline Int GetNumMsgToSend(){return GetDestCount();}
      inline void SetDataReady(bool rdy){ 
        isReady_=rdy;
        //numRecv_ = rdy?1:0;
      }
      inline void SetTag(Int tag){ tag_ = tag;}
      inline int GetTag(){ return tag_;}

      bool IsDone(){return done_;}
      bool IsDataReady(){return isReady_;}
      bool IsDataReceived(){return numRecv_==1;}

      Int * GetDests(){ return &myDests_[0];}
      Int GetDest(Int i){ return myDests_[i];}
      Int GetDestCount(){ return myDests_.size();}
      Int GetRoot(){ return myRoot_;}

      bool IsRoot(){ return myRoot_==myRank_;}
      Int GetMsgSize(){ return msgSize_;}

      void ForwardMessage( ){
        if(myRequests_.m()!=GetDestCount()){
          myRequests_.Resize(GetDestCount());
          SetValue(myRequests_,MPI_REQUEST_NULL);
        }
        for( Int idxRecv = 0; idxRecv < myDests_.size(); ++idxRecv ){
          Int iProc = myDests_[idxRecv];
          // Use Isend to send to multiple targets
          MPI_Isend( myData_, msgSize_, MPI_BYTE, 
              iProc, tag_,comm_, &myRequests_[idxRecv] );

#if ( _DEBUGlevel_ >= 1 ) || defined(BCAST_VERBOSE)
          statusOFS<<myRank_<<" FWD to "<<iProc<<" on tag "<<tag_<<std::endl;
#endif
#ifdef COMM_PROFILE_BCAST
//        statusOFS<<idxRecv<<std::endl;
//        statusOFS<<myDests_<<std::endl;
//        statusOFS<<Granks_<<std::endl;
          //PROFILE_COMM(myGRank_,Granks_[idxRecv],tag_,msgSize_);
          PROFILE_COMM(myGRank_,commGlobRanks[comm_][iProc],tag_,msgSize_);
#endif
        } // for (iProc)
        fwded_ = true;
      }

      void CleanupBuffers(){
        myRequests_.Clear();
        myStatuses_.Clear();
        myRecvBuffer_.Clear();
      }


      void SetLocalBuffer(T * locBuffer){
        if(myData_!=NULL && myData_!=locBuffer){
          if(numRecv_>0){
            CopyLocalBuffer(locBuffer);
          }
          if(!fwded_){
            myRecvBuffer_.Clear(); 
          }
        }

        myData_ = locBuffer;
      }

      //async wait and forward
      virtual bool Progress(){
        if(done_){
          return true;
        }

        bool done = false;

        if (myRank_==myRoot_){
          if(isReady_){
            if(!fwded_){
#if ( _DEBUGlevel_ >= 1 ) || defined(BCAST_VERBOSE)
              statusOFS<<myRank_<<" FORWARDING on tag "<<tag_<<std::endl;
#endif
              ForwardMessage();
            }
            else{

              if(myStatuses_.m()!=GetDestCount()){
                myStatuses_.Resize(GetDestCount());
                recvRequest_ = MPI_REQUEST_NULL;
              }
              //test the send requests
              int flag = 0;
              int reqCnt = GetDestCount();
              if(reqCnt>0){
                assert(reqCnt == myRequests_.m());
                MPI_Testall(reqCnt,&myRequests_[0],&flag,&myStatuses_[0]);
                done = flag==1;
              }
              else{
                done=true;
              }
            }
          }
        }
        else{
          bool received = (numRecv_==1);
          
          if(!received){
            if(recvRequest_ == MPI_REQUEST_NULL ){
#if ( _DEBUGlevel_ >= 1 ) || defined(BCAST_VERBOSE)
              statusOFS<<myRank_<<" POSTING RECV on tag "<<tag_<<std::endl;
#endif
              //post the recv
              PostRecv();
            }
            else{

              if(myStatuses_.m()!=GetDestCount()){
                myStatuses_.Resize(GetDestCount());
                recvRequest_ = MPI_REQUEST_NULL;
              }
#if ( _DEBUGlevel_ >= 1 ) || defined(BCAST_VERBOSE)
              statusOFS<<myRank_<<" TESTING RECV on tag "<<tag_<<std::endl;
#endif
              //test
              int flag = 0;
              MPI_Status stat;
              int test = MPI_Test(&recvRequest_,&flag,&stat);
              assert(test==MPI_SUCCESS);

              if(flag==1){
                numRecv_=1;
                received = true;

                if(!fwded_){
#if ( _DEBUGlevel_ >= 1 ) || defined(BCAST_VERBOSE)
                  statusOFS<<myRank_<<" FORWARDING on tag "<<tag_<<std::endl;
#endif
                  ForwardMessage();
                }
              }
            }
          }
          else {
            assert(fwded_);
            //test the send requests
            int flag = 0;
            int reqCnt = GetDestCount();
            if(reqCnt>0){
              assert(reqCnt == myRequests_.m());
              MPI_Testall(reqCnt,&myRequests_[0],&flag,&myStatuses_[0]);
              done = flag==1;
            }
            else{
              done=true;
            }
          }
        }

        if(done){
          //free the unnecessary arrays
          myRequests_.Clear();
          myStatuses_.Clear();
#if ( _DEBUGlevel_ >= 1 ) || defined(BCAST_VERBOSE)
          statusOFS<<myRank_<<" EVERYTHING COMPLETED on tag "<<tag_<<std::endl;
#endif
        }

        done_ = done;

        return done;
      }

      //blocking wait
      void Wait(){
        if(!done_){
          while(!Progress());
        }
      }

      T * GetLocalBuffer(){
        return myData_;
      }

      virtual void PostRecv()
      {
        if(this->numRecv_<1 && this->recvRequest_==MPI_REQUEST_NULL && myRank_!=myRoot_){

          if(myData_==NULL){
            myRecvBuffer_.Resize(msgSize_);
            myData_ = (T*)&myRecvBuffer_[0];
          }
          MPI_Irecv( (char*)this->myData_, this->msgSize_, MPI_BYTE, 
              this->myRoot_, this->tag_,this->comm_, &this->recvRequest_ );
        }
      }



      void CopyLocalBuffer(T* destBuffer){
        std::copy((char*)myData_,(char*)myData_+GetMsgSize(),(char*)destBuffer);
      }




  };


template< typename T>
class FTreeBcast2: public TreeBcast2<T>{
  protected:
    virtual void buildTree(Int * ranks, Int rank_cnt){

      Int idxStart = 0;
      Int idxEnd = rank_cnt;



      this->myRoot_ = ranks[0];

      if(this->myRank_==this->myRoot_){
        this->myDests_.insert(this->myDests_.begin(),&ranks[1],&ranks[0]+rank_cnt);
      }

#if (defined(BCAST_VERBOSE)) 
      statusOFS<<"My root is "<<this->myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<this->myDests_.size();++i){statusOFS<<this->myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }



  public:
    FTreeBcast2(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeBcast2<T>(pComm,ranks,rank_cnt,msgSize){
      //build the binary tree;
      buildTree(ranks,rank_cnt);
    }


    virtual FTreeBcast2 * clone() const{
      FTreeBcast2 * out = new FTreeBcast2(*this);
      return out;
    } 
};

template< typename T>
class BTreeBcast2: public TreeBcast2<T>{
  protected:
    virtual void buildTree(Int * ranks, Int rank_cnt){

      Int idxStart = 0;
      Int idxEnd = rank_cnt;



      Int prevRoot = ranks[0];
      while(idxStart<idxEnd){
        Int curRoot = ranks[idxStart];
        Int listSize = idxEnd - idxStart;

        if(listSize == 1){
          if(curRoot == this->myRank_){
            this->myRoot_ = prevRoot;
            break;
          }
        }
        else{
          Int halfList = floor(ceil(double(listSize) / 2.0));
          Int idxStartL = idxStart+1;
          Int idxStartH = idxStart+halfList;

          if(curRoot == this->myRank_){
            if ((idxEnd - idxStartH) > 0 && (idxStartH - idxStartL)>0){
              Int childL = ranks[idxStartL];
              Int childR = ranks[idxStartH];

              this->myDests_.push_back(childL);
              this->myDests_.push_back(childR);
            }
            else if ((idxEnd - idxStartH) > 0){
              Int childR = ranks[idxStartH];
              this->myDests_.push_back(childR);
            }
            else{
              Int childL = ranks[idxStartL];
              this->myDests_.push_back(childL);
            }
            this->myRoot_ = prevRoot;
            break;
          } 

          if( this->myRank_ < ranks[idxStartH]){
            idxStart = idxStartL;
            idxEnd = idxStartH;
          }
          else{
            idxStart = idxStartH;
          }
          prevRoot = curRoot;
        }

      }

#if (defined(BCAST_VERBOSE))
      statusOFS<<"My root is "<<myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<this->myDests_.size();++i){statusOFS<<this->myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }



  public:
    BTreeBcast2(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeBcast2<T>(pComm,ranks,rank_cnt,msgSize){
      //build the binary tree;
      buildTree(ranks,rank_cnt);
    }

    virtual BTreeBcast2<T> * clone() const{
      BTreeBcast2<T> * out = new BTreeBcast2<T>(*this);
      return out;
    }

};



template< typename T>
class ModBTreeBcast2: public TreeBcast2<T>{
  protected:
    double rseed_;

    virtual void buildTree(Int * ranks, Int rank_cnt){

      Int idxStart = 0;
      Int idxEnd = rank_cnt;

      //sort the ranks with the modulo like operation
      if(rank_cnt>1){
        //Int new_idx = (int)((rand()+1.0) * (double)rank_cnt / ((double)RAND_MAX+1.0));

//        srand(ranks[0]+rank_cnt);
        Int new_idx = (Int)rseed_ % (rank_cnt - 1) + 1; 
        //Int new_idx = (int)((rank_cnt - 0) * ( (double)this->rseed_ / (double)RAND_MAX ) + 0);// (this->rseed_)%(rank_cnt-1)+1;
        //statusOFS<<"NEW IDX: "<<new_idx<<endl;



        Int * new_start = &ranks[new_idx];

        //for(int i =0;i<rank_cnt;++i){statusOFS<<ranks[i]<<" ";} statusOFS<<std::endl;

//        Int * new_start = std::lower_bound(&ranks[1],&ranks[0]+rank_cnt,ranks[0]);
        //just swap the two chunks   r[0] | r[1] --- r[new_start-1] | r[new_start] --- r[end]
        // becomes                   r[0] | r[new_start] --- r[end] | r[1] --- r[new_start-1] 
        std::rotate(&ranks[1], new_start, &ranks[0]+rank_cnt);

        //for(int i =0;i<rank_cnt;++i){statusOFS<<ranks[i]<<" ";} statusOFS<<std::endl;
      }

      Int prevRoot = ranks[0];
      while(idxStart<idxEnd){
        Int curRoot = ranks[idxStart];
        Int listSize = idxEnd - idxStart;

        if(listSize == 1){
          if(curRoot == this->myRank_){
            this->myRoot_ = prevRoot;
            break;
          }
        }
        else{
          Int halfList = floor(ceil(double(listSize) / 2.0));
          Int idxStartL = idxStart+1;
          Int idxStartH = idxStart+halfList;

          if(curRoot == this->myRank_){
            if ((idxEnd - idxStartH) > 0 && (idxStartH - idxStartL)>0){
              Int childL = ranks[idxStartL];
              Int childR = ranks[idxStartH];

              this->myDests_.push_back(childL);
              this->myDests_.push_back(childR);
            }
            else if ((idxEnd - idxStartH) > 0){
              Int childR = ranks[idxStartH];
              this->myDests_.push_back(childR);
            }
            else{
              Int childL = ranks[idxStartL];
              this->myDests_.push_back(childL);
            }
            this->myRoot_ = prevRoot;
            break;
          } 

          //not true anymore ?
          //first half to 
TIMER_START(FIND_RANK);
          Int * pos = std::find(&ranks[idxStartL], &ranks[idxStartH], this->myRank_);
TIMER_STOP(FIND_RANK);
          if( pos != &ranks[idxStartH]){
            idxStart = idxStartL;
            idxEnd = idxStartH;
          }
          else{
            idxStart = idxStartH;
          }
          prevRoot = curRoot;
        }

      }

#if (defined(REDUCE_VERBOSE))
      statusOFS<<"My root is "<<myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<this->myDests_.size();++i){statusOFS<<this->myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }



  public:
    ModBTreeBcast2(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize, double rseed):TreeBcast2<T>(pComm,ranks,rank_cnt,msgSize){
      //build the binary tree;
      rseed_ = rseed;
      buildTree(ranks,rank_cnt);
    }

    //virtual void Copy(const ModBTreeBcast & Tree){
    //  comm_ = Tree.comm_;
    //  myRank_ = Tree.myRank_;
    //  myRoot_ = Tree.myRoot_; 
    //  msgSize_ = Tree.msgSize_;

    //  numRecv_ = Tree.numRecv_;
    //  tag_= Tree.tag_;
    //  mainRoot_= Tree.mainRoot_;
    //  isReady_ = Tree.isReady_;
    //  myDests_ = Tree.myDests_;

    //  rseed_ = Tree.rseed_;
    //  myRank_ = Tree.myRank_;
    //  myRoot_ = Tree.myRoot_; 
    //  msgSize_ = Tree.msgSize_;

    //  numRecv_ = Tree.numRecv_;
    //  tag_= Tree.tag_;
    //  mainRoot_= Tree.mainRoot_;
    //  isReady_ = Tree.isReady_;
    //  myDests_ = Tree.myDests_;
    //}
 
    virtual ModBTreeBcast2 * clone() const{
      ModBTreeBcast2 * out = new ModBTreeBcast2(*this);
      return out;
    }

};



#endif

































class TreeBcast{
  protected:
    Int myRoot_;
    MPI_Comm comm_;
    vector<Int> myDests_;
    Int myRank_;
    Int msgSize_;
    bool isReady_;
    Int mainRoot_;
    Int tag_;
    Int numRecv_;


#if defined(COMM_PROFILE_BCAST) || defined(COMM_PROFILE)
protected:
    Int myGRank_;
    Int myGRoot_;
    //vector<int> Granks_;
public:
    void SetGlobalComm(const MPI_Comm & pGComm){
      if(commGlobRanks.count(comm_)==0){
        MPI_Group group2 = MPI_GROUP_NULL;
        MPI_Comm_group(pGComm, &group2);
        MPI_Group group1 = MPI_GROUP_NULL;
        MPI_Comm_group(comm_, &group1);

        Int size;
        MPI_Comm_size(comm_,&size);
        vector<int> globRanks(size);
        vector<int> Lranks(size);
        for(int i = 0; i<size;++i){Lranks[i]=i;}
        MPI_Group_translate_ranks(group1, size, &Lranks[0],group2, &globRanks[0]);
        commGlobRanks[comm_] = globRanks;
      }
      myGRoot_ = commGlobRanks[comm_][myRoot_];
      myGRank_ = commGlobRanks[comm_][myRank_];
     //   Granks_.resize(myDests_.size());
     //   for(int i = 0; i<myDests_.size();++i){
     //     Granks_[i] = globRanks[myDests_[i]];
     //   }
        //statusOFS<<myDests_<<std::endl;
        //statusOFS<<Granks_<<std::endl;
    }
#endif


  
    virtual void buildTree(Int * ranks, Int rank_cnt)=0;
  public:
    TreeBcast(){
      comm_ = MPI_COMM_WORLD;
      myRank_=-1;
      myRoot_ = -1; 
      msgSize_ = -1;
      numRecv_ = -1;
      tag_=-1;
      mainRoot_=-1;
      isReady_ = false;
    }


    TreeBcast(const MPI_Comm & pComm, Int * ranks, Int rank_cnt,Int msgSize){
      comm_ = pComm;
      MPI_Comm_rank(comm_,&myRank_);
      myRoot_ = -1; 
      msgSize_ = msgSize;

      numRecv_ = 0;
      tag_=-1;
      mainRoot_=ranks[0];
      isReady_ = false;
    }

    TreeBcast(const TreeBcast & Tree){
      Copy(Tree); 
    }
 
    virtual void Copy(const TreeBcast & Tree){
      comm_ = Tree.comm_;
      myRank_ = Tree.myRank_;
      myRoot_ = Tree.myRoot_; 
      msgSize_ = Tree.msgSize_;

      tag_= Tree.tag_;
      mainRoot_= Tree.mainRoot_;
      myDests_ = Tree.myDests_;

      //numRecv_ = Tree.numRecv_;
      //isReady_ = Tree.isReady_;
      isReady_ = false;
      numRecv_ = 0;
    }

    virtual TreeBcast * clone() const = 0; 

    void Reset(){
//statusOFS<<"RESET CALLED"<<std::endl;
      this->numRecv_ = 0;
      this->isReady_ = false;
    }


    
    static TreeBcast * Create(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize,double rseed);

    virtual inline Int GetNumRecvMsg(){return numRecv_;}
    virtual inline Int GetNumMsgToRecv(){return 1;}
 inline void SetDataReady(bool rdy){ isReady_=rdy; }
    inline void SetTag(Int tag){ tag_ = tag;}
      inline int GetTag(){ return tag_;}


    Int * GetDests(){ return &myDests_[0];}
    Int GetDest(Int i){ return myDests_[i];}
    Int GetDestCount(){ return myDests_.size();}
    Int GetRoot(){ return myRoot_;}
    Int GetMsgSize(){ return msgSize_;}

    void ForwardMessage( char * data, size_t size, int tag, MPI_Request * requests ){
        tag_ = tag;
                  for( Int idxRecv = 0; idxRecv < myDests_.size(); ++idxRecv ){
                    Int iProc = myDests_[idxRecv];
                    // Use Isend to send to multiple targets
                    MPI_Isend( data, size, MPI_BYTE, 
                        iProc, tag,comm_, &requests[2*iProc+1] );

#if defined(COMM_PROFILE_BCAST) || defined(COMM_PROFILE)
//        statusOFS<<idxRecv<<std::endl;
//        statusOFS<<Granks_<<std::endl;
      //PROFILE_COMM(myGRank_,Granks_[idxRecv],tag,msgSize_);
          PROFILE_COMM(myGRank_,commGlobRanks[comm_][iProc],tag_,msgSize_);
#endif
                  } // for (iProc)
    }


};

class FTreeBcast: public TreeBcast{
  protected:
    virtual void buildTree(Int * ranks, Int rank_cnt){

      Int idxStart = 0;
      Int idxEnd = rank_cnt;



      myRoot_ = ranks[0];

      if(myRank_==myRoot_){
        myDests_.insert(myDests_.begin(),&ranks[1],&ranks[0]+rank_cnt);
      }

#if (defined(BCAST_VERBOSE))
      statusOFS<<"My root is "<<myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<myDests_.size();++i){statusOFS<<myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }



  public:
    FTreeBcast(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeBcast(pComm,ranks,rank_cnt,msgSize){
      //build the binary tree;
      buildTree(ranks,rank_cnt);
    }


    virtual FTreeBcast * clone() const{
      FTreeBcast * out = new FTreeBcast(*this);
      return out;
    } 
};



class BTreeBcast: public TreeBcast{
  protected:
    virtual void buildTree(Int * ranks, Int rank_cnt){

      Int idxStart = 0;
      Int idxEnd = rank_cnt;



      Int prevRoot = ranks[0];
      while(idxStart<idxEnd){
        Int curRoot = ranks[idxStart];
        Int listSize = idxEnd - idxStart;

        if(listSize == 1){
          if(curRoot == myRank_){
            myRoot_ = prevRoot;
            break;
          }
        }
        else{
          Int halfList = floor(ceil(double(listSize) / 2.0));
          Int idxStartL = idxStart+1;
          Int idxStartH = idxStart+halfList;

          if(curRoot == myRank_){
            if ((idxEnd - idxStartH) > 0 && (idxStartH - idxStartL)>0){
              Int childL = ranks[idxStartL];
              Int childR = ranks[idxStartH];

              myDests_.push_back(childL);
              myDests_.push_back(childR);
            }
            else if ((idxEnd - idxStartH) > 0){
              Int childR = ranks[idxStartH];
              myDests_.push_back(childR);
            }
            else{
              Int childL = ranks[idxStartL];
              myDests_.push_back(childL);
            }
            myRoot_ = prevRoot;
            break;
          } 

          if( myRank_ < ranks[idxStartH]){
            idxStart = idxStartL;
            idxEnd = idxStartH;
          }
          else{
            idxStart = idxStartH;
          }
          prevRoot = curRoot;
        }

      }

#if (defined(BCAST_VERBOSE))
      statusOFS<<"My root is "<<myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<myDests_.size();++i){statusOFS<<myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }



  public:
    BTreeBcast(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeBcast(pComm,ranks,rank_cnt,msgSize){
      //build the binary tree;
      buildTree(ranks,rank_cnt);
    }

    virtual BTreeBcast * clone() const{
      BTreeBcast * out = new BTreeBcast(*this);
      return out;
    }

};



class ModBTreeBcast: public TreeBcast{
  protected:
    double rseed_;

    virtual void buildTree(Int * ranks, Int rank_cnt){

      Int idxStart = 0;
      Int idxEnd = rank_cnt;

      //sort the ranks with the modulo like operation
      if(rank_cnt>1){
        //Int new_idx = (int)((rand()+1.0) * (double)rank_cnt / ((double)RAND_MAX+1.0));

//        srand(ranks[0]+rank_cnt);
        //Int new_idx = (Int)rseed_ % (rank_cnt - 1) + 1; 
        Int new_idx = (int)((rank_cnt - 0) * ( (double)this->rseed_ / (double)RAND_MAX ) + 0);// (this->rseed_)%(rank_cnt-1)+1;
        //statusOFS<<"NEW IDX: "<<new_idx<<endl;



        Int * new_start = &ranks[new_idx];

        //for(int i =0;i<rank_cnt;++i){statusOFS<<ranks[i]<<" ";} statusOFS<<std::endl;

//        Int * new_start = std::lower_bound(&ranks[1],&ranks[0]+rank_cnt,ranks[0]);
        //just swap the two chunks   r[0] | r[1] --- r[new_start-1] | r[new_start] --- r[end]
        // becomes                   r[0] | r[new_start] --- r[end] | r[1] --- r[new_start-1] 
        std::rotate(&ranks[1], new_start, &ranks[0]+rank_cnt);

        //for(int i =0;i<rank_cnt;++i){statusOFS<<ranks[i]<<" ";} statusOFS<<std::endl;
      }

      Int prevRoot = ranks[0];
      while(idxStart<idxEnd){
        Int curRoot = ranks[idxStart];
        Int listSize = idxEnd - idxStart;

        if(listSize == 1){
          if(curRoot == myRank_){
            myRoot_ = prevRoot;
            break;
          }
        }
        else{
          Int halfList = floor(ceil(double(listSize) / 2.0));
          Int idxStartL = idxStart+1;
          Int idxStartH = idxStart+halfList;

          if(curRoot == myRank_){
            if ((idxEnd - idxStartH) > 0 && (idxStartH - idxStartL)>0){
              Int childL = ranks[idxStartL];
              Int childR = ranks[idxStartH];

              myDests_.push_back(childL);
              myDests_.push_back(childR);
            }
            else if ((idxEnd - idxStartH) > 0){
              Int childR = ranks[idxStartH];
              myDests_.push_back(childR);
            }
            else{
              Int childL = ranks[idxStartL];
              myDests_.push_back(childL);
            }
            myRoot_ = prevRoot;
            break;
          } 

          //not true anymore ?
          //first half to 
TIMER_START(FIND_RANK);
          Int * pos = std::find(&ranks[idxStartL], &ranks[idxStartH], myRank_);
TIMER_STOP(FIND_RANK);
          if( pos != &ranks[idxStartH]){
            idxStart = idxStartL;
            idxEnd = idxStartH;
          }
          else{
            idxStart = idxStartH;
          }
          prevRoot = curRoot;
        }

      }

#if (defined(REDUCE_VERBOSE))
      statusOFS<<"My root is "<<myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<myDests_.size();++i){statusOFS<<myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }



  public:
    ModBTreeBcast(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize, double rseed):TreeBcast(pComm,ranks,rank_cnt,msgSize){
      //build the binary tree;
      rseed_ = rseed;
      buildTree(ranks,rank_cnt);
    }

    virtual void Copy(const ModBTreeBcast & Tree){
      ((TreeBcast*)this)->Copy(*((const TreeBcast*)&Tree));


      rseed_ = Tree.rseed_;

    }
 
    virtual ModBTreeBcast * clone() const{
      ModBTreeBcast * out = new ModBTreeBcast(*this);
      return out;
    }

};


class RandBTreeBcast: public TreeBcast{
  protected:
    virtual void buildTree(Int * ranks, Int rank_cnt){

      Int idxStart = 0;
      Int idxEnd = rank_cnt;

      //random permute ranks
      if(rank_cnt>1){
        for(int i =0;i<rank_cnt;++i){statusOFS<<ranks[i]<<" ";} statusOFS<<std::endl;
        srand(ranks[0]);
        std::random_shuffle(&ranks[1],&ranks[0]+rank_cnt);
        for(int i =0;i<rank_cnt;++i){statusOFS<<ranks[i]<<" ";} statusOFS<<std::endl;

      }

      Int prevRoot = ranks[0];
      while(idxStart<idxEnd){
        Int curRoot = ranks[idxStart];
        Int listSize = idxEnd - idxStart;

        if(listSize == 1){
          if(curRoot == myRank_){
            myRoot_ = prevRoot;
            break;
          }
        }
        else{
          Int halfList = floor(ceil(double(listSize) / 2.0));
          Int idxStartL = idxStart+1;
          Int idxStartH = idxStart+halfList;

          if(curRoot == myRank_){
            if ((idxEnd - idxStartH) > 0 && (idxStartH - idxStartL)>0){
              Int childL = ranks[idxStartL];
              Int childR = ranks[idxStartH];

              myDests_.push_back(childL);
              myDests_.push_back(childR);
            }
            else if ((idxEnd - idxStartH) > 0){
              Int childR = ranks[idxStartH];
              myDests_.push_back(childR);
            }
            else{
              Int childL = ranks[idxStartL];
              myDests_.push_back(childL);
            }
            myRoot_ = prevRoot;
            break;
          } 

          //not true anymore ?
          //first half to 
          Int * pos = std::find(&ranks[idxStartL], &ranks[idxStartH], myRank_);
          if( pos != &ranks[idxStartH]){
            idxStart = idxStartL;
            idxEnd = idxStartH;
          }
          else{
            idxStart = idxStartH;
          }
          prevRoot = curRoot;
        }

      }

#if (defined(REDUCE_VERBOSE))
      statusOFS<<"My root is "<<myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<myDests_.size();++i){statusOFS<<myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }



  public:
    RandBTreeBcast(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeBcast(pComm,ranks,rank_cnt,msgSize){
      //build the binary tree;
      buildTree(ranks,rank_cnt);
    }

    virtual RandBTreeBcast * clone() const{
      RandBTreeBcast * out = new RandBTreeBcast(*this);
      return out;
    }

};




class PalmTreeBcast: public TreeBcast{
  protected:
  virtual void buildTree(Int * ranks, Int rank_cnt){
          Int numLevel = floor(log2(rank_cnt));
          Int numRoots = 0;
          for(Int level=0;level<numLevel;++level){
            numRoots = std::min( rank_cnt, numRoots + (Int)pow(2.0,level));
            Int numNextRoots = std::min(rank_cnt,numRoots + (Int)pow(2.0,(level+1)));
            Int numReceivers = numNextRoots - numRoots;
            for(Int ip = 0; ip<numRoots;++ip){
              Int p = ranks[ip];
              for(Int ir = ip; ir<numReceivers;ir+=numRoots){
                Int r = ranks[numRoots+ir];
                if(r==myRank_){
                  myRoot_ = p;
                }

                if(p==myRank_){
                  myDests_.push_back(r);
                }
              }
            }
          }

#if (defined(BCAST_VERBOSE))
      statusOFS<<"My root is "<<myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<myDests_.size();++i){statusOFS<<myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }



  public:
    PalmTreeBcast(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeBcast(pComm,ranks,rank_cnt,msgSize){
      //build the binary tree;
      buildTree(ranks,rank_cnt);
    }

    virtual PalmTreeBcast * clone() const{
      PalmTreeBcast * out = new PalmTreeBcast(*this);
      return out;
    }

};


template< typename T>
class TreeReduce: public TreeBcast{
  protected:
    T * myData_;
    MPI_Request sendRequest_;
    NumVec<char> myLocalBuffer_;
    NumVec<char> myRecvBuffers_;
    NumVec<T *> remoteData_;
    NumVec<MPI_Request> myRequests_;
    NumVec<MPI_Status> myStatuses_;
    NumVec<int> recvIdx_;

    bool fwded_;
    bool done_;
    bool isAllocated_;
    Int numRecvPosted_;

  public:
    TreeReduce(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeBcast(pComm,ranks,rank_cnt,msgSize){
      myData_ = NULL;
      sendRequest_ = MPI_REQUEST_NULL;
      fwded_=false;
      done_=false;
      isAllocated_=false;
      numRecvPosted_= 0;
    }


    virtual TreeReduce * clone() const = 0; 

    TreeReduce(const TreeReduce & Tree){
      this->Copy(Tree);
    }

    virtual void Copy(const TreeReduce & Tree){
      ((TreeBcast*)this)->Copy(*(const TreeBcast*)&Tree);

//      this->comm_ = Tree.comm_;
//      this->myRank_ = Tree.myRank_;
//      this->myRoot_ = Tree.myRoot_; 
//      this->msgSize_ = Tree.msgSize_;
//      this->numRecv_ = Tree.numRecv_;
//      this->tag_= Tree.tag_;
//      this->mainRoot_= Tree.mainRoot_;
//      this->isReady_ = Tree.isReady_;
//      this->myDests_ = Tree.myDests_;


      this->myData_ = NULL;
      this->sendRequest_ = MPI_REQUEST_NULL;
      this->fwded_= false;
      this->done_= false;
      this->isAllocated_= Tree.isAllocated_;
      this->numRecvPosted_= 0;

      //this->myLocalBuffer_.resize(Tree.myLocalBuffer_.size());
      //this->remoteData_ = Tree.remoteData_;
      //this->recvIdx_ = Tree.recvIdx_;

      CleanupBuffers();
    }
 


    bool IsAllocated(){return isAllocated_;}

    virtual ~TreeReduce(){
      CleanupBuffers();
    }


    static TreeReduce<T> * Create(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize,double rseed);

    virtual inline Int GetNumMsgToRecv(){return GetDestCount();}

    virtual void AllocRecvBuffers(){
      remoteData_.Resize(GetDestCount());
      //SetValue(remoteData_,(T*)NULL);

      //assert(myRecvBuffers_==NULL);
      //myRecvBuffers_ = new char[GetDestCount()*msgSize_];


      myRecvBuffers_.Resize(GetDestCount()*msgSize_);
      //SetValue(myRecvBuffers_,(char)0);

      for( Int idxRecv = 0; idxRecv < GetDestCount(); ++idxRecv ){
        remoteData_[idxRecv] = (T*)&(myRecvBuffers_[idxRecv*msgSize_]);
        //Int nelem = msgSize_ / sizeof(T);
        //std::fill(remoteData_[idxRecv],remoteData_[idxRecv]+nelem,ZERO<T>());
      }

      myRequests_.Resize(GetDestCount());
      SetValue(myRequests_,MPI_REQUEST_NULL);
      myStatuses_.Resize(GetDestCount());
      recvIdx_.Resize(GetDestCount());

      sendRequest_ = MPI_REQUEST_NULL;

      isAllocated_ = true;
    }

    void CleanupBuffers(){
      myLocalBuffer_.Clear();
//      if(myLocalBuffer_!=NULL){
//        delete []myLocalBuffer_;
//        myLocalBuffer_=NULL;
//      }


      remoteData_.Clear();
//      myRecvBuffers_.Clear();
//      if(myRecvBuffers_!=NULL){
//        delete []myRecvBuffers_;
//        myRecvBuffers_=NULL;
//      }


      myRequests_.Clear();
      myStatuses_.Clear();
      recvIdx_.Clear();


//              if(myLocalBuffer_!=NULL){
//                delete [] myLocalBuffer_;
//              }
//              myLocalBuffer_=NULL;


    }

      void Reset(){
//        assert(done_ || myDests_.size()==0);
        CleanupBuffers();
        done_=false;

      myData_ = NULL;
      sendRequest_ = MPI_REQUEST_NULL;
      fwded_=false;
      isAllocated_=false;
      isReady_=false;
      numRecv_ = 0;
      numRecvPosted_= 0;
      }



    void SetLocalBuffer(T * locBuffer){
      if(myData_!=NULL && myData_!=locBuffer){

//statusOFS<<"DOING SUM"<<std::endl;
//gdb_lock();
        blas::Axpy(msgSize_/sizeof(T), ONE<T>(), myData_, 1, locBuffer, 1 );
        myLocalBuffer_.Clear(); 
      }

      myData_ = locBuffer;
    }

    inline bool AccumulationDone(){
      if(myRank_==myRoot_ && isAllocated_){
        isReady_=true;
      }
      return isReady_ && (numRecv_ == GetDestCount());
    }


    inline bool IsDone(){
      if(myRank_==myRoot_ && isAllocated_){
        isReady_=true;
      }

      bool retVal = AccumulationDone();
      if(myRoot_ != myRank_ && !fwded_){
        retVal = false;
      }

      if (retVal && myRoot_ != myRank_ && fwded_){
        //test the send request
        int flag = 0;
        MPI_Test(&sendRequest_,&flag,MPI_STATUS_IGNORE);
        retVal = flag==1;
      }
      
      return retVal;
    }

    //async wait and forward
    virtual bool Progress(){

        if(done_){
          return true;
        }
      if(!isAllocated_){
        return false;
      }

      if(myRank_==myRoot_ && isAllocated_){
        isReady_=true;
      }

//      if(this->numRecvPosted_==0){ 
//        this->PostFirstRecv();
//      }

      bool retVal = AccumulationDone();
      if(isReady_ && !retVal){

        //assert(isAllocated_);

        //mpi_test_some on my requests
        int recvCount = -1;
        int reqCnt = GetDestCount();
        assert(reqCnt == myRequests_.m());
        MPI_Testsome(reqCnt,&myRequests_[0],&recvCount,&recvIdx_[0],&myStatuses_[0]);
        //if something has been received, accumulate and potentially forward it
        for(Int i = 0;i<recvCount;++i ){
          Int idx = recvIdx_[i];

         if(idx!=MPI_UNDEFINED){

          Int size = 0;
          MPI_Get_count(&myStatuses_[i], MPI_BYTE, &size);


#if ( _DEBUGlevel_ >= 1 ) || defined(REDUCE_VERBOSE)
        statusOFS<<myRank_<<" RECVD from "<<myStatuses_[i].MPI_SOURCE<<" on tag "<<tag_<<std::endl;
#endif
          if(size>0){
            //If myData is 0, allocate to the size of what has been received
            if(myData_==NULL){
              //assert(size==msgSize_);
              myLocalBuffer_.Resize(msgSize_);
    
              myData_ = (T*)&myLocalBuffer_[0];
              Int nelem = +msgSize_/sizeof(T);
              std::fill(myData_,myData_+nelem,ZERO<T>());
            }

            Reduce(idx,i);

          }

          numRecv_++;
          //MPI_Request_free(&myRequests_[idx]);
          }
        }

      }
      else if (isReady_ && sendRequest_ == MPI_REQUEST_NULL && myRoot_ != myRank_ && !fwded_){
        //free the unnecessary arrays
        myRecvBuffers_.Clear();
        myRequests_.Clear();
        myStatuses_.Clear();
        recvIdx_.Clear();

        //assert(isAllocated_);

        //Forward
        Forward();
        retVal = false;
      }
      else{
        retVal = IsDone();
        if(retVal){
          //free the unnecessary arrays
          myRecvBuffers_.Clear();
          myRequests_.Clear();
          myStatuses_.Clear();
          recvIdx_.Clear();
        }
      }

      
      done_ = retVal;
      return retVal;
    }

    //blocking wait
    void Wait(){
      if(!done_){
        while(!Progress());
      }
    }

    T * GetLocalBuffer(){
       return myData_;
    }



    void CopyLocalBuffer(T* destBuffer){
       std::copy((char*)myData_,(char*)myData_+GetMsgSize(),(char*)destBuffer);
    }


    virtual void PostFirstRecv()
    {
      if(this->GetDestCount()>this->numRecvPosted_){
        for( Int idxRecv = 0; idxRecv < myDests_.size(); ++idxRecv ){
          Int iProc = myDests_[idxRecv];
          //assert(msgSize_>=0);
          MPI_Irecv( (char*)remoteData_[idxRecv], msgSize_, MPI_BYTE, 
              iProc, tag_,comm_, &myRequests_[idxRecv] );
          this->numRecvPosted_++;
        } // for (iProc)
      }
    }




    protected:
    virtual void Reduce( Int idxRecv, Int idReq){
      //add thing to my data
      blas::Axpy(msgSize_/sizeof(T), ONE<T>(), remoteData_[idxRecv], 1, myData_, 1 );
    }

    void Forward(){ 
      //forward to my root if I have reseived everything
      Int iProc = myRoot_;
      // Use Isend to send to multiple targets
      if(myData_==NULL){
        MPI_Isend( NULL, 0, MPI_BYTE, 
            iProc, tag_,comm_, &sendRequest_ );
#ifdef COMM_PROFILE
      PROFILE_COMM(myGRank_,myGRoot_,tag_,0);
#endif
      }
      else{
        MPI_Isend( (char*)myData_, msgSize_, MPI_BYTE, 
            iProc, tag_,comm_, &sendRequest_ );
#ifdef COMM_PROFILE
      PROFILE_COMM(myGRank_,myGRoot_,tag_,msgSize_);
#endif
      }

#if ( _DEBUGlevel_ >= 1 ) || defined(REDUCE_VERBOSE)
        statusOFS<<myRank_<<" FWD to "<<iProc<<" on tag "<<tag_<<std::endl;
#endif

      fwded_ = true;
      
    }

};


template< typename T>
class FTreeReduce: public TreeReduce<T>{
    protected:
    virtual void buildTree(Int * ranks, Int rank_cnt){

      Int idxStart = 0;
      Int idxEnd = rank_cnt;



      this->myRoot_ = ranks[0];

      if(this->myRank_==this->myRoot_){
        this->myDests_.insert(this->myDests_.begin(),&ranks[1],&ranks[0]+rank_cnt);
      }

#if (defined(REDUCE_VERBOSE))
      statusOFS<<"My root is "<<this->myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<this->myDests_.size();++i){statusOFS<<this->myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }

    virtual void Reduce( ){
      //add thing to my data
      blas::Axpy(this->msgSize_/sizeof(T), ONE<T>(), this->remoteData_[0], 1, this->myData_, 1 );


#if (defined(REDUCE_DEBUG))
     statusOFS << std::endl << /*"["<<snode.Index<<"]*/" Recv contrib"<< std::endl; 
      for(int i = 0; i < this->msgSize_/sizeof(T); ++i){
        statusOFS<< this->remoteData_[0][i]<< " ";
        if(i%3==0){statusOFS<<std::endl;}
      }
      statusOFS<<std::endl;

     statusOFS << std::endl << /*"["<<snode.Index<<"]*/" Reduce buffer now is"<< std::endl; 
      for(int i = 0; i < this->msgSize_/sizeof(T); ++i){
        statusOFS<< this->myData_[i]<< " ";
        if(i%3==0){statusOFS<<std::endl;}
      }
      statusOFS<<std::endl;
#endif

    }



    public:
    FTreeReduce(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeReduce<T>(pComm, ranks, rank_cnt, msgSize){
      buildTree(ranks,rank_cnt);
    }

    virtual void PostFirstRecv()
    {
//        if(!this->isAllocated_){
//          this->AllocRecvBuffers();
//        }
        if(this->isAllocated_ && this->GetDestCount()>this->numRecvPosted_){
          MPI_Irecv( (char*)this->remoteData_[0], this->msgSize_, MPI_BYTE, 
            MPI_ANY_SOURCE, this->tag_,this->comm_, &this->myRequests_[0] );
          this->numRecvPosted_++;
        }
    }

    virtual void AllocRecvBuffers(){
      if(this->GetDestCount()>0){
        this->remoteData_.Resize(1);

        this->myRecvBuffers_.Resize(this->msgSize_);

        this->remoteData_[0] = (T*)&(this->myRecvBuffers_[0]);

        this->myRequests_.Resize(1);
        SetValue(this->myRequests_,MPI_REQUEST_NULL);
        this->myStatuses_.Resize(1);
        this->recvIdx_.Resize(1);
      }

      this->sendRequest_ = MPI_REQUEST_NULL;

      this->isAllocated_ = true;
    }

    virtual bool Progress(){

      if(!this->isAllocated_){
        return true;
      }


      if(this->myRank_==this->myRoot_ && this->isAllocated_){
        this->isReady_=true;
      }
     
//      if(this->numRecvPosted_==0){ 
//        this->PostFirstRecv();
//      }

      bool retVal = this->AccumulationDone();
      if(this->isReady_ && !retVal){

        //assert(this->isAllocated_);

        //mpi_test_some on my requests
        int recvCount = -1;
        int reqCnt = 1;

        MPI_Testsome(reqCnt,&this->myRequests_[0],&recvCount,&this->recvIdx_[0],&this->myStatuses_[0]);
        //MPI_Waitsome(reqCnt,&myRequests_[0],&recvCount,&recvIdx_[0],&myStatuses_[0]);
        //if something has been received, accumulate and potentially forward it
        for(Int i = 0;i<recvCount;++i ){
          Int idx = this->recvIdx_[i];

         if(idx!=MPI_UNDEFINED){

          Int size = 0;
          MPI_Get_count(&this->myStatuses_[i], MPI_BYTE, &size);


#if ( _DEBUGlevel_ >= 1 ) || defined(REDUCE_VERBOSE)

        statusOFS<<this->myRank_<<" RECVD from "<<this->myStatuses_[i].MPI_SOURCE<<" on tag "<<this->tag_<<std::endl;
#endif
          if(size>0){
            //If myData is 0, allocate to the size of what has been received
            if(this->myData_==NULL){
              //assert(size==this->msgSize_);
              this->myLocalBuffer_.Resize(this->msgSize_);
    
              this->myData_ = (T*)&this->myLocalBuffer_[0];
              Int nelem = this->msgSize_/sizeof(T);
              std::fill(this->myData_,this->myData_+nelem,ZERO<T>());
            }

            this->Reduce();
          }

          this->numRecv_++;
          }
        }

        if(recvCount>0){
          this->PostFirstRecv();
        }
      }
      else if (this->isReady_ && this->sendRequest_ == MPI_REQUEST_NULL && this->myRoot_ != this->myRank_ && !this->fwded_){
        //free the unnecessary arrays
        this->myRecvBuffers_.Clear();
        this->myRequests_.Clear();
        this->myStatuses_.Clear();
        this->recvIdx_.Clear();

        //Forward
        this->Forward();
        retVal = false;
      }
      else{
        retVal = this->IsDone();
        if(retVal){
        //free the unnecessary arrays
        this->myRecvBuffers_.Clear();
        this->myRequests_.Clear();
        this->myStatuses_.Clear();
        this->recvIdx_.Clear();
        }
      }

      return retVal;
    }


    virtual FTreeReduce * clone() const{
      FTreeReduce * out = new FTreeReduce(*this);
      return out;
    }



};



template< typename T>
class BTreeReduce: public TreeReduce<T>{
    protected:
    virtual void buildTree(Int * ranks, Int rank_cnt){
      Int idxStart = 0;
      Int idxEnd = rank_cnt;



      Int prevRoot = ranks[0];
      while(idxStart<idxEnd){
        Int curRoot = ranks[idxStart];
        Int listSize = idxEnd - idxStart;

        if(listSize == 1){
          if(curRoot == this->myRank_){
            this->myRoot_ = prevRoot;
            break;
          }
        }
        else{
          Int halfList = floor(ceil(double(listSize) / 2.0));
          Int idxStartL = idxStart+1;
          Int idxStartH = idxStart+halfList;

          if(curRoot == this->myRank_){
            if ((idxEnd - idxStartH) > 0 && (idxStartH - idxStartL)>0){
              Int childL = ranks[idxStartL];
              Int childR = ranks[idxStartH];

              this->myDests_.push_back(childL);
              this->myDests_.push_back(childR);
            }
            else if ((idxEnd - idxStartH) > 0){
              Int childR = ranks[idxStartH];
              this->myDests_.push_back(childR);
            }
            else{
              Int childL = ranks[idxStartL];
              this->myDests_.push_back(childL);
            }
            this->myRoot_ = prevRoot;
            break;
          } 

          if( this->myRank_ < ranks[idxStartH]){
            idxStart = idxStartL;
            idxEnd = idxStartH;
          }
          else{
            idxStart = idxStartH;
          }
          prevRoot = curRoot;
        }

      }

#if (defined(REDUCE_VERBOSE))
      statusOFS<<"My root is "<<this->myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<this->myDests_.size();++i){statusOFS<<this->myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }
    public:
    BTreeReduce(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeReduce<T>(pComm, ranks, rank_cnt, msgSize){
      buildTree(ranks,rank_cnt);
    }

    virtual BTreeReduce * clone() const{
      BTreeReduce * out = new BTreeReduce(*this);
      return out;
    }
};


template< typename T>
class ModBTreeReduce: public TreeReduce<T>{
    protected:
    double rseed_;
    virtual void buildTree(Int * ranks, Int rank_cnt){

      Int idxStart = 0;
      Int idxEnd = rank_cnt;

      //sort the ranks with the modulo like operation
      if(rank_cnt>1){
        //generate a random position in [1 .. rand_cnt]
        //Int new_idx = (int)((rand()+1.0) * (double)rank_cnt / ((double)RAND_MAX+1.0));
        //srand(ranks[0]+rank_cnt);
        //Int new_idx = rseed_%(rank_cnt-1)+1;

        //Int new_idx = (int)((rank_cnt - 0) * ( (double)this->rseed_ / (double)RAND_MAX ) + 0);// (this->rseed_)%(rank_cnt-1)+1;
        //Int new_idx = (Int)rseed_ % (rank_cnt - 1) + 1; 
        Int new_idx = (int)((rank_cnt - 0) * ( (double)this->rseed_ / (double)RAND_MAX ) + 0);// (this->rseed_)%(rank_cnt-1)+1;
        //Int new_idx = (int)((rank_cnt - 0) * ( (double)this->rseed_ / (double)RAND_MAX ) + 0);// (this->rseed_)%(rank_cnt-1)+1;


        Int * new_start = &ranks[new_idx];
//        for(int i =0;i<rank_cnt;++i){statusOFS<<ranks[i]<<" ";} statusOFS<<std::endl;
        
//        Int * new_start = std::lower_bound(&ranks[1],&ranks[0]+rank_cnt,ranks[0]);
        //just swap the two chunks   r[0] | r[1] --- r[new_start-1] | r[new_start] --- r[end]
        // becomes                   r[0] | r[new_start] --- r[end] | r[1] --- r[new_start-1] 
        std::rotate(&ranks[1], new_start, &ranks[0]+rank_cnt);
//        for(int i =0;i<rank_cnt;++i){statusOFS<<ranks[i]<<" ";} statusOFS<<std::endl;
      }

      Int prevRoot = ranks[0];
      while(idxStart<idxEnd){
        Int curRoot = ranks[idxStart];
        Int listSize = idxEnd - idxStart;

        if(listSize == 1){
          if(curRoot == this->myRank_){
            this->myRoot_ = prevRoot;
            break;
          }
        }
        else{
          Int halfList = floor(ceil(double(listSize) / 2.0));
          Int idxStartL = idxStart+1;
          Int idxStartH = idxStart+halfList;

          if(curRoot == this->myRank_){
            if ((idxEnd - idxStartH) > 0 && (idxStartH - idxStartL)>0){
              Int childL = ranks[idxStartL];
              Int childR = ranks[idxStartH];

              this->myDests_.push_back(childL);
              this->myDests_.push_back(childR);
            }
            else if ((idxEnd - idxStartH) > 0){
              Int childR = ranks[idxStartH];
              this->myDests_.push_back(childR);
            }
            else{
              Int childL = ranks[idxStartL];
              this->myDests_.push_back(childL);
            }
            this->myRoot_ = prevRoot;
            break;
          } 

          //not true anymore ?
          //first half to 
TIMER_START(FIND_RANK);
          Int * pos = std::find(&ranks[idxStartL], &ranks[idxStartH], this->myRank_);
TIMER_STOP(FIND_RANK);
          if( pos != &ranks[idxStartH]){
            idxStart = idxStartL;
            idxEnd = idxStartH;
          }
          else{
            idxStart = idxStartH;
          }
          prevRoot = curRoot;
        }

      }

#if (defined(REDUCE_VERBOSE)) 
      statusOFS<<"My root is "<<this->myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<this->myDests_.size();++i){statusOFS<<this->myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }
    public:
    ModBTreeReduce(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize, double rseed):TreeReduce<T>(pComm, ranks, rank_cnt, msgSize){
      this->rseed_ = rseed;
      buildTree(ranks,rank_cnt);
    }

    virtual void Copy(const ModBTreeReduce & Tree){
      ((TreeReduce<T>*)this)->Copy(*((const TreeReduce<T>*)&Tree));
      //this->comm_ = Tree.comm_;
      //this->myRank_ = Tree.myRank_;
      //this->myRoot_ = Tree.myRoot_; 
      //this->msgSize_ = Tree.msgSize_;

      //this->numRecv_ = Tree.numRecv_;
      //this->tag_= Tree.tag_;
      //this->mainRoot_= Tree.mainRoot_;
      //this->isReady_ = Tree.isReady_;
      //this->myDests_ = Tree.myDests_;


      //this->myData_ = Tree.myData_;
      //this->sendRequest_ = Tree.sendRequest_;
      //this->fwded_= Tree.fwded_;
      //this->isAllocated_= Tree.isAllocated_;
      //this->numRecvPosted_= Tree.numRecvPosted_;

      //this->myLocalBuffer_ = Tree.myLocalBuffer_;
      //this->myRecvBuffers_ = Tree.myRecvBuffers_;
      //this->remoteData_ = Tree.remoteData_;
      //this->myRequests_ = Tree.myRequests_;
      //this->myStatuses_ = Tree.myStatuses_;
      //this->recvIdx_ = Tree.recvIdx_;
      this->rseed_ = Tree.rseed_;
    }
 



    virtual ModBTreeReduce * clone() const{
      ModBTreeReduce * out = new ModBTreeReduce(*this);
      return out;
    }

};


template< typename T>
class PalmTreeReduce: public TreeReduce<T>{
    protected:

  virtual void buildTree(Int * ranks, Int rank_cnt){
          Int numLevel = floor(log2(rank_cnt));
          Int numRoots = 0;
          for(Int level=0;level<numLevel;++level){
            numRoots = std::min( rank_cnt, numRoots + (Int)pow(2,level));
            Int numNextRoots = std::min(rank_cnt,numRoots + (Int)pow(2,(level+1)));
            Int numReceivers = numNextRoots - numRoots;
            for(Int ip = 0; ip<numRoots;++ip){
              Int p = ranks[ip];
              for(Int ir = ip; ir<numReceivers;ir+=numRoots){
                Int r = ranks[numRoots+ir];
                if(r==this->myRank_){
                  this->myRoot_ = p;
                }

                if(p==this->myRank_){
                  this->myDests_.push_back(r);
                }
              }
            }
          }

#if (defined(BCAST_VERBOSE))
      statusOFS<<"My root is "<<this->myRoot_<<std::endl;
      statusOFS<<"My dests are ";
      for(int i =0;i<this->myDests_.size();++i){statusOFS<<this->myDests_[i]<<" ";}
      statusOFS<<std::endl;
#endif
    }



  public:
    PalmTreeReduce(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize):TreeReduce<T>(pComm,ranks,rank_cnt,msgSize){
      //build the binary tree;
      buildTree(ranks,rank_cnt);
    }



    virtual void Copy(const PalmTreeReduce & Tree){
      ((TreeReduce<T>*)this)->Copy(*((const TreeReduce<T>*)&Tree));
      //this->comm_ = Tree.comm_;
      //this->myRank_ = Tree.myRank_;
      //this->myRoot_ = Tree.myRoot_; 
      //this->msgSize_ = Tree.msgSize_;

      //this->numRecv_ = Tree.numRecv_;
      //this->tag_= Tree.tag_;
      //this->mainRoot_= Tree.mainRoot_;
      //this->isReady_ = Tree.isReady_;
      //this->myDests_ = Tree.myDests_;


      //this->myData_ = Tree.myData_;
      //this->sendRequest_ = Tree.sendRequest_;
      //this->fwded_= Tree.fwded_;
      //this->isAllocated_= Tree.isAllocated_;
      //this->numRecvPosted_= Tree.numRecvPosted_;

      //this->myLocalBuffer_ = Tree.myLocalBuffer_;
      //this->myRecvBuffers_ = Tree.myRecvBuffers_;
      //this->remoteData_ = Tree.remoteData_;
      //this->myRequests_ = Tree.myRequests_;
      //this->myStatuses_ = Tree.myStatuses_;
      //this->recvIdx_ = Tree.recvIdx_;
      //this->rseed_ = Tree.rseed_;
    }
 



    virtual PalmTreeReduce * clone() const{
      PalmTreeReduce * out = new PalmTreeReduce(*this);
      return out;
    }

};



    inline TreeBcast * TreeBcast::Create(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize, double rseed){
      //get communicator size
      Int nprocs = 0;
      MPI_Comm_size(pComm, &nprocs);


#if defined(FTREE)
        return new FTreeBcast(pComm,ranks,rank_cnt,msgSize);
#elif defined(MODBTREE)
        return new ModBTreeBcast(pComm,ranks,rank_cnt,msgSize, rseed);
#elif defined(BTREE)
        return new BTreeBcast(pComm,ranks,rank_cnt,msgSize);
#elif defined(PALMTREE)
        return new PalmTreeBcast(pComm,ranks,rank_cnt,msgSize);
#endif


//      return new PalmTreeBcast(pComm,ranks,rank_cnt,msgSize);
//      return new ModBTreeBcast(pComm,ranks,rank_cnt,msgSize, rseed);
//      return new RandBTreeBcast(pComm,ranks,rank_cnt,msgSize);

      if(nprocs<=FTREE_LIMIT){
        return new FTreeBcast(pComm,ranks,rank_cnt,msgSize);
      }
      else{
        return new ModBTreeBcast(pComm,ranks,rank_cnt,msgSize, rseed);
      }




    }




template< typename T>
    inline TreeReduce<T> * TreeReduce<T>::Create(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize, double rseed){
      //get communicator size
      Int nprocs = 0;
      MPI_Comm_size(pComm, &nprocs);

#if defined(FTREE)
        return new FTreeReduce<T>(pComm,ranks,rank_cnt,msgSize);
#elif defined(MODBTREE)
        return new ModBTreeReduce<T>(pComm,ranks,rank_cnt,msgSize, rseed);
#elif defined(BTREE)
        return new BTreeReduce<T>(pComm,ranks,rank_cnt,msgSize);
#elif defined(PALMTREE)
        return new PalmTreeReduce<T>(pComm,ranks,rank_cnt,msgSize);
#endif


      if(nprocs<=FTREE_LIMIT){
#if ( _DEBUGlevel_ >= 1 ) || defined(REDUCE_VERBOSE)
statusOFS<<"FLAT TREE USED"<<endl;
#endif
        return new FTreeReduce<T>(pComm,ranks,rank_cnt,msgSize);
      }
      else{
#if ( _DEBUGlevel_ >= 1 ) || defined(REDUCE_VERBOSE)
statusOFS<<"BINARY TREE USED"<<endl;
#endif
        return new ModBTreeReduce<T>(pComm,ranks,rank_cnt,msgSize, rseed);
        //return new BTreeReduce<T>(pComm,ranks,rank_cnt,msgSize);
      }
    }


#ifdef NEW_BCAST
template< typename T>
    inline TreeBcast2<T> * TreeBcast2<T>::Create(const MPI_Comm & pComm, Int * ranks, Int rank_cnt, Int msgSize, double rseed){
      //get communicator size
      Int nprocs = 0;
      MPI_Comm_size(pComm, &nprocs);




#if defined(FTREE)
        return new FTreeBcast2<T>(pComm,ranks,rank_cnt,msgSize);
#elif defined(MODBTREE)
        return new ModBTreeBcast2<T>(pComm,ranks,rank_cnt,msgSize,rseed);
#elif defined(BTREE)
        return new BTreeBcast2<T>(pComm,ranks,rank_cnt,msgSize);
#endif


      if(nprocs<=FTREE_LIMIT){
#if ( _DEBUGlevel_ >= 1 ) || defined(REDUCE_VERBOSE)
statusOFS<<"FLAT TREE USED"<<endl;
#endif

        return new FTreeBcast2<T>(pComm,ranks,rank_cnt,msgSize);

      }
      else{
#if ( _DEBUGlevel_ >= 1 ) || defined(REDUCE_VERBOSE)
statusOFS<<"BINARY TREE USED"<<endl;
#endif
        return new ModBTreeBcast2<T>(pComm,ranks,rank_cnt,msgSize, rseed);
//        //return new BTreeReduce<T>(pComm,ranks,rank_cnt,msgSize);
      }
    }
#endif






}

#endif