rheolef  6.5
geo_mpi_dual.cc
Go to the documentation of this file.
1 #include "rheolef/config.h"
2 #ifdef _RHEOLEF_HAVE_MPI
3 #include "geo_partition_scotch.h"
4 #include <algorithm> // fill, copy
5 
6 namespace rheolef {
7 using namespace std;
8 
9 static void ikeysort(int total_elems, KeyValueType *pbase);
10 
11 void geo_dual (
12  my_idxtype *elmdist,
13  my_idxtype *eptr,
14  vector<my_idxtype>& eind,
15  int *ncommonnodes,
16  vector<my_idxtype>& xadj,
17  vector<my_idxtype>& adjncy,
18  const mpi::communicator& comm)
19 {
20  int i, j, jj, k, kk, m;
21  int pe, count, mask, pass;
22  int lnns, my_nns, node;
23  int firstelm, firstnode, lnode, nrecv, nsend;
24  my_idxtype maxcount;
25 
26  int npes = comm.size();
27  int mype = comm.rank();
28  srand (mype);
29 
30  int nelms = elmdist[mype+1]-elmdist[mype];
31  check_macro(nelms >= 0, "unexpected nelms <= 0");
32 
33  mask = (1<<11)-1;
34 
35  int gminnode = mpi::all_reduce (comm, eind[idxamin(eptr[nelms], eind)], mpi::minimum<int>());
36  for (i=0; i<eptr[nelms]; i++)
37  eind[i] -= gminnode;
38 
39  int gmaxnode = mpi::all_reduce (comm, eind[idxamax(eptr[nelms], eind)], mpi::maximum<int>());
40  vector<my_idxtype> nodedist (npes+1, 0);
41  for (nodedist[0]=0, i=0, j=gmaxnode+1; i<npes; i++) {
42  k = j/(npes-i);
43  nodedist[i+1] = nodedist[i]+k;
44  j -= k;
45  }
46  my_nns = nodedist[mype+1]-nodedist[mype];
47  firstnode = nodedist[mype];
48 
49  vector<KeyValueType> nodelist (eptr[nelms]);
50  vector<my_idxtype> auxarray (eptr[nelms]);
51  vector<my_idxtype> htable (amax(my_nns, mask+1), -1);
52  vector<int> scounts (4*npes+2);
53  int *rcounts = scounts.begin().operator->() + npes;
54  int *sdispl = scounts.begin().operator->() + 2*npes;
55  int *rdispl = scounts.begin().operator->() + 3*npes+1;
56  for (i=0; i<nelms; i++) {
57  for (j=eptr[i]; j<eptr[i+1]; j++) {
58  nodelist[j].key = eind[j];
59  nodelist[j].val = j;
60  auxarray[j] = i; /* remember the local element ID that uses this node */
61  }
62  }
63  ikeysort(eptr[nelms], nodelist.begin().operator->());
64 
65  for (count=1, i=1; i<eptr[nelms]; i++) {
66  if (nodelist[i].key > nodelist[i-1].key)
67  count++;
68  }
69 
70  lnns = count;
71  vector<my_idxtype> nmap (lnns);
72  count = 1;
73  nmap[0] = nodelist[0].key;
74  eind[nodelist[0].val] = 0;
75  nodelist[0].val = auxarray[nodelist[0].val]; /* Store the local element ID */
76  for (i=1; i<eptr[nelms]; i++) {
77  if (nodelist[i].key > nodelist[i-1].key) {
78  nmap[count] = nodelist[i].key;
79  count++;
80  }
81  eind[nodelist[i].val] = count-1;
82  nodelist[i].val = auxarray[nodelist[i].val]; /* Store the local element ID */
83  }
84  comm.barrier();
85  fill (scounts.begin(), scounts.begin() + npes, 0);
86  for (pe=i=0; i<eptr[nelms]; i++) {
87  while (nodelist[i].key >= nodedist[pe+1])
88  pe++;
89  scounts[pe] += 2;
90  }
91  check_macro(pe < npes, "unexpected pe");
92 
93  mpi::all_to_all (comm, scounts.begin().operator->(), rcounts);
94 
95  icopy(npes, scounts.begin().operator->(), sdispl);
96  init_csr_ptr(npes, sdispl);
97 
98  icopy(npes, rcounts, rdispl);
99  init_csr_ptr(npes, rdispl);
100 
101  check_macro(sdispl[npes] == eptr[nelms]*2, "unexpected sdispl[]");
102 
103  nrecv = rdispl[npes]/2;
104  vector<KeyValueType> recvbuffer (amax(1, nrecv));
105 
106  MPI_Alltoallv(nodelist.begin().operator->(), scounts.begin().operator->(), sdispl, IDX_DATATYPE, recvbuffer.begin().operator->(), rcounts, rdispl, IDX_DATATYPE, comm);
107  vector<my_idxtype> gnptr (my_nns+1, 0);
108  for (i=0; i<npes; i++) {
109  for (j=rdispl[i]/2; j<rdispl[i+1]/2; j++) {
110  lnode = recvbuffer[j].key-firstnode;
111  check_macro(lnode >= 0 && lnode < my_nns, "unexpected lnode range")
112 
113  gnptr[lnode]++;
114  }
115  }
116  init_csr_ptr (my_nns, gnptr.begin());
117 
118  vector<my_idxtype> gnind (amax(1, gnptr[my_nns]));
119  for (pe=0; pe<npes; pe++) {
120  firstelm = elmdist[pe];
121  for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
122  lnode = recvbuffer[j].key-firstnode;
123  gnind[gnptr[lnode]++] = recvbuffer[j].val+firstelm;
124  }
125  }
126  SHIFTCSR(i, my_nns, gnptr);
127  fill (scounts.begin(), scounts.begin() + npes, 0);
128 
129  /* use a hash table to ensure that each node is sent to a proc only once */
130  for (pe=0; pe<npes; pe++) {
131  for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
132  lnode = recvbuffer[j].key-firstnode;
133  if (htable[lnode] == -1) {
134  scounts[pe] += gnptr[lnode+1]-gnptr[lnode];
135  htable[lnode] = 1;
136  }
137  }
138 
139  /* now reset the hash table */
140  for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
141  lnode = recvbuffer[j].key-firstnode;
142  htable[lnode] = -1;
143  }
144  }
145  mpi::all_to_all (comm, scounts.begin().operator->(), rcounts);
146  icopy(npes, scounts.begin().operator->(), sdispl);
147  init_csr_ptr(npes, sdispl);
148 
149  nsend = sdispl[npes];
150  nodelist.clear();
151  vector<my_idxtype> sbuffer (amax(1, nsend));
152  count = 0;
153  for (pe=0; pe<npes; pe++) {
154  for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
155  lnode = recvbuffer[j].key-firstnode;
156  if (htable[lnode] == -1) {
157  for (k=gnptr[lnode]; k<gnptr[lnode+1]; k++) {
158  if (k == gnptr[lnode])
159  sbuffer[count++] = -1*(gnind[k]+1);
160  else
161  sbuffer[count++] = gnind[k];
162  }
163  htable[lnode] = 1;
164  }
165  }
166  check_macro(count == sdispl[pe+1], "unexpected count");
167 
168  /* now reset the hash table */
169  for (j=rdispl[pe]/2; j<rdispl[pe+1]/2; j++) {
170  lnode = recvbuffer[j].key-firstnode;
171  htable[lnode] = -1;
172  }
173  }
174 
175  icopy(npes, rcounts, rdispl);
176  init_csr_ptr(npes, rdispl);
177 
178  nrecv = rdispl[npes];
179  recvbuffer.clear();
180  vector<my_idxtype> rbuffer (amax(1, nrecv));
181 
182  MPI_Alltoallv(sbuffer.begin().operator->(), scounts.begin().operator->(), sdispl, IDX_DATATYPE, rbuffer.begin().operator->(), rcounts, rdispl, IDX_DATATYPE, comm);
183 
184  k = -1;
185  vector<my_idxtype> nptr (lnns+1, 0);
186  my_idxtype *nind = rbuffer.begin().operator->(); // QUOI ??? TODO: comprendre .... re-utilisation de la mem ?
187  for (pe=0; pe<npes; pe++) {
188  for (j=rdispl[pe]; j<rdispl[pe+1]; j++) {
189  if (nind[j] < 0) {
190  k++;
191  nind[j] = (-1*nind[j])-1;
192  }
193  nptr[k]++;
194  }
195  }
196  init_csr_ptr(lnns, nptr.begin());
197 
198  check_macro(k+1 == lnns, "unexpected k+1");
199  check_macro(nptr[lnns] == nrecv, "unexpected nptr[]")
200 
201  xadj.resize(nelms+1);
202  fill (xadj.begin(), xadj.end(), 0);
203  my_idxtype *myxadj = xadj.begin().operator->(); // TODO: suppress ptrs !
204  fill (htable.begin(), htable.begin() + mask+1, -1);
205  firstelm = elmdist[mype];
206 
207  maxcount = 200;
208  vector<my_idxtype> ind (maxcount);
209  vector<my_idxtype> wgt (maxcount);
210  my_idxtype *myadjncy = NULL;
211 
212  for (pass=0; pass<2; pass++) {
213  for (i=0; i<nelms; i++) {
214  for (count=0, j=eptr[i]; j<eptr[i+1]; j++) {
215  node = eind[j];
216 
217  for (k=nptr[node]; k<nptr[node+1]; k++) {
218  if ((kk=nind[k]) == firstelm+i) continue;
219  m = htable[(kk&mask)];
220  if (m == -1) {
221  ind[count] = kk;
222  wgt[count] = 1;
223  htable[(kk&mask)] = count++;
224  } else {
225  if (ind[m] == kk) {
226  wgt[m]++;
227  } else {
228  for (jj=0; jj<count; jj++) {
229  if (ind[jj] == kk) {
230  wgt[jj]++;
231  break;
232  }
233  }
234  if (jj == count) {
235  ind[count] = kk;
236  wgt[count++] = 1;
237  }
238  }
239  }
240  if (count == maxcount-1) {
241  ind.resize (2*maxcount);
242  wgt.resize (2*maxcount);
243  maxcount *= 2;
244  }
245  }
246  }
247  for (j=0; j<count; j++) {
248  htable[(ind[j]&mask)] = -1;
249  if (wgt[j] >= *ncommonnodes) {
250  if (pass == 0)
251  myxadj[i]++;
252  else
253  myadjncy[myxadj[i]++] = ind[j];
254  }
255  }
256  }
257 
258  if (pass == 0) {
259  init_csr_ptr(nelms, myxadj);
260  adjncy.resize (myxadj[nelms]);
261  myadjncy = adjncy.begin().operator->(); // TODO: avoid pointers !
262  }
263  else {
264  SHIFTCSR(i, nelms, myxadj);
265  }
266  }
267  // 10) correctly renumber the elements array */
268  for (i=0; i<eptr[nelms]; i++)
269  eind[i] = nmap[eind[i]] + gminnode;
270 }
271 
272 /* Discontinue quicksort algorithm when partition gets below this size.
273  This particular magic number was chosen to work best on a Sun 4/260. */
274 #define MAX_THRESH 20
275 
276 /* Byte-wise swap two items of size SIZE. */
277 #define QSSWAP(a, b, stmp) do { stmp = (a); (a) = (b); (b) = stmp; } while (0)
278 
279 /* Stack node declarations used to store unfulfilled partition obligations. */
280 typedef struct {
283 } stack_node;
284 
285 
286 /* The next 4 #defines implement a very fast in-line stack abstraction. */
287 #define STACK_SIZE (8 * sizeof(unsigned long int))
288 #define PUSH(low, high) ((void) ((top->lo = (low)), (top->hi = (high)), ++top))
289 #define POP(low, high) ((void) (--top, (low = top->lo), (high = top->hi)))
290 #define STACK_NOT_EMPTY (stack < top)
291 
292 
293 static
294 void
295 ikeysort(int total_elems, KeyValueType *pbase) {
296  KeyValueType pivot, stmp;
297  if (total_elems == 0)
298  /* Avoid lossage with unsigned arithmetic below. */
299  return;
300 
301  if (total_elems > MAX_THRESH) {
302  KeyValueType *lo = pbase;
303  KeyValueType *hi = &lo[total_elems - 1];
304  stack_node stack[STACK_SIZE]; /* Largest size needed for 32-bit int!!! */
305  stack_node *top = stack + 1;
306 
307  while (STACK_NOT_EMPTY) {
308  KeyValueType *left_ptr;
309  KeyValueType *right_ptr;
310  KeyValueType *mid = lo + ((hi - lo) >> 1);
311 
312  if (mid->key < lo->key)
313  QSSWAP(*mid, *lo, stmp);
314  if (hi->key < mid->key)
315  QSSWAP(*mid, *hi, stmp);
316  else
317  goto jump_over;
318  if (mid->key < lo->key)
319  QSSWAP(*mid, *lo, stmp);
320 
321 jump_over:;
322  pivot = *mid;
323  left_ptr = lo + 1;
324  right_ptr = hi - 1;
325 
326  /* Here's the famous ``collapse the walls'' section of quicksort.
327  Gotta like those tight inner loops! They are the main reason
328  that this algorithm runs much faster than others. */
329  do {
330  while (left_ptr->key < pivot.key)
331  left_ptr++;
332 
333  while (pivot.key < right_ptr->key)
334  right_ptr--;
335 
336  if (left_ptr < right_ptr) {
337  QSSWAP (*left_ptr, *right_ptr, stmp);
338  left_ptr++;
339  right_ptr--;
340  }
341  else if (left_ptr == right_ptr) {
342  left_ptr++;
343  right_ptr--;
344  break;
345  }
346  } while (left_ptr <= right_ptr);
347 
348  /* Set up pointers for next iteration. First determine whether
349  left and right partitions are below the threshold size. If so,
350  ignore one or both. Otherwise, push the larger partition's
351  bounds on the stack and continue sorting the smaller one. */
352 
353  if ((size_t) (right_ptr - lo) <= MAX_THRESH) {
354  if ((size_t) (hi - left_ptr) <= MAX_THRESH)
355  /* Ignore both small partitions. */
356  POP (lo, hi);
357  else
358  /* Ignore small left partition. */
359  lo = left_ptr;
360  }
361  else if ((size_t) (hi - left_ptr) <= MAX_THRESH)
362  /* Ignore small right partition. */
363  hi = right_ptr;
364  else if ((right_ptr - lo) > (hi - left_ptr)) {
365  /* Push larger left partition indices. */
366  PUSH (lo, right_ptr);
367  lo = left_ptr;
368  }
369  else {
370  /* Push larger right partition indices. */
371  PUSH (left_ptr, hi);
372  hi = right_ptr;
373  }
374  }
375  }
376  /* Once the BASE_PTR array is partially sorted by quicksort the rest
377  is completely sorted using insertion sort, since this is efficient
378  for partitions below MAX_THRESH size. BASE_PTR points to the beginning
379  of the array to sort, and END_PTR points at the very last element in
380  the array (*not* one beyond it!). */
381  {
382  KeyValueType *end_ptr = &pbase[total_elems - 1];
383  KeyValueType *tmp_ptr = pbase;
384  KeyValueType *thresh = (end_ptr < pbase + MAX_THRESH ? end_ptr : pbase + MAX_THRESH);
385  register KeyValueType *run_ptr;
386 
387  /* Find smallest element in first threshold and place it at the
388  array's beginning. This is the smallest array element,
389  and the operation speeds up insertion sort's inner loop. */
390 
391  for (run_ptr = tmp_ptr + 1; run_ptr <= thresh; run_ptr++)
392  if (run_ptr->key < tmp_ptr->key)
393  tmp_ptr = run_ptr;
394 
395  if (tmp_ptr != pbase)
396  QSSWAP(*tmp_ptr, *pbase, stmp);
397 
398  /* Insertion sort, running from left-hand-side up to right-hand-side. */
399  run_ptr = pbase + 1;
400  while (++run_ptr <= end_ptr) {
401  tmp_ptr = run_ptr - 1;
402  while (run_ptr->key < tmp_ptr->key)
403  tmp_ptr--;
404 
405  tmp_ptr++;
406  if (tmp_ptr != run_ptr) {
407  KeyValueType elmnt = *run_ptr;
408  KeyValueType *mptr;
409 
410  for (mptr=run_ptr; mptr>tmp_ptr; mptr--)
411  *mptr = *(mptr-1);
412  *mptr = elmnt;
413  }
414  }
415  }
416 }
417 
418 } // namespace rheolef
419 #endif // _RHEOLEF_HAVE_MPI
420