Commit 026853e6 authored by COULAUD Olivier's avatar COULAUD Olivier
Browse files

Fix bugs

parent dc431542
......@@ -65,8 +65,8 @@ std::size_t read_data(parallel_manager& para, const std::string& filename, Conta
scalfmm::tools::DistFmaGenericLoader<ValueType, Dimension> loader(filename, para, para.io_master());
const int number_of_particles = loader.getNumberOfParticles();
const int local_number_of_particles = loader.getMyNumberOfParticles();
const std::int64_t number_of_particles = loader.getNumberOfParticles();
std::int64_t local_number_of_particles = loader.getMyNumberOfParticles();
width = loader.getBoxWidth();
Centre = loader.getBoxCenter();
......@@ -96,9 +96,12 @@ template<typename ContainerType, typename Box_type, typename Distrib_type>
void build_distrib(parallel_manager& para, ContainerType& particle_container, std::size_t& number_of_particles,
const Box_type& box, const int& level, Distrib_type& index_blocs)
{
auto rank = para.get_process_id();
// using morton_type = std::size_t;
using index_type = std::int64_t;
using morton_distrib_type = std::array<index_type, 2>;
///
/// sort particles according to their Morton index
scalfmm::tree::let::sort_particles_inplace(para, particle_container, box, level);
......@@ -107,23 +110,34 @@ void build_distrib(parallel_manager& para, ContainerType& particle_container, st
std::vector<index_type> leafMortonIdx(localNumberOfParticles);
// As the particles are sorted the leafMortonIdx is sorted too
#pragma omp parallel for shared(localNumberOfParticles, box, level)
for(std::size_t part = 0; part < localNumberOfParticles; ++part)
for(std::size_t part = 0; part < particle_container.size(); ++part)
{
leafMortonIdx[part] = scalfmm::index::get_morton_index(particle_container[part].position(), box, level);
}
/// construct the particle distribution particles_distrib and in
/// leafMortonIdx whe have the lorton index of the leaf = blocs
/// leafMortonIdx whe have the morton index of the leaf = blocs
std::vector<morton_distrib_type> particles_distrib(para.get_num_processes());
#ifdef LEAVES_DIST
particles_distrib = std::move(scalfmm::tree::distrib::balanced_leaves(para, leafMortonIdx));
#else
particles_distrib = std::move(
scalfmm::tree::distrib::balanced_particles(para, particle_container, leafMortonIdx, number_of_particles));
// the morton index of the leaf
auto last = std::unique(leafMortonIdx.begin(), leafMortonIdx.end());
leafMortonIdx.erase(last, leafMortonIdx.end());
#endif
/// Put the particles on the good processor
scalfmm::tree::distrib::fit_particles_in_distrib(para, particle_container, leafMortonIdx, particles_distrib, box,
level, number_of_particles);
///
// //#pragma omp parallel for shared(localNumberOfParticles, box, level)
// for(std::size_t part = 0; part < particle_container.size(); ++part)
// {
// std::cout << " new (" << rank << ")" << part << " " << particle_container[part]
// << " m= " << scalfmm::index::get_morton_index(particle_container[part].position(), box, level)
// << std::endl;
// }
///
/// build the distribution of rows/column or columns per leaf or bloc )
///
auto nbBlocs = leafMortonIdx.size();
......@@ -131,7 +145,6 @@ void build_distrib(parallel_manager& para, ContainerType& particle_container, st
std::size_t row{0};
index_blocs[nbBlocs] = particle_container.size();
auto rank = para.get_process_id();
auto myDistrib = particles_distrib[rank];
// scalfmm::out::print("rank(" + std::to_string(rank) + ") particles_distrib: ", particles_distrib);
// scalfmm::out::print("rank(" + std::to_string(rank) + ") leafMortonIdx: ", leafMortonIdx);
......@@ -158,7 +171,27 @@ void build_distrib(parallel_manager& para, ContainerType& particle_container, st
// std::cout << "Bloc distribution of rows:\n ";
// scalfmm::out::print("rank(" + std::to_string(rank) + ") blocs: ", blocs);
}
void build_matrix() {}
template<typename ContainerType, typename MatrixKernelType>
void build_save_matrix(const std::string& output_file, ContainerType container, MatrixKernelType& mk)
{
std::ofstream file(output_file);
std::cout << "writing matrix in file " << output_file << "\n";
file.setf(std::ios::fixed);
file.precision(14);
auto NNZ = container.size() * (container.size() - 1) / 2;
file << "%%MatrixMarket matrix coordinate real symmetric\n%\n";
file << container.size() << " " << container.size() << " " << NNZ << std::endl;
for(std::size_t i = 1; i < container.size(); ++i)
{
for(std::size_t j = 0; j < i; ++j)
{
auto val = std::get<0>(mk.evaluate(container[i].position(), container[j].position()));
file << i << " " << j << " " << val << std::endl;
}
}
file.close();
}
template<std::size_t Dimension, typename... Parameters>
auto run(inria::tcli::parser<Parameters...> const& parser, parallel_manager& para) -> int
{
......@@ -221,9 +254,14 @@ auto run(inria::tcli::parser<Parameters...> const& parser, parallel_manager& par
scalfmm::out::print("rank(" + std::to_string(rank) + ") blocs distrib: ", index_dist);
/// Construct the bloc matrix
build_matrix();
matrix_kernel_type mk;
if(para.get_num_processes() == 1)
{
build_save_matrix(output_file, container, mk);
}
////
std::cout << " End run \n";
return 0;
}
......@@ -246,6 +284,7 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
constexpr std::size_t dim = 2;
return run<dim>(parser, para);
std::cerr << " out run \n";
break;
}
case 3:
......@@ -261,5 +300,4 @@ auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
break;
}
para.end();
}
......@@ -236,14 +236,19 @@ namespace scalfmm::tree
auto nb_proc = manager.get_num_processes();
std::vector<std::array<Morton_type, 2>> morton_distrib(nb_proc, {0, 0});
auto UmortonArray(morton_array);
auto last = std::unique(UmortonArray.begin(), UmortonArray.end());
UmortonArray.erase(last, UmortonArray.end());
std::vector<MortonDistrib_type> weight(UmortonArray.size(), {0, 0});
auto bad_index = morton_array.back() + 1;
auto LeafMortonIndex(morton_array);
auto last = std::unique(LeafMortonIndex.begin(), LeafMortonIndex.end());
LeafMortonIndex.erase(last, LeafMortonIndex.end());
/// LeafMortonIndex has the size of the number of leaves
/// weight = ({morton index, number of particles}) for each leaf
std::vector<MortonDistrib_type> weight(LeafMortonIndex.size(), {bad_index, 0});
std::size_t pos = 0;
weight[pos][0] = UmortonArray[pos];
{
weight[pos][0] = LeafMortonIndex[pos];
std::cout << scalfmm::colors::red
<< "leaf size: " << LeafMortonIndex.size() << std::endl;
{ // loop on the number of particles
for(std::size_t part = 0; part < morton_array.size(); ++part)
{
// std::cout << "part " << part << " " <<
......@@ -251,14 +256,14 @@ namespace scalfmm::tree
// leafMortonIdx[pos]
// << " " << weight[pos] <<
// std::endl;
while(morton_array[part] != UmortonArray[pos])
while(morton_array[part] != LeafMortonIndex[pos])
{
// std::cout << " new pos " << pos <<
// std::endl;
pos++;
}
weight[pos][1] += 1;
weight[pos][0] = UmortonArray[pos];
weight[pos][0] = LeafMortonIndex[pos];
}
/// out::print("rank(" + std::to_string(rank) + ") weight: ", weight);
}
......@@ -268,6 +273,8 @@ namespace scalfmm::tree
// // [min, max] On two consecutive proceeses we may have
// max[p] = min[p+1]
// // we remove such case
/// Find the minimal and maximal index for the current process
/// this index should be unique (check with the neigboors left and right)
auto minIndex{weight[0]};
auto maxIndex{weight[weight.size() - 1]};
if(nb_proc == 1)
......@@ -275,6 +282,7 @@ namespace scalfmm::tree
morton_distrib[0] = {minIndex[0], maxIndex[0]};
// return morton_distrib;
}
out::print("rank(" + std::to_string(rank) + ") weight initial: ", weight);
#ifdef SCALFMM_USE_MPI
......@@ -292,15 +300,15 @@ namespace scalfmm::tree
weight[0][1] += weight_prev[1];
}
out::print("rank(" + std::to_string(rank) + ") weight final: ", weight);
// out::print("rank(" + std::to_string(rank) + ") weight final: ", weight);
///
/// compute the number of particles in the leaves
int nb_part = 0;
for(int i = 0; i < weight.size(); ++i)
{
nb_part += weight[i][1];
}
// Now the Morton indexex are uniq on all processes
// Now the Morton indexes are unique on all processes
//////////////////////////////////////////////////////////////////
/// Construct a uniform distribution of the morton index
///
......@@ -316,23 +324,40 @@ namespace scalfmm::tree
part_distrib[0] = local;
/// print("rank(" + std::to_string(rank) + ") Distrib cells Index: ", part_distrib);
// std::cout << "rank(" << rank << ") Distrib Leaf Index: " <<
// nb_part << std::endl;
out::print("rank(" + std::to_string(rank) + ") 0 Distrib cells Index: ", part_distrib);
// std::cout << "rank(" << rank << ") local: " <<local[0]<<" " <<local[1]<<" " <<local[2] <<std::endl;
/// share the distribution on all processors
auto nb_elt = sizeof(local);
conf.comm.allgather(local.data(), nb_elt, MPI_CHAR, part_distrib[0].data(), nb_elt, MPI_CHAR /*, 0*/);
// print("rank(" + std::to_string(rank) + ") Distrib cells Index: ", part_distrib);
out::print("rank(" + std::to_string(rank) + ") Distrib cells Index: ", part_distrib);
///
/// Try to have the same number of particles on a processor
///
block = number_of_particles / nb_proc;
std::vector<int> tomove(nb_proc, 0);
std::vector<int> tosendL(nb_proc, 0), tosendR(nb_proc, 0);
std::vector<int> toreceiv(nb_proc, 0);
for(int i = 0; i < nb_proc; ++i)
{
tomove[i] = part_distrib[i][2] - block;
std::vector<Morton_type> numberLeaves(nb_proc, 0);
Morton_type maxLeaves{0}, minLeaves{0};
for (int i = 0; i < nb_proc; ++i) {
/// Prevent to have 0 cell on a processor.
numberLeaves[i] = part_distrib[i][1] - part_distrib[i][0] +1 ;
maxLeaves = std::max(numberLeaves[i], maxLeaves);
}
// print("rank(" + std::to_string(rank) + ") initial tomove: ", tomove);
out::print("rank(" + std::to_string(rank) + ") numberLeaves: ", numberLeaves);
std::cout << "rank(" + std::to_string(rank) + ") initial tomove: "
<< maxLeaves << std::endl;
/// Prevent to have 0 cell on a processor.
if (maxLeaves > 1) {
for(int i = 0; i < nb_proc; ++i)
{
tomove[i] = part_distrib[i][2] - block;
}
}
out::print("rank(" + std::to_string(rank) + ") initial tomove: ", tomove);
// if(rank == 0)
for(int i = 0; i < nb_proc - 1; ++i)
......@@ -358,10 +383,10 @@ namespace scalfmm::tree
}
tosendR[nb_proc - 1] = 0;
// print("rank(" + std::to_string(rank) + ") tomove: ", tomove);
// print("rank(" + std::to_string(rank) + ") tosendR: ", tosendR);
// print("rank(" + std::to_string(rank) + ") tosendRL: ", tosendL);
// ///
// out::print("rank(" + std::to_string(rank) + ") tomove: ", tomove);
// out::print("rank(" + std::to_string(rank) + ") tosendR: ", tosendR);
// out::print("rank(" + std::to_string(rank) + ") tosendRL: ", tosendL);
///
// std::cout << "tosendL(" + std::to_string(rank) + "): " << tosendL[rank] << std::endl;
// std::cout << "tosendR(" + std::to_string(rank) + "): " << tosendR[rank] << std::endl;
// if(rank > 0)
......@@ -379,7 +404,10 @@ namespace scalfmm::tree
///
int nb_leaf_to_left{0}, nb_leaf_to_right{0}, nb_part_to_left{0}, nb_part_to_right{0};
Morton_type morton_to_left{0}, morton_to_right{0};
MortonDistrib_type MortonPart_to_left{0}, MortonPart_to_right{0};
MortonDistrib_type MortonPart_to_left{{0, 0}},
MortonPart_to_right{{0, 0}};
// std::cout << rank << " Morton [ " << MortonPart_to_left<< ", " << MortonPart_to_right << "]" << std::endl;
if(tosendL[rank] > 0)
{
int leaf_idx = 0;
......@@ -426,7 +454,7 @@ namespace scalfmm::tree
// std::cout << rank << "send morton_to_right " << morton_to_right << std::endl;
}
local[3] = 0;
// std::cout << rank << " local partition [ " << local[0] << ", " << local[1] << "]" << std::endl;
std::cout << rank << " local partition [ " << local[0] << ", " << local[1] << "]" << std::endl;
/// Send the number
/// send to left and right
......@@ -490,7 +518,7 @@ namespace scalfmm::tree
// std::array<Morton_type, 3> local{weight[0][0], weight[weight.size() - 1][0], nb_part};
std::array<Morton_type, 2> local1 = {morton_from_left, morton_from_right};
// std::cout << rank << " final local 1 [ " << local1[0] << ", " << local1[1] << "]" << std::endl;
std::cout << rank << " final local 1 [ " << local1[0] << ", " << local1[1] << "]" << std::endl;
morton_distrib[0] = local1;
......@@ -624,18 +652,20 @@ namespace scalfmm::tree
int my_rank = manager.get_process_id();
int nb_proc = manager.get_num_processes();
#ifdef SCALFMM_USE_MPI
std::cout << " (" << my_rank << ") size " << particles.size() << " "
<< morton_array.size() << std::endl;
auto comm = manager.get_communicator();
// std::cout << "\n------------- fit_particles_in_distrib -------------" << std::endl;
// out::print("rank(" + std::to_string(my_rank) + ") morton_array: ", morton_array);
out::print("rank(" + std::to_string(my_rank) + ") morton_array: ", morton_array);
// get the min and the max morton index of the particles own by the
// process
// send the number of communication we will receive
auto mortonMin = morton_dist[my_rank][0];
auto mortonMax = morton_dist[my_rank][1];
auto to_comm = std::move(compute_communications(my_rank, particles, morton_array, morton_dist));
auto to_comm = std::move(compute_communications(
my_rank, particles, morton_array, morton_dist));
// std::cout << " (" << my_rank << ") " << std::get<0>(to_comm) << std::endl;
// Send these numbers
auto nb_message = std::get<0>(to_comm);
auto nb_length_left = std::get<1>(to_comm);
......@@ -766,14 +796,17 @@ namespace scalfmm::tree
delete[] buffer_right;
}
particles = std::move(newArray);
std::cout << my_rank << " local Particles size: " << particles.size() << std::endl;
// std::cout << my_rank << " local Particles size: " << particles.size() << std::endl;
///
/// Check if we still have th good number of particles
using int_length = std::int64_t;
int_length new_num_particles{0};
int_length local_num_particles{static_cast<int_length>(particles.size())};
auto mpi_type = inria::mpi::get_datatype<int_length>();
conf.comm.allreduce(&local_num_particles, &new_num_particles, 1, mpi_type, MPI_SUM);
conf.comm.allreduce(&local_num_particles, &new_num_particles, 1,
mpi_type, MPI_SUM);
// std::cout << " (" << my_rank << ") " << local_num_particles << " "
// << new_num_particles << std::endl;
if(new_num_particles - total_num_particles != 0)
{
......@@ -1425,9 +1458,14 @@ namespace scalfmm::tree
const int& leaf_level)
{
// out::print("rank(" + std::to_string(rank) + ") partArray ori: ", particle_container);
// auto rank = manager.get_process_id();
#ifdef SCALFMM_USE_MPI
inria::mpi_config conf(manager.get_communicator());
// inria::mpi_config conf(manager.get_communicator());
// std::int64_t new_num_particles, local_num_particles{static_cast<int64_t>(particle_container.size())};
// auto mpi_type = inria::mpi::get_datatype< std::int64_t>();
// conf.comm.allreduce(&local_num_particles, &new_num_particles, 1, mpi_type, MPI_SUM);
// std::cout << " (" << rank << ") size " << local_num_particles << " " << new_num_particles << std::endl;
inria::sort(manager.get_communicator(), particle_container,
[&box, &leaf_level](const auto& p1, const auto& p2) {
......@@ -1436,6 +1474,10 @@ namespace scalfmm::tree
return m1 < m2;
});
// local_num_particles = static_cast<size_t>(particle_container.size());
// conf.comm.allreduce(&local_num_particles, &new_num_particles, 1, mpi_type, MPI_SUM);
// std::cout << " (" << rank << ") size " << local_num_particles << " " << new_num_particles << std::endl;
#else
std::sort(particle_container.begin(), particle_container.end(),
[&box, &leaf_level](const auto& p1, const auto& p2) {
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment