From: Anthony Geay <anthony.geay@edf.fr>
Date: Wed, 29 Apr 2020 09:25:15 +0000 (+0200)
Subject: WIP
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=e48540f06dae28dbc59aa3be9de1a514267275db;p=tools%2Fmedcoupling.git

WIP
---

diff --git a/src/ParaMEDMEM/CommInterface.cxx b/src/ParaMEDMEM/CommInterface.cxx
index 6fda70ddd..b61c05bd7 100644
--- a/src/ParaMEDMEM/CommInterface.cxx
+++ b/src/ParaMEDMEM/CommInterface.cxx
@@ -62,22 +62,18 @@ namespace MEDCoupling
     \endverbatim
   */
 
+  void CommInterface::gatherArrays(MPI_Comm comm, int root, const DataArrayIdType *array, std::vector< MCAuto<DataArrayIdType> >& arraysOut) const
+  {
+    this->gatherArraysT2<mcIdType>(comm,root,array,arraysOut);
+  }
+
   /*!
    * Generalized AllGather collective communication.
    * This method send input \a array to all procs.
    */
   void CommInterface::allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::vector< MCAuto<DataArrayIdType> >& arraysOut) const
   {
-    std::unique_ptr<mcIdType[]> result, resultIndex;
-    int size(this->allGatherArrays(comm,array,result,resultIndex));
-    arraysOut.resize(size);
-    for(int i = 0 ; i < size ; ++i)
-    {
-      arraysOut[i] = DataArrayIdType::New();
-      mcIdType nbOfEltPack(resultIndex[i+1]-resultIndex[i]);
-      arraysOut[i]->alloc(nbOfEltPack,1);
-      std::copy(result.get()+resultIndex[i],result.get()+resultIndex[i+1],arraysOut[i]->getPointer());
-    }
+    this->allGatherArraysT2<mcIdType>(comm,array,arraysOut);
   }
 
   /*!
@@ -86,18 +82,7 @@ namespace MEDCoupling
    */
   int CommInterface::allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::unique_ptr<mcIdType[]>& result, std::unique_ptr<mcIdType[]>& resultIndex) const
   {
-    int size;
-    this->commSize(comm,&size);
-    std::unique_ptr<mcIdType[]> nbOfElems(new mcIdType[size]);
-    mcIdType nbOfCellsRequested(array->getNumberOfTuples());
-    this->allGather(&nbOfCellsRequested,1,MPI_ID_TYPE,nbOfElems.get(),1,MPI_ID_TYPE,comm);
-    mcIdType nbOfCellIdsSum(std::accumulate(nbOfElems.get(),nbOfElems.get()+size,0));
-    result.reset(new mcIdType[nbOfCellIdsSum]);
-    std::unique_ptr<int[]> nbOfElemsInt( CommInterface::ToIntArray<mcIdType>(nbOfElems,size) );
-    std::unique_ptr<int[]> offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) );
-    this->allGatherV(array->begin(),nbOfCellsRequested,MPI_ID_TYPE,result.get(),nbOfElemsInt.get(),offsetsIn.get(),MPI_ID_TYPE,comm);
-    resultIndex = ComputeOffsetFull<mcIdType>(nbOfElems,size);
-    return size;
+    return this->allGatherArraysT<mcIdType>(comm,array,result,resultIndex);
   }
 
   void CommInterface::allToAllArrays(MPI_Comm comm, const std::vector< MCAuto<DataArrayDouble> >& arrays, std::vector< MCAuto<DataArrayDouble> >& arraysOut) const
diff --git a/src/ParaMEDMEM/CommInterface.hxx b/src/ParaMEDMEM/CommInterface.hxx
index e969f9394..a59fce0b9 100644
--- a/src/ParaMEDMEM/CommInterface.hxx
+++ b/src/ParaMEDMEM/CommInterface.hxx
@@ -101,6 +101,8 @@ namespace MEDCoupling
     int getCount(MPI_Status *status, MPI_Datatype datatype, int *count) const { return MPI_Get_count(status, datatype, count); }
 
     int broadcast(void* buffer, int count, MPI_Datatype datatype, int root, MPI_Comm comm) const { return MPI_Bcast(buffer, count,  datatype, root, comm); }
