Update/debug input output (re-sort locally, update interval maker)

bfps/cpp/particles/abstract_particles_distr.hpp

@@ -467,7 +467,7 @@ public:
 
         const int offesetOutLow = (current_partition_size==1? nbOutLower : 0);
 
-        const int nbOutUpper = current_my_nb_particles_per_partition[current_partition_size-1] - offesetOutLow - (particles_utils::partition_extra<size_particle_positions>(
+        const int nbOutUpper = current_my_nb_particles_per_partition[current_partition_size-1] - offesetOutLow - particles_utils::partition_extra<size_particle_positions>(
                     &(*inout_positions_particles)[(current_offset_particles_for_partition[current_partition_size-1]+offesetOutLow)*size_particle_positions],
                     myTotalNbParticles - (current_offset_particles_for_partition[current_partition_size-1]+offesetOutLow),
                     [&](const real_number val[]){
@@ -485,7 +485,7 @@ public:
                               inout_rhs_particles[idx_rhs][idx2*size_particle_rhs + idx_val]);
                 }
             }
-        }));
+        }, (current_offset_particles_for_partition[current_partition_size-1]+offesetOutLow));
         DEBUG_MSG("[%d] nbOutUpper %d\n", my_rank, nbOutUpper);
 
         // Exchange number

bfps/cpp/particles/abstract_particles_output.hpp

@@ -30,6 +30,7 @@ class abstract_particles_output {
 
     std::unique_ptr<real_number[]> buffer_particles_positions_recv;
     std::vector<std::unique_ptr<real_number[]>> buffer_particles_rhs_recv;
+    std::unique_ptr<int[]> buffer_indexes_recv;
     int size_buffers_recv;
 
 
@@ -64,6 +65,7 @@ public:
         buffer_indexes_send.release();
         buffer_particles_positions_send.release();
         size_buffers_send = -1;
+        buffer_indexes_recv.release();
         buffer_particles_positions_recv.release();
         size_buffers_recv = -1;
         for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
@@ -80,7 +82,7 @@ public:
         DEBUG_MSG("[%d] nb_particles %d to distribute for saving \n", my_rank, nb_particles);
 
         {
-            TIMEZONE("sort");
+            TIMEZONE("sort-to-distribute");
 
             if(size_buffers_send < nb_particles && nb_particles){
                 buffer_indexes_send.reset(new std::pair<int,int>[nb_particles]);
@@ -120,17 +122,21 @@ public:
         const particles_utils::IntervalSplitter<int> particles_splitter(total_nb_particles, nb_processes, my_rank);
         DEBUG_MSG("[%d] nb_particles_per_proc %d for saving\n", my_rank, particles_splitter.getMySize());
 
+        int* buffer_indexes_send_tmp = reinterpret_cast<int*>(buffer_indexes_send.get());// trick re-use buffer_indexes_send memory
         std::vector<int> nb_particles_to_send(nb_processes, 0);
         for(int idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
             nb_particles_to_send[particles_splitter.getOwner(buffer_indexes_send[idx_part].second)] += 1;
+            buffer_indexes_send_tmp[idx_part] = buffer_indexes_send[idx_part].second;
         }
 
         alltoall_exchanger exchanger(mpi_com, std::move(nb_particles_to_send));
         // nb_particles_to_send is invalid after here
 
         const int nb_to_receive = exchanger.getTotalToRecv();
+        assert(nb_to_receive == particles_splitter.getMySize());
 
         if(size_buffers_recv < nb_to_receive && nb_to_receive){
+            buffer_indexes_recv.reset(new int[nb_to_receive]);
             buffer_particles_positions_recv.reset(new real_number[nb_to_receive*size_particle_positions]);
             for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
                 buffer_particles_rhs_recv[idx_rhs].reset(new real_number[nb_to_receive*size_particle_rhs]);
@@ -138,13 +144,47 @@ public:
             size_buffers_recv = nb_to_receive;
         }
 
+        {
+            TIMEZONE("exchange");
             // Could be done with multiple asynchronous coms
-        exchanger.alltoallv(buffer_particles_positions_send.get(), buffer_particles_positions_recv.get(), size_particle_positions);
+            exchanger.alltoallv<int>(buffer_indexes_send_tmp, buffer_indexes_recv.get());
+            exchanger.alltoallv<real_number>(buffer_particles_positions_send.get(), buffer_particles_positions_recv.get(), size_particle_positions);
+            for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
+                exchanger.alltoallv<real_number>(buffer_particles_rhs_send[idx_rhs].get(), buffer_particles_rhs_recv[idx_rhs].get(), size_particle_rhs);
+            }
+        }
+
+        if(size_buffers_send < nb_to_receive && nb_to_receive){
+            buffer_indexes_send.reset(new std::pair<int,int>[nb_to_receive]);
+            buffer_particles_positions_send.reset(new real_number[nb_to_receive*size_particle_positions]);
             for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
-            exchanger.alltoallv(buffer_particles_rhs_send[idx_rhs].get(), buffer_particles_rhs_recv[idx_rhs].get(), size_particle_rhs);
+                buffer_particles_rhs_send[idx_rhs].reset(new real_number[nb_to_receive*size_particle_rhs]);
+            }
+            size_buffers_send = nb_to_receive;
+        }
+
+        {
+            TIMEZONE("copy-local-order");
+            for(int idx_part = 0 ; idx_part < nb_to_receive ; ++idx_part){
+                const int src_idx = idx_part;
+                const int dst_idx = buffer_indexes_recv[idx_part]-particles_splitter.getMyOffset();
+                assert(0 <= dst_idx);
+                assert(dst_idx < particles_splitter.getMySize());
+
+                for(int idx_val = 0 ; idx_val < size_particle_positions ; ++idx_val){
+                    buffer_particles_positions_send[dst_idx*size_particle_positions + idx_val]
+                            = buffer_particles_positions_recv[src_idx*size_particle_positions + idx_val];
+                }
+                for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
+                    for(int idx_val = 0 ; idx_val < int(size_particle_rhs) ; ++idx_val){
+                        buffer_particles_rhs_send[idx_rhs][dst_idx*size_particle_rhs + idx_val]
+                                = buffer_particles_rhs_recv[idx_rhs][src_idx*size_particle_rhs + idx_val];
+                    }
+                }
+            }
         }
 
