src/utility/parallel_tools.cpp

Go to the documentation of this file.

#include <numeric>
#include "parallel_tools.h"

namespace parallel {

const MPI_Datatype Mpi_DType<int>::mpi_dtype(MPI_INT);
const MPI_Datatype Mpi_DType<unsigned int>::mpi_dtype(MPI_UNSIGNED);
const MPI_Datatype Mpi_DType<unsigned long>::mpi_dtype(MPI_UNSIGNED_LONG);

template<>
int determine_proc_of_this_element<std::pair<id_t, id_t> >(std::pair<id_t, id_t> nodes, std::vector<int>& procs) {
    std::sort(procs.begin(), procs.end());
    id_t target = (nodes.first + nodes.second) % procs.size();
    assert (target <= procs.size());
    return procs[target];
}

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) {
    std::sort(procs.begin(), procs.end());
    id_t target = (nodes.first + nodes.second + nodes.third) % procs.size();
    assert (target <= procs.size());
    return procs[target];
}

template<>
int determine_proc_of_this_element<id_t>(id_t node, std::vector<int>& procs) {
    std::sort(procs.begin(), procs.end());
    id_t target = node % procs.size();
    assert (target <= procs.size());
    return procs[target];
}

std::ostream& operator << (std::ostream& ostr, const utility::triple<id_t, id_t, id_t>& x) {
    ostr << "Triple (face) (";
    ostr << x.first << ", ";
    ostr << x.second << ", ";
    ostr << x.third << ")";
    return ostr;
}

std::ostream& operator << (std::ostream& ostr, const std::pair<id_t, id_t>& x) {
    ostr << "Pair (edge) ("; 
    ostr << x.first << ", ";
    ostr << x.second << ")";
    return ostr;
}

int set_ownership(const std::vector<id_t>& gids, std::vector<int>& owned_by, MPI_Comm comm) {
    int mpi_rank, mpi_size; 
    MPI_Comm_rank(comm, &mpi_rank);      
    MPI_Comm_size(comm, &mpi_size);

    // Gather gids from all processors
    std::vector<int> recv_counts, recv_displs;
    unsigned int recv_size = set_counts_displs(gids.size(), recv_counts, recv_displs, comm);                     
    std::vector<id_t> recvbuf(recv_size);
    
    int* recvcounts_bytes = new int[mpi_size];
    for (int i = 0; i < mpi_size; i ++)
        recvcounts_bytes[i] = recv_counts[i] * sizeof(id_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((void *)&gids[0], gids.size()*sizeof(id_t), MPI_BYTE, 
                   &recvbuf[0], recvcounts_bytes, displs_bytes, MPI_BYTE, comm);
    delete[] recvcounts_bytes;
    delete[] displs_bytes;
    
    // Create multimap key:gid, value:ids of processors claiming gid
    std::multimap<id_t,int> procs_claiming_gid;
    std::vector<id_t>::iterator current_elem = recvbuf.begin();
    for (int proc_id = 0; proc_id < mpi_size; ++ proc_id) {
        for (int elem_lid = 0; elem_lid < recv_counts[proc_id]; ++ elem_lid) {
            std::pair<id_t, int> p(*current_elem, proc_id);
            procs_claiming_gid.insert(p);
            ++ current_elem;
        }
    }
    
    // Determine owning processors for all gids used locally
    int num_owned = 0;
    owned_by.clear();
    for (std::vector<id_t>::const_iterator it = gids.begin(); it != gids.end(); ++ it) {
        std::vector<int> claimers;
        std::pair<std::multimap<id_t,int>::iterator, std::multimap<id_t,int>::iterator> range = procs_claiming_gid.equal_range(*it);
        for (std::multimap<id_t,int>::iterator j = range.first; j != range.second; ++j)
            claimers.push_back((*j).second);
        int owner = determine_proc_of_this_element(*it, claimers);
        owned_by.push_back(owner);
        if (owner == mpi_rank)
            num_owned += 1;
    }
    
    return num_owned;
}

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 mpi_size;
    MPI_Comm_size(comm, &mpi_size);

    // Allgatherv eliminated_gids
    int num_my_eliminated_gids = eliminated_gids.size();
    std::vector<int> recvcnts(mpi_size);
    std::vector<int> displs(mpi_size);
    int sum_eliminated_gids = set_counts_displs(num_my_eliminated_gids, recvcnts, displs, comm);
    std::vector<id_t> eliminated_gid_recv(sum_eliminated_gids);
    MPI_Allgatherv(&eliminated_gids[0], num_my_eliminated_gids, MPI_INT,
                   &eliminated_gid_recv[0], &recvcnts[0], &displs[0], MPI_INT, comm);
#if 0                   
    for (int i = 0; i < sum_eliminated_gids; ++ i)
        rDebug("eliminated_gid_recv[%02d]        = %2d", i, eliminated_gid_recv[i]);
#endif
    
    // Sort, unique
    std::sort(eliminated_gid_recv.begin(), eliminated_gid_recv.end());
#if 0                   
    for (int i = 0; i < sum_eliminated_gids; ++ i)
        rDebug("eliminated_gid_recv_sort[%02d]   = %2d", i, eliminated_gid_recv[i]);
#endif
    std::vector<id_t>::iterator new_last;
    new_last = std::unique(eliminated_gid_recv.begin(), eliminated_gid_recv.end());
#if 0                   
    for (int i = 0; i < new_last - eliminated_gid_recv.begin(); ++ i)
        rDebug("eliminated_gid_recv_unique[%02d] = %2d", i, eliminated_gid_recv[i]);
#endif

    // Compute processor-specific mgid map containing only gids required by the processor
    mgid_map.clear();
    std::vector<id_t>::iterator insert_pos;
    for (std::vector<id_t>::const_iterator gid_it=gids.begin(); gid_it != gids.end(); ++ gid_it) {
        id_t gid = *gid_it;
        insert_pos = std::lower_bound(eliminated_gid_recv.begin(), new_last, gid);
        if (insert_pos == new_last || *insert_pos != *gid_it) {
            // gid not eliminated
            // FIXME: this should also work if insert_pos == new_last
            mgid_map[gid] = gid - (insert_pos - eliminated_gid_recv.begin());
        }
    }
    *num_global_mgids = num_global_gids - (new_last - eliminated_gid_recv.begin());
}

int set_counts_displs(int my_count, 
                      std::vector<int>& counts, 
                      std::vector<int>& displs, 
                      MPI_Comm comm) 
{
    int mpi_size; 
    MPI_Comm_size(comm, &mpi_size);
    counts.resize(mpi_size);
    displs.resize(mpi_size);
    MPI_Allgather(&my_count, 1, MPI_INT, &counts[0], 1, MPI_INT, comm);
    
    int size = 0;
    for (int i = 0; i < mpi_size; ++ i)
        size += counts[i];    
    displs[0] = 0;
    for (int i = 1; i < mpi_size; ++ i)
        displs[i] = displs[i-1] + counts[i-1];
    return size;
}

} // namespace parallel

---