+    int gather(const void* sendbuf, int sendcount, MPI_Datatype sendtype, void* recvbuf, int recvcount, MPI_Datatype recvtype, int root, MPI_Comm comm) const { return MPI_Gather(sendbuf,sendcount,sendtype,recvbuf,recvcount,recvtype,root,comm); }
+    int gatherV(const void* sendbuf, int sendcount, MPI_Datatype sendtype, void* recvbuf, const int recvcounts[], const int displs[], MPI_Datatype recvtype, int root, MPI_Comm comm) const { return MPI_Gatherv(sendbuf,sendcount,sendtype,recvbuf,recvcounts,displs,recvtype,root,comm); }
     int allGather(void* sendbuf, int sendcount, MPI_Datatype sendtype,
                   void* recvbuf, int recvcount, MPI_Datatype recvtype,
                   MPI_Comm comm) const { return MPI_Allgather(sendbuf,sendcount, sendtype, recvbuf, recvcount, recvtype, comm); }
@@ -118,11 +120,96 @@ namespace MEDCoupling
                MPI_Op op, int root, MPI_Comm comm) const { return MPI_Reduce(sendbuf, recvbuf, count, datatype, op, root, comm); }
     int allReduce(void* sendbuf, void* recvbuf, int count, MPI_Datatype datatype, MPI_Op op, MPI_Comm comm) const { return MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm); }
   public:
+    void gatherArrays(MPI_Comm comm, int root, const DataArrayIdType *array, std::vector< MCAuto<DataArrayIdType> >& arraysOut) const;
     void allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::vector< MCAuto<DataArrayIdType> >& arraysOut) const;
     int allGatherArrays(MPI_Comm comm, const DataArrayIdType *array, std::unique_ptr<mcIdType[]>& result, std::unique_ptr<mcIdType[]>& resultIndex) const;
     void allToAllArrays(MPI_Comm comm, const std::vector< MCAuto<DataArrayIdType> >& arrays, std::vector< MCAuto<DataArrayIdType> >& arraysOut) const;
     void allToAllArrays(MPI_Comm comm, const std::vector< MCAuto<DataArrayDouble> >& arrays, std::vector< MCAuto<DataArrayDouble> >& arraysOut) const;
     void allToAllArrays(MPI_Comm comm, const std::vector< MCAuto<DataArrayDouble> >& arrays, MCAuto<DataArrayDouble>& arraysOut) const;
