src/utility/parallel_tools.h

Go to the documentation of this file.

#ifndef ParallelTools_h
#define ParallelTools_h
#include <cassert>
#include <algorithm>
#include <iostream>
#include <map>
#include <utility>
#include <sstream>
#include <vector>
#include "mpi.h"
#include "utility/triple.h"
#include "tetmesh/entity.h"
#include "myrlog.h"

namespace parallel {

  //using namespace std;
using mesh::id_t;

template <typename T>
class Mpi_DType {
};

template <>
class Mpi_DType<int> {
    static const MPI_Datatype mpi_dtype;
};

template <>
class Mpi_DType<unsigned int> {
    static const MPI_Datatype mpi_dtype;
};

template <>
class Mpi_DType<unsigned long> {
    static const MPI_Datatype mpi_dtype;
};

/*
 * Template specialisations for determine_proc_of_this_element.
 */

template<typename T>
int determine_proc_of_this_element(T nodes, std::vector<int>& procs);

template<>
int determine_proc_of_this_element<std::pair<id_t, id_t> >(std::pair<id_t, id_t> nodes, std::vector<int>& procs);

template<>
int determine_proc_of_this_element<utility::triple<id_t, id_t, id_t> >(utility::triple<id_t, id_t, id_t> nodes, std::vector<int>& procs);

template<>
int determine_proc_of_this_element<id_t>(id_t node, std::vector<int>& procs);

template <typename T>
int assign_gids(std::vector<T>& elements, 
                std::vector<id_t>& gids, 
                std::vector<int>& elem_counts,
                std::vector<int>& elem_displs,
                MPI_Comm comm) 
{
    gids.clear();
        
    int mpi_rank, mpi_size; 
    MPI_Comm_rank(comm, &mpi_rank);      
    MPI_Comm_size(comm, &mpi_size);
    
    // Sort elements to facilitate a binary search later on
    std::sort(elements.begin(), elements.end());
    
    // Obtain recvbuf containing all T's on all processors. The processor ID of the processor 
    // that sent a particular T is availible from MPI support stuff
    std::vector<int> recvcounts, recvbuf_displs;
    unsigned int recvbuf_size = set_counts_displs(elements.size(), recvcounts, recvbuf_displs, comm);                     
    std::vector<T> recvbuf(recvbuf_size);
    
    // Byte stuff.
    int sendcount_bytes = elements.size() * sizeof(T);
    int* recvcounts_bytes = new int[mpi_size];
    for (int i = 0; i < mpi_size; i ++)
        recvcounts_bytes[i] = recvcounts[i] * sizeof(T);   
    
    int* displs_bytes = new int[mpi_size];
    displs_bytes[0] = 0;
    for (int i = 1; i < mpi_size; i++)
        displs_bytes[i] = displs_bytes[i-1] + recvcounts_bytes[i-1];            
    
    MPI_Allgatherv (&elements[0], sendcount_bytes, MPI_BYTE, &recvbuf[0], recvcounts_bytes, displs_bytes, MPI_BYTE, comm);
    delete[] recvcounts_bytes;
    delete[] displs_bytes;
    
    gids.resize(recvbuf_size);
    
    // Construct multimap containing as keys the T objects 
    // and as values the processor ID's of processors claiming that object
    std::multimap<T,int> procs_claiming_element;
    typename std::vector<T>::iterator current_elem = recvbuf.begin();
    for (int proc_id = 0; proc_id < mpi_size; proc_id++) {
        for (int elem_lid = 0; elem_lid < recvcounts[proc_id]; elem_lid++) {
            std::pair<T, int> p(*current_elem, proc_id);
            procs_claiming_element.insert(p);
            ++current_elem;
        }
    }
    
    /*  Determine processor for each unique T by creating a vector of size 
    recvbuf containing the target processor of the elements in recvbuf */     
    std::vector<int> proc_owning_element;
    for (typename std::vector<T>::iterator it = recvbuf.begin(); it != recvbuf.end(); ++it) {
        std::vector<int> claimers;
        std::pair<typename std::multimap<T,int>::iterator, typename std::multimap<T,int>::iterator> range = procs_claiming_element.equal_range(*it);
        for (typename std::multimap<T,int>::iterator j = range.first; j != range.second; ++j)
            claimers.push_back((*j).second);
        assert (claimers.size() >= 2);
        proc_owning_element.push_back(determine_proc_of_this_element<T>(*it, claimers));
    }
    
#if 0
    for (typename std::multimap<T,int>::iterator it = procs_claiming_element.begin(); it != procs_claiming_element.end(); ++it) {
        ostringstream buf;
        buf << (*it).first;
        rDebugAll("elem %s claimed by p%02d", buf.str().c_str(), (*it).second);
    }
    for (unsigned int i = 0; i < proc_owning_element.size(); i++) 
        cout << proc_owning_element[i] << " "; cout << endl;
#endif
    
    // Calculate number of elements per processor and global ID offset for each processor
    elem_counts.clear();
    elem_counts.resize(mpi_size);
    id_t element_counter = 0;
    for (int i_proc = 0; i_proc < mpi_size; i_proc++) {
        elem_counts[i_proc] = 0;
        for (int j_elem = 0; j_elem < recvcounts[i_proc]; j_elem++) {            
            if (proc_owning_element[element_counter] == i_proc)
                elem_counts[i_proc]++;
            element_counter ++;
        }
    }
    elem_displs.resize(mpi_size);
    elem_displs[0] = 0;
    for (int i = 1; i < mpi_size; i++)
        elem_displs[i] = elem_displs[i-1] + elem_counts[i-1];
    int num_global_ids = elem_displs[mpi_size - 1] + elem_counts[mpi_size - 1];

#if 0
    cout << "Recvcounts: "; for(int i = 0; i < mpi_size; i++) cout << recvcounts[i] << " "; cout << "." << endl;
    cout << "Elements_per_proc: "; for(int i = 0; i < mpi_size; i++) cout << elem_counts[i] << " "; cout << "." << endl;
    cout << "Offsets: "; for(int i = 0; i < mpi_size; i++) cout << elem_displs[i] << " "; cout << "." << endl;
#endif
    
    // Assign local ID's to local elements in recvbuf in a consecutive manner
    std::vector<id_t> unmoved_local_ids(recvbuf_size);    
    element_counter = 0;
    for (int i_proc = 0; i_proc < mpi_size; i_proc++) {
        int id_counter = 0;
        for (int j_elems = 0; j_elems < recvcounts[i_proc]; j_elems++) {
            if (proc_owning_element[element_counter] == i_proc) {
                unmoved_local_ids[element_counter] = id_counter;
                id_counter++;
            } else {
                unmoved_local_ids[element_counter] = mesh::ID_NONE;
            }
            element_counter++;
        }
        assert(id_counter == elem_counts[i_proc]);
    }
    assert (element_counter == recvbuf_size);
    
    // Assign global ID's of local elements in recvbuf (store in gids vector)
    element_counter = 0;
    int nof_assigned_elements = 0;
    std::vector<id_t> unassigned_elements;
    for (int i_proc = 0; i_proc < mpi_size; i_proc++) {
        for (int j_elems = 0; j_elems < recvcounts[i_proc]; j_elems++) {
            if (unmoved_local_ids[element_counter] != mesh::ID_NONE) {
                gids[element_counter] = unmoved_local_ids[element_counter] + elem_displs[i_proc];
                nof_assigned_elements++;
            } else {
                gids[element_counter] = mesh::ID_NONE;
                unassigned_elements.push_back(element_counter);
            }
        element_counter++;
        }
    }
    assert (element_counter == recvbuf_size);
    
    // Assign gids of remote elements in recvbuf (using binary search)
    for(std::vector<id_t>::iterator it = unassigned_elements.begin(); it != unassigned_elements.end(); ++it) {
       //for (typename std::vector<T>::iterator it = recvbuf.begin(); it != recvbuf.end(); ++it) {
        T elem = recvbuf[*it];
        int owner = proc_owning_element[*it];
        typename std::vector<T>::iterator search_start_point = recvbuf.begin() + recvbuf_displs[owner];
        typename std::vector<T>::iterator search_end_point = search_start_point + recvcounts[owner];
        std::pair<typename std::vector<T>::iterator, typename std::vector<T>::iterator> result = equal_range(search_start_point, search_end_point, elem);
        assert (result.second - result.first == 1); //Want to find exactly one element.
        gids[*it] = gids[result.first - recvbuf.begin()];
    }    
#if 0
    element_counter = 0;
    for (int i_proc = 0; i_proc < mpi_size; i_proc++) {
        for (int j_elems = 0; j_elems < recvcounts[i_proc]; j_elems++) {
            if (mpi_rank == 1)
                rDebugAll("iproc=%02d, gids[%2d]=%2d", i_proc, element_counter, gids[element_counter]);
            element_counter ++;
        }
    }    
#endif
    if (recvbuf_displs[mpi_rank] != 0) {
        for (int i = 0; i < recvcounts[mpi_rank]; ++ i)
            gids[i] = gids[i + recvbuf_displs[mpi_rank]];
        gids.resize(recvcounts[mpi_rank]);
    }
            
    return num_global_ids;
}

template <typename T>
int assign_gids2(bool* my_claims,
                 std::vector<T>& elements,
                 std::vector<id_t>& gids, 
                 std::vector<int>& elem_counts,
                 std::vector<int>& elem_displs,
                 MPI_Comm comm) 
{
    gids.clear();
        
    int mpi_rank, mpi_size; 
    MPI_Comm_rank(comm, &mpi_rank);      
    MPI_Comm_size(comm, &mpi_size);
        
    // ----------------------------------------------------------------------------------------
    // Gather "elements" from all processors:
    std::vector<int> recvcounts, recvbuf_displs;
    unsigned int recvbuf_size = set_counts_displs(elements.size(), recvcounts, recvbuf_displs, comm);                     
    std::vector<T> recvbuf(recvbuf_size);
    
    // Byte stuff.
    int sendcount_bytes = elements.size() * sizeof(T);
    int* recvcounts_bytes = new int[mpi_size];
    int* displs_bytes = new int[mpi_size];
    for (int i = 0; i < mpi_size; i ++) {
        recvcounts_bytes[i] = recvcounts[i] * sizeof(T);
        displs_bytes[i] = recvbuf_displs[i] * sizeof(T);
    }
            
    MPI_Allgatherv (&elements[0], sendcount_bytes, MPI_BYTE, &recvbuf[0], recvcounts_bytes, displs_bytes, MPI_BYTE, comm);
    delete[] recvcounts_bytes;
    delete[] displs_bytes;
    
    // ----------------------------------------------------------------------------------------
    // Gather claim vectors from all processors
    bool* all_claims = new bool[recvbuf_size];
    assert(sizeof(bool) == 1);
    MPI_Allgatherv(my_claims, elements.size(), MPI_BYTE, all_claims, &recvcounts[0], &recvbuf_displs[0], MPI_BYTE, comm);

    // ----------------------------------------------------------------------------------------
    // Construct multimap: (element) -> (ids claiming processors...)
    std::multimap<T,int> procs_claiming_element;
    typename std::vector<T>::iterator current_elem = recvbuf.begin();
    bool* current_claim = all_claims;
    for (int proc_id = 0; proc_id < mpi_size; proc_id++) {
        for (int elem_lid = 0; elem_lid < recvcounts[proc_id]; elem_lid++) {
            if (*current_claim) {
                std::pair<T, int> p(*current_elem, proc_id);
                procs_claiming_element.insert(p);
            }
            ++ current_elem;
            ++ current_claim;
        }
    }
    delete all_claims;

    // ----------------------------------------------------------------------------------------
    // Sort elements for each processor to facilitate a binary search later on.
    for (int proc_id = 0; proc_id < mpi_size; ++ proc_id) {
        typename std::vector<T>::iterator sort_start_point = recvbuf.begin() + recvbuf_displs[proc_id];
        typename std::vector<T>::iterator sort_end_point = sort_start_point + recvcounts[proc_id];
        std::sort(sort_start_point, sort_end_point);
    }
    // elements must be ordered in the same way as the corresponding block in recvbuf.
    std::sort(elements.begin(), elements.end());
    
    // ----------------------------------------------------------------------------------------
    // Determine owning processor for each element in revbuf.
    
    std::vector<int> proc_owning_element;
    for (typename std::vector<T>::iterator it = recvbuf.begin(); it != recvbuf.end(); ++it) {
        std::vector<int> claimers;
        std::pair<typename std::multimap<T,int>::iterator, typename std::multimap<T,int>::iterator> range = procs_claiming_element.equal_range(*it);
        for (typename std::multimap<T,int>::iterator j = range.first; j != range.second; ++j)
            claimers.push_back((*j).second);
        // At least one processor must claim the element.
#if 1
        if (claimers.size() == 0) {
            rErrorAll("edge (%d, %d) is not claimed", (*it).first, (*it).second);   
        }
#endif
        assert (claimers.size() >= 1); 
        proc_owning_element.push_back(determine_proc_of_this_element<T>(*it, claimers));
    }
    
#if 0
    for (typename std::multimap<T,int>::iterator it = procs_claiming_element.begin(); it != procs_claiming_element.end(); ++it) {
        ostringstream buf;
        buf << (*it).first;
        rDebugAll("elem %s claimed by p%02d", buf.str().c_str(), (*it).second);
    }
    for (unsigned int i = 0; i < proc_owning_element.size(); i++) 
        cout << proc_owning_element[i] << " "; cout << endl;
#endif
    
    // Calculate number of elements per processor and global ID offset for each processor
    elem_counts.clear();
    elem_counts.resize(mpi_size);
    id_t element_counter = 0;
    for (int i_proc = 0; i_proc < mpi_size; i_proc++) {
        elem_counts[i_proc] = 0;
        for (int j_elem = 0; j_elem < recvcounts[i_proc]; j_elem++) {            
            if (proc_owning_element[element_counter] == i_proc)
                elem_counts[i_proc] ++;
            element_counter ++;
        }
    }
    elem_displs.resize(mpi_size);
    elem_displs[0] = 0;
    for (int i = 1; i < mpi_size; i++)
        elem_displs[i] = elem_displs[i-1] + elem_counts[i-1];
    int num_global_ids = elem_displs[mpi_size - 1] + elem_counts[mpi_size - 1];

#if 0
    cout << "Recvcounts: "; for(int i = 0; i < mpi_size; i++) cout << recvcounts[i] << " "; cout << "." << endl;
    cout << "Elements_per_proc: "; for(int i = 0; i < mpi_size; i++) cout << elem_counts[i] << " "; cout << "." << endl;
    cout << "Offsets: "; for(int i = 0; i < mpi_size; i++) cout << elem_displs[i] << " "; cout << "." << endl;
#endif
    
    // Assign local ID's to local elements in recvbuf in a consecutive manner
    std::vector<id_t> unmoved_local_ids(recvbuf_size);    
    element_counter = 0;
    for (int i_proc = 0; i_proc < mpi_size; i_proc++) {
        int id_counter = 0;
        for (int j_elems = 0; j_elems < recvcounts[i_proc]; j_elems++) {
            if (proc_owning_element[element_counter] == i_proc) {
                unmoved_local_ids[element_counter] = id_counter;
                id_counter++;
            } else {
                unmoved_local_ids[element_counter] = mesh::ID_NONE;
            }
            element_counter++;
        }
        assert(id_counter == elem_counts[i_proc]);
    }
    assert (element_counter == recvbuf_size);
    
    // Assign global ID's of local elements in recvbuf (store in gids vector)
    gids.resize(recvbuf_size);
    element_counter = 0;
    int nof_assigned_elements = 0;
    std::vector<id_t> unassigned_elements;
    for (int i_proc = 0; i_proc < mpi_size; i_proc++) {
        for (int j_elems = 0; j_elems < recvcounts[i_proc]; j_elems++) {
            if (unmoved_local_ids[element_counter] != mesh::ID_NONE) {
                gids[element_counter] = unmoved_local_ids[element_counter] + elem_displs[i_proc];
                nof_assigned_elements++;
            } else {
                gids[element_counter] = mesh::ID_NONE;
                unassigned_elements.push_back(element_counter);
            }
        element_counter++;
        }
    }
    assert (element_counter == recvbuf_size);
    
    // Assign gids of remote elements in recvbuf (using binary search)
    for(std::vector<id_t>::iterator it = unassigned_elements.begin(); it != unassigned_elements.end(); ++it) {
       //for (typename std::vector<T>::iterator it = recvbuf.begin(); it != recvbuf.end(); ++it) {
        T elem = recvbuf[*it];
        int owner = proc_owning_element[*it];
        typename std::vector<T>::iterator search_start_point = recvbuf.begin() + recvbuf_displs[owner];
        typename std::vector<T>::iterator search_end_point = search_start_point + recvcounts[owner];
        std::pair<typename std::vector<T>::iterator, typename std::vector<T>::iterator> result = equal_range(search_start_point, search_end_point, elem);
        assert(result.second - result.first == 1); // Want to find exactly one element.
        gids[*it] = gids[result.first - recvbuf.begin()];
    }    
#if 0
    element_counter = 0;
    for (int i_proc = 0; i_proc < mpi_size; i_proc++) {
        for (int j_elems = 0; j_elems < recvcounts[i_proc]; j_elems++) {
            if (mpi_rank == 1)
                rDebugAll("iproc=%02d, gids[%2d]=%2d", i_proc, element_counter, gids[element_counter]);
            element_counter ++;
        }
    }    
#endif
    if (recvbuf_displs[mpi_rank] != 0) {
        for (int i = 0; i < recvcounts[mpi_rank]; ++ i)
            gids[i] = gids[i + recvbuf_displs[mpi_rank]];
        gids.resize(recvcounts[mpi_rank]);
    }
            
    return num_global_ids;
}

int set_ownership(const std::vector<id_t>& gids, std::vector<int>& owned_by, MPI_Comm comm);

void remove_gids(std::vector<id_t> gids, std::vector<id_t> eliminated_gids, std::vector<id_t> map_elimated);

void compute_mgid_map(id_t num_global_gids,
                      const std::vector<id_t>& gids,
                      std::vector<id_t>& eliminated_gids,
                      id_t* num_global_mgids,
                      std::map<id_t,id_t>& mgid_map,
                      MPI_Comm comm);

int set_counts_displs(int my_count, 
                      std::vector<int>& counts, 
                      std::vector<int>& displs, 
                      MPI_Comm comm);
                      
} // namespace parallel

std::ostream& operator << (std::ostream& ostr, const utility::triple<mesh::id_t, mesh::id_t, mesh::id_t>& x);

std::ostream& operator << (std::ostream& ostr, const std::pair<mesh::id_t, mesh::id_t>& x);

#endif

---