geo_mpi_dual.cc
Go to the documentation of this file.
#include "rheolef/config.h"
#ifdef _RHEOLEF_HAVE_MPI
// mesh2dual
#include "geo_partition_scotch.h"
#include <algorithm> // fill, copy
namespace rheolef {
using namespace std;
static void ikeysort(int total_elems, KeyValueType *pbase);
// -------------------------------------------------------------------------
// This function converts a mesh into a dual graph
// -------------------------------------------------------------------------
void geo_dual (
my_idxtype *elmdist,
my_idxtype *eptr,
vector<my_idxtype>& eind,
int *ncommonnodes,
vector<my_idxtype>& xadj,
vector<my_idxtype>& adjncy,
const mpi::communicator& comm)
{
int i, j, jj, k, kk, m;
int pe, count, mask, pass;
int lnns, my_nns, node;
int firstelm, firstnode, lnode, nrecv, nsend;
my_idxtype maxcount;
int npes = comm.size();
int mype = comm.rank();
srand (mype);
int nelms = elmdist[mype+1]-elmdist[mype];
check_macro(nelms >= 0, "unexpected nelms <= 0");
mask = (1<<11)-1;
// ----------------------------
// 1) determine number of nodes
// ----------------------------
int gminnode = mpi::all_reduce (comm, eind[idxamin(eptr[nelms], eind)], mpi::minimum<int>());
for (i=0; i<eptr[nelms]; i++)
eind[i] -= gminnode;
int gmaxnode = mpi::all_reduce (comm, eind[idxamax(eptr[nelms], eind)], mpi::maximum<int>());
// ----------------------------------------
// 2) construct node distribution array
// ----------------------------------------
vector<my_idxtype> nodedist (npes+1, 0);
for (nodedist[0]=0, i=0, j=gmaxnode+1; i<npes; i++) {
k = j/(npes-i);
nodedist[i+1] = nodedist[i]+k;
j -= k;
}
my_nns = nodedist[mype+1]-nodedist[mype];
firstnode = nodedist[mype];
vector<KeyValueType> nodelist (eptr[nelms]);
vector<my_idxtype> auxarray (eptr[nelms]);
vector<my_idxtype> htable (amax(my_nns, mask+1), -1);
vector<int> scounts (4*npes+2);
int *rcounts = scounts.begin().operator->() + npes;
int *sdispl = scounts.begin().operator->() + 2*npes;
int *rdispl = scounts.begin().operator->() + 3*npes+1;
// -----------------------------------------------
// 3) first find a local numbering of the nodes
// -----------------------------------------------
for (i=0; i<nelms; i++) {
for (j=eptr[i]; j<eptr[i+1]; j++) {
nodelist[j].key = eind[j];
nodelist[j].val = j;
auxarray[j] = i; /* remember the local element ID that uses this node */
}
}
ikeysort(eptr[nelms], nodelist.begin().operator->());
for (count=1, i=1; i<eptr[nelms]; i++) {
if (nodelist[i].key > nodelist[i-1].key)
count++;
}
lnns = count;
vector<my_idxtype> nmap (lnns);
// -----------------------------------------------
// 4) renumber the nodes of the elements array
// -----------------------------------------------
count = 1;
nmap[0] = nodelist[0].key;
eind[nodelist[0].val] = 0;
nodelist[0].val = auxarray[nodelist[0].val]; /* Store the local element ID */
for (i=1; i<eptr[nelms]; i++) {
if (nodelist[i].key > nodelist[i-1].key) {
nmap[count] = nodelist[i].key;
count++;
}
eind[nodelist[i].val] = count-1;
nodelist[i].val = auxarray[nodelist[i].val]; /* Store the local element ID */
}
comm.barrier();
// --------------------------------------------------------
// 5) perform comms necessary to construct node-element list
// --------------------------------------------------------
fill (scounts.begin(), scounts.begin() + npes, 0);
for (pe=i=0; i<eptr[nelms]; i++) {
while (nodelist[i].key >= nodedist[pe+1])
pe++;
scounts[pe] += 2;
}
check_macro(pe < npes, "unexpected pe");
mpi::all_to_all (comm, scounts.begin().operator->(), rcounts);
icopy(npes, scounts.begin().operator->(), sdispl);
init_csr_ptr(npes, sdispl);
icopy(npes, rcounts, rdispl);
init_csr_ptr(npes, rdispl);
check_macro(sdispl[npes] == eptr[nelms]*2, "unexpected sdispl[]");
nrecv = rdispl[npes]/2;
vector<KeyValueType> recvbuffer (amax(1, nrecv));
// no boost::mpi equivalent:
MPI_Alltoallv(nodelist.begin().operator->(), scounts.begin().operator->(), sdispl, IDX_DATATYPE, recvbuffer.begin().operator->(), rcounts, rdispl, IDX_DATATYPE, comm);
// --------------------------------------------------------
// 6) construct global node-element list
// --------------------------------------------------------
vector<my_idxtype> gnptr (my_nns+1, 0);
for (i=0; i<npes; i++) {
for (j=rdispl[i]/2; j<rdispl[i+1]/2; j++) {
lnode = recvbuffer[j].key-firstnode;
check_macro(lnode >= 0 && lnode < my_nns, "unexpected lnode range")
gnptr[lnode]++;
}
}
init_csr_ptr (my_nns, gnptr.begin());
vector<my_idxtype> gnind (amax(1, gnptr[my_nns]));
for (pe=0; pe<npes; pe++) {
firstelm = elmdist[pe];
for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
lnode = recvbuffer[j].key-firstnode;
gnind[gnptr[lnode]++] = recvbuffer[j].val+firstelm;
}
}
SHIFTCSR(i, my_nns, gnptr);
// --------------------------------------------------------
// 7) send the node-element info to the relevant processors
// --------------------------------------------------------
fill (scounts.begin(), scounts.begin() + npes, 0);
/* use a hash table to ensure that each node is sent to a proc only once */
for (pe=0; pe<npes; pe++) {
for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
lnode = recvbuffer[j].key-firstnode;
if (htable[lnode] == -1) {
scounts[pe] += gnptr[lnode+1]-gnptr[lnode];
htable[lnode] = 1;
}
}
/* now reset the hash table */
for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
lnode = recvbuffer[j].key-firstnode;
htable[lnode] = -1;
}
}
mpi::all_to_all (comm, scounts.begin().operator->(), rcounts);
icopy(npes, scounts.begin().operator->(), sdispl);
init_csr_ptr(npes, sdispl);
// --------------------------------------------------------
// 8) create the send buffer
// --------------------------------------------------------
nsend = sdispl[npes];
nodelist.clear();
vector<my_idxtype> sbuffer (amax(1, nsend));
count = 0;
for (pe=0; pe<npes; pe++) {
for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
lnode = recvbuffer[j].key-firstnode;
if (htable[lnode] == -1) {
for (k=gnptr[lnode]; k<gnptr[lnode+1]; k++) {
if (k == gnptr[lnode])
sbuffer[count++] = -1*(gnind[k]+1);
else
sbuffer[count++] = gnind[k];
}
htable[lnode] = 1;
}
}
check_macro(count == sdispl[pe+1], "unexpected count");
/* now reset the hash table */
for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
lnode = recvbuffer[j].key-firstnode;
htable[lnode] = -1;
}
}
icopy(npes, rcounts, rdispl);
init_csr_ptr(npes, rdispl);
nrecv = rdispl[npes];
recvbuffer.clear();
vector<my_idxtype> rbuffer (amax(1, nrecv));
// no boost::mpi equivalent:
MPI_Alltoallv(sbuffer.begin().operator->(), scounts.begin().operator->(), sdispl, IDX_DATATYPE, rbuffer.begin().operator->(), rcounts, rdispl, IDX_DATATYPE, comm);
k = -1;
vector<my_idxtype> nptr (lnns+1, 0);
my_idxtype *nind = rbuffer.begin().operator->(); // QUOI ??? TODO: comprendre .... re-utilisation de la mem ?
for (pe=0; pe<npes; pe++) {
for (j=rdispl[pe]; j<rdispl[pe+1]; j++) {
if (nind[j] < 0) {
k++;
nind[j] = (-1*nind[j])-1;
}
nptr[k]++;
}
}
init_csr_ptr(lnns, nptr.begin());
check_macro(k+1 == lnns, "unexpected k+1");
check_macro(nptr[lnns] == nrecv, "unexpected nptr[]")
xadj.resize(nelms+1);
fill (xadj.begin(), xadj.end(), 0);
my_idxtype *myxadj = xadj.begin().operator->(); // TODO: suppress ptrs !
fill (htable.begin(), htable.begin() + mask+1, -1);
firstelm = elmdist[mype];
// ---------------------------------------------------------------------------
// 9) Two passes -- in first pass, simply find out the memory requirements
// ---------------------------------------------------------------------------
maxcount = 200;
vector<my_idxtype> ind (maxcount);
vector<my_idxtype> wgt (maxcount);
my_idxtype *myadjncy = NULL;
for (pass=0; pass<2; pass++) {
for (i=0; i<nelms; i++) {
for (count=0, j=eptr[i]; j<eptr[i+1]; j++) {
node = eind[j];
for (k=nptr[node]; k<nptr[node+1]; k++) {
if ((kk=nind[k]) == firstelm+i) continue;
m = htable[(kk&mask)];
if (m == -1) {
ind[count] = kk;
wgt[count] = 1;
htable[(kk&mask)] = count++;
} else {
if (ind[m] == kk) {
wgt[m]++;
} else {
for (jj=0; jj<count; jj++) {
if (ind[jj] == kk) {
wgt[jj]++;
break;
}
}
if (jj == count) {
ind[count] = kk;
wgt[count++] = 1;
}
}
}
// Adjust the memory.
if (count == maxcount-1) {
ind.resize (2*maxcount);
wgt.resize (2*maxcount);
maxcount *= 2;
}
}
}
for (j=0; j<count; j++) {
htable[(ind[j]&mask)] = -1;
if (wgt[j] >= *ncommonnodes) {
if (pass == 0)
myxadj[i]++;
else
myadjncy[myxadj[i]++] = ind[j];
}
}
}
if (pass == 0) {
init_csr_ptr(nelms, myxadj);
adjncy.resize (myxadj[nelms]);
myadjncy = adjncy.begin().operator->(); // TODO: avoid pointers !
}
else {
SHIFTCSR(i, nelms, myxadj);
}
}
// ---------------------------------------------------------------------------
// 10) correctly renumber the elements array */
// ---------------------------------------------------------------------------
for (i=0; i<eptr[nelms]; i++)
eind[i] = nmap[eind[i]] + gminnode;
}
// -------------------------------------------------------------------------
// quick sort algorithm
// -------------------------------------------------------------------------
/* Discontinue quicksort algorithm when partition gets below this size.
This particular magic number was chosen to work best on a Sun 4/260. */
#define MAX_THRESH 20
/* Byte-wise swap two items of size SIZE. */
#define QSSWAP(a, b, stmp) do { stmp = (a); (a) = (b); (b) = stmp; } while (0)
/* Stack node declarations used to store unfulfilled partition obligations. */
typedef struct {
KeyValueType *lo;
KeyValueType *hi;
} stack_node;
/* The next 4 #defines implement a very fast in-line stack abstraction. */
#define STACK_SIZE (8 * sizeof(unsigned long int))
#define PUSH(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
#define POP(low, high) ((void) (--top, (low = top->lo), (high = top->hi)))
#define STACK_NOT_EMPTY (stack < top)
static
void
ikeysort(int total_elems, KeyValueType *pbase) {
KeyValueType pivot, stmp;
if (total_elems == 0)
/* Avoid lossage with unsigned arithmetic below. */
return;
if (total_elems > MAX_THRESH) {
KeyValueType *lo = pbase;
KeyValueType *hi = &lo[total_elems - 1];
stack_node stack[STACK_SIZE]; /* Largest size needed for 32-bit int!!! */
stack_node *top = stack + 1;
while (STACK_NOT_EMPTY) {
KeyValueType *left_ptr;
KeyValueType *right_ptr;
KeyValueType *mid = lo + ((hi - lo) >> 1);
if (mid->key < lo->key)
QSSWAP(*mid, *lo, stmp);
if (hi->key < mid->key)
QSSWAP(*mid, *hi, stmp);
else
goto jump_over;
if (mid->key < lo->key)
QSSWAP(*mid, *lo, stmp);
jump_over:;
pivot = *mid;
left_ptr = lo + 1;
right_ptr = hi - 1;
/* Here's the famous ``collapse the walls'' section of quicksort.
Gotta like those tight inner loops! They are the main reason
that this algorithm runs much faster than others. */
do {
while (left_ptr->key < pivot.key)
left_ptr++;
while (pivot.key < right_ptr->key)
right_ptr--;
if (left_ptr < right_ptr) {
QSSWAP (*left_ptr, *right_ptr, stmp);
left_ptr++;
right_ptr--;
}
else if (left_ptr == right_ptr) {
left_ptr++;
right_ptr--;
break;
}
} while (left_ptr <= right_ptr);
/* Set up pointers for next iteration. First determine whether
left and right partitions are below the threshold size. If so,
ignore one or both. Otherwise, push the larger partition's
bounds on the stack and continue sorting the smaller one. */
if ((size_t) (right_ptr - lo) <= MAX_THRESH) {
if ((size_t) (hi - left_ptr) <= MAX_THRESH)
/* Ignore both small partitions. */
POP (lo, hi);
else
/* Ignore small left partition. */
lo = left_ptr;
}
else if ((size_t) (hi - left_ptr) <= MAX_THRESH)
/* Ignore small right partition. */
hi = right_ptr;
else if ((right_ptr - lo) > (hi - left_ptr)) {
/* Push larger left partition indices. */
PUSH (lo, right_ptr);
lo = left_ptr;
}
else {
/* Push larger right partition indices. */
PUSH (left_ptr, hi);
hi = right_ptr;
}
}
}
/* Once the BASE_PTR array is partially sorted by quicksort the rest
is completely sorted using insertion sort, since this is efficient
for partitions below MAX_THRESH size. BASE_PTR points to the beginning
of the array to sort, and END_PTR points at the very last element in
the array (*not* one beyond it!). */
{
KeyValueType *end_ptr = &pbase[total_elems - 1];
KeyValueType *tmp_ptr = pbase;
KeyValueType *thresh = (end_ptr < pbase + MAX_THRESH ? end_ptr : pbase + MAX_THRESH);
KeyValueType *run_ptr;
/* Find smallest element in first threshold and place it at the
array's beginning. This is the smallest array element,
and the operation speeds up insertion sort's inner loop. */
for (run_ptr = tmp_ptr + 1; run_ptr <= thresh; run_ptr++)
if (run_ptr->key < tmp_ptr->key)
tmp_ptr = run_ptr;
if (tmp_ptr != pbase)
QSSWAP(*tmp_ptr, *pbase, stmp);
/* Insertion sort, running from left-hand-side up to right-hand-side. */
run_ptr = pbase + 1;
while (++run_ptr <= end_ptr) {
tmp_ptr = run_ptr - 1;
while (run_ptr->key < tmp_ptr->key)
tmp_ptr--;
tmp_ptr++;
if (tmp_ptr != run_ptr) {
KeyValueType elmnt = *run_ptr;
KeyValueType *mptr;
for (mptr=run_ptr; mptr>tmp_ptr; mptr--)
*mptr = *(mptr-1);
*mptr = elmnt;
}
}
}
}
} // namespace rheolef
#endif // _RHEOLEF_HAVE_MPI
rheolef::my_idxtype
int my_idxtype
Definition: geo_partition_scotch.h:28
rheolef::KeyValueType
struct KeyValueType KeyValueType
Definition: geo_partition_scotch.h:122
MAX_THRESH
#define MAX_THRESH
Definition: geo_mpi_dual.cc:333
STACK_SIZE
#define STACK_SIZE
Definition: geo_mpi_dual.cc:346
PUSH
#define PUSH(low, high)
Definition: geo_mpi_dual.cc:347
POP
#define POP(low, high)
Definition: geo_mpi_dual.cc:348
QSSWAP
#define QSSWAP(a, b, stmp)
Definition: geo_mpi_dual.cc:336
STACK_NOT_EMPTY
#define STACK_NOT_EMPTY
Definition: geo_mpi_dual.cc:349
geo_partition_scotch.h
icopy
#define icopy(n, a, b)
Definition: geo_partition_scotch.h:125
IDX_DATATYPE
#define IDX_DATATYPE
Definition: geo_partition_scotch.h:29
SHIFTCSR
#define SHIFTCSR(i, n, a)
Definition: geo_partition_scotch.h:127
amax
#define amax(a, b)
Definition: geo_partition_scotch.h:124
mkgeo_sector.m
m
Definition: mkgeo_sector.sh:118
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::idxamax
int idxamax(int n, const std::vector< my_idxtype > &x)
Definition: geo_partition_scotch.h:160
rheolef::geo_dual
void geo_dual(my_idxtype *elmdist, my_idxtype *eptr, vector< my_idxtype > &eind, int *ncommonnodes, vector< my_idxtype > &xadj, vector< my_idxtype > &adjncy, const mpi::communicator &comm)
Definition: geo_mpi_dual.cc:35
rheolef::idxamin
int idxamin(int n, const std::vector< my_idxtype > &x)
Definition: geo_partition_scotch.h:150