+    
+    template<class T>
+    int gatherArraysT(MPI_Comm comm, int root, const typename Traits<T>::ArrayType *array, std::unique_ptr<T[]>& result, std::unique_ptr<mcIdType[]>& resultIndex, int& rank) const
+    {
+      using DataArrayT = typename Traits<T>::ArrayType;
+      int size;
+      this->commSize(comm,&size);
+      rank = -1;
+      this->commRank(comm,&rank);
+      std::unique_ptr<mcIdType[]> nbOfElems;
+      if(rank==0)
+        nbOfElems.reset(new mcIdType[size]);
+      mcIdType nbOfCellsRequested(array->getNumberOfTuples());
+      this->gather(&nbOfCellsRequested,1,MPI_ID_TYPE,nbOfElems.get(),1,MPI_ID_TYPE,root,comm);
+      mcIdType nbOfCellIdsSum(std::accumulate(nbOfElems.get(),nbOfElems.get()+size,0));
+      if(rank==0)
+        result.reset(new T[nbOfCellIdsSum]);
+      std::unique_ptr<int[]> nbOfElemsInt,offsetsIn;
+      if(rank==0)
+      {
+        nbOfElemsInt = CommInterface::ToIntArray<mcIdType>(nbOfElems,size);
+        offsetsIn = CommInterface::ComputeOffset(nbOfElemsInt,size);
+      }
+      this->gatherV(array->begin(),nbOfCellsRequested,ParaTraits<T>::MPIDataType,result.get(),nbOfElemsInt.get(),offsetsIn.get(),ParaTraits<T>::MPIDataType,root,comm);
+      resultIndex = ComputeOffsetFull<mcIdType>(nbOfElems,size);
+      return size;
+    }
+
+    template<class T>
+    void gatherArraysT2(MPI_Comm comm, int root, const typename Traits<T>::ArrayType *array, std::vector< MCAuto<typename Traits<T>::ArrayType> >& arraysOut) const
+    {
+      using DataArrayT = typename Traits<T>::ArrayType;
+      std::unique_ptr<T[]> result;
+      std::unique_ptr<mcIdType[]> resultIndex;
+      int rank(-1);
+      int size(this->gatherArraysT<T>(comm,root,array,result,resultIndex,rank));
+      arraysOut.resize(size);
+      for(int i = 0 ; i < size ; ++i)
+      {
+        arraysOut[i] = DataArrayT::New();
+        if(rank == 0)
+        {
+          mcIdType nbOfEltPack(resultIndex[i+1]-resultIndex[i]);
+          arraysOut[i]->alloc(nbOfEltPack,1);
+          std::copy(result.get()+resultIndex[i],result.get()+resultIndex[i+1],arraysOut[i]->getPointer());
+        }
+      }
+    }
+
+    template<class T>
+    int allGatherArraysT(MPI_Comm comm, const typename Traits<T>::ArrayType *array, std::unique_ptr<T[]>& result, std::unique_ptr<mcIdType[]>& resultIndex) const
+    {
+      using DataArrayT = typename Traits<T>::ArrayType;
+      int size;
+      this->commSize(comm,&size);
+      std::unique_ptr<mcIdType[]> nbOfElems(new mcIdType[size]);
+      mcIdType nbOfCellsRequested(array->getNumberOfTuples());
+      this->allGather(&nbOfCellsRequested,1,MPI_ID_TYPE,nbOfElems.get(),1,MPI_ID_TYPE,comm);
+      mcIdType nbOfCellIdsSum(std::accumulate(nbOfElems.get(),nbOfElems.get()+size,0));
+      result.reset(new T[nbOfCellIdsSum]);
+      std::unique_ptr<int[]> nbOfElemsInt( CommInterface::ToIntArray<mcIdType>(nbOfElems,size) );
+      std::unique_ptr<int[]> offsetsIn( CommInterface::ComputeOffset(nbOfElemsInt,size) );
+      this->allGatherV(array->begin(),nbOfCellsRequested,ParaTraits<T>::MPIDataType,result.get(),nbOfElemsInt.get(),offsetsIn.get(),ParaTraits<T>::MPIDataType,comm);
+      resultIndex = ComputeOffsetFull<mcIdType>(nbOfElems,size);
+      return size;
+    }
+
+    template<class T>
+    void allGatherArraysT2(MPI_Comm comm, const typename Traits<T>::ArrayType *array, std::vector< MCAuto<typename Traits<T>::ArrayType> >& arraysOut) const
+    {
+      using DataArrayT = typename Traits<T>::ArrayType;
+      std::unique_ptr<T[]> result;
+      std::unique_ptr<mcIdType[]> resultIndex;
+      int size(this->allGatherArraysT<T>(comm,array,result,resultIndex));
+      arraysOut.resize(size);
+      for(int i = 0 ; i < size ; ++i)
+      {
+        arraysOut[i] = DataArrayT::New();
+        mcIdType nbOfEltPack(resultIndex[i+1]-resultIndex[i]);
+        arraysOut[i]->alloc(nbOfEltPack,1);
+        std::copy(result.get()+resultIndex[i],result.get()+resultIndex[i+1],arraysOut[i]->getPointer());
+      }
+    }
+
     template<class T>
     int allToAllArraysT2(MPI_Comm comm, const std::vector< MCAuto<typename Traits<T>::ArrayType> >& arrays, MCAuto<typename Traits<T>::ArrayType>& arrayOut, std::unique_ptr<mcIdType[]>& nbOfElems2, mcIdType& nbOfComponents) const
     {
diff --git a/src/ParaMEDMEM/ParaDataArray.hxx b/src/ParaMEDMEM/ParaDataArray.hxx
index a6052f1a6..af7700738 100644
--- a/src/ParaMEDMEM/ParaDataArray.hxx
+++ b/src/ParaMEDMEM/ParaDataArray.hxx
@@ -42,7 +42,7 @@ namespace MEDCoupling
     std::size_t getHeapMemorySizeWithoutChildren() const override;
     std::vector<const BigMemoryObject *> getDirectChildrenWithNull() const override;
     void checkOKOneComponent(const std::string& msg);
-  private:
+  protected:
     MCAuto<typename Traits<T>::ArrayType> _seq_da;
   };
 