-        write(idx_time_step, buffer_particles_positions_recv.get(), buffer_particles_rhs_recv.data(),
+        write(idx_time_step, buffer_particles_positions_send.get(), buffer_particles_rhs_send.data(),
               nb_to_receive, particles_splitter.getMyOffset());
     }
 

bfps/cpp/particles/alltoall_exchanger.hpp

@@ -27,6 +27,7 @@ class alltoall_exchanger {
 public:
     alltoall_exchanger(const MPI_Comm& in_mpi_com, std::vector<int>/*no ref to move here*/ in_nb_items_to_send)
         :mpi_com(in_mpi_com), nb_items_to_send(std::move(in_nb_items_to_send)), total_to_recv(0){
+        TIMEZONE("alltoall_exchanger::constructor");
 
         AssertMpi(MPI_Comm_rank(mpi_com, &my_rank));
         AssertMpi(MPI_Comm_size(mpi_com, &nb_processes));

bfps/cpp/particles/particles_input_hdf5.hpp

@@ -237,7 +237,7 @@ public:
                                       split_particles_rhs[idx_rhs][idx2*size_particle_rhs + idx_val]);
                         }
                     }
-                });
+                }, previousOffset);
 
                 nb_particles_per_proc[idx_proc] = localOffset;
                 previousOffset += localOffset;

bfps/cpp/particles/particles_utils.hpp

@@ -47,7 +47,7 @@ inline int partition(real_number* array, const int size, Predicate pdc)
     for(int idx = idxInsert ; idx < size ; ++idx){
         if(pdc(&array[idx*nb_values])){
             for(int idxVal = 0 ; idxVal < nb_values ; ++idxVal){
-                std::iter_swap(array[idx*nb_values + idxVal], array[idxInsert*nb_values + idxVal]);
+                std::swap(array[idx*nb_values + idxVal], array[idxInsert*nb_values + idxVal]);
             }
             idxInsert += 1;
         }
@@ -188,8 +188,7 @@ public:
         else{
             step_split = double(nb_items)/double(nb_intervals);
             offset_mine = NumType(step_split*double(my_idx));
-            size_mine = NumType(step_split*double(my_idx+1))-offset_mine;
-            assert(my_idx != nb_intervals-1 || (offset_mine+size_mine) == nb_items);
+            size_mine = (my_idx != nb_intervals-1 ? NumType(step_split*double(my_idx+1)) : nb_items) -offset_mine;
         }
     }
 
@@ -210,8 +209,13 @@ public:
     }
 
     NumType getOwner(const NumType in_item_idx) const {
-        const NumType owner = NumType(double(in_item_idx)/step_split);
+        NumType owner = NumType(double(in_item_idx)/step_split);
+        if(owner != nb_intervals-1 && NumType(step_split*double(owner+1)) <= in_item_idx){
+            owner += 1;
+        }
         assert(owner < nb_intervals);
+        assert(IntervalSplitter(nb_items, nb_intervals, owner).getMyOffset() <= in_item_idx);
+        assert(in_item_idx < IntervalSplitter(nb_items, nb_intervals, owner).getMySize()+IntervalSplitter(nb_items, nb_intervals, owner).getMyOffset());
         return owner;
     }
 };