@@ -50,7 +50,7 @@ namespace MEDCoupling
   class ParaDataArrayDiscrete : public ParaDataArrayTemplate<T>
   {
   public:
-    typename Traits<T>::ArrayType *buildComplement() const;
+    typename Traits<T>::ArrayType *buildComplement(T nbOfElems) const;
   protected:
     ParaDataArrayDiscrete(typename Traits<T>::ArrayType *seqDa):ParaDataArrayTemplate<T>(seqDa) { }
   };
diff --git a/src/ParaMEDMEM/ParaDataArray.txx b/src/ParaMEDMEM/ParaDataArray.txx
index 83598a43f..898a41946 100644
--- a/src/ParaMEDMEM/ParaDataArray.txx
+++ b/src/ParaMEDMEM/ParaDataArray.txx
@@ -21,6 +21,7 @@
 #pragma once
 
 #include "ParaDataArray.hxx"
+#include "CommInterface.hxx"
 
 #include <sstream>
 
@@ -61,11 +62,40 @@ namespace MEDCoupling
   }
 
   /*!
-    Parallel version of DataArrayInt::buildComplement.
+    Parallel version of DataArrayInt::buildComplement. Returns result on proc 0. Not allocated DataArrayT is returned for all procs.
    */
   template<class T>
-  typename Traits<T>::ArrayType *ParaDataArrayDiscrete<T>::buildComplement() const
+  typename Traits<T>::ArrayType *ParaDataArrayDiscrete<T>::buildComplement(T nbOfElems) const
   {
+    using DataArrayT = typename Traits<T>::ArrayType;
+    this->checkOKOneComponent("ParaDataArray::buildComplement");
+    MPI_Comm comm(MPI_COMM_WORLD);
+    CommInterface ci;
+    int size;
+    ci.commSize(comm,&size);
+    std::vector< MCAuto<DataArrayT> > idsCaptured(size);
+    for(int curRk = 0 ; curRk < size ; ++curRk)
+    {
+      mcIdType curStart(0),curEnd(0);
+      DataArrayTools<T>::GetSlice(0,nbOfElems,1,static_cast<T>(curRk),static_cast<T>(size),curStart,curEnd);
+      MCAuto<DataArrayIdType> idsInGlobalIds(this->_seq_da->findIdsInRange(curStart,curEnd));
+      idsCaptured[curRk] = this->_seq_da->selectByTupleIdSafe(idsInGlobalIds->begin(),idsInGlobalIds->end());
+    }
+    // communication : 1 arrays are going to be all2allized : ids
+    MCAuto<DataArrayIdType> aggregatedIds;
+    {
+      std::vector< MCAuto<DataArrayT> > myRkIdsCaptured;
+      ci.allToAllArraysT<T>(comm,idsCaptured,myRkIdsCaptured);
+      aggregatedIds = DataArrayIdType::Aggregate(FromVecAutoToVecOfConst<DataArrayIdType>(myRkIdsCaptured));
+    }
+    aggregatedIds->sort();
+    aggregatedIds = aggregatedIds->buildUnique();
+    int rank(-1);
+    ci.commRank(comm,&rank);
+    T vmin(std::numeric_limits<T>::max()),vmax(-std::numeric_limits<T>::max());
+    DataArrayTools<T>::GetSlice(0,nbOfElems,1,static_cast<T>(rank),static_cast<T>(size),vmin,vmax);
+    MCAuto<DataArrayIdType> retIds(aggregatedIds->findIdsNotInRange(vmin,vmax));
+    MCAuto<DataArrayT> ret(aggregatedIds->selectByTupleIdSafe(retIds->begin(),retIds->end()));
     return nullptr;
   }