MFEM  v3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
communication.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include <mpi.h>
17 #ifdef __bgq__
18 #include <mpix.h>
19 #endif
20 
21 #include "array.hpp"
22 #include "table.hpp"
23 #include "sets.hpp"
24 #include "communication.hpp"
25 #include <iostream>
26 #include <map>
27 
28 using namespace std;
29 
30 namespace mfem
31 {
32 
33 GroupTopology::GroupTopology(const GroupTopology &gt)
34  : MyComm(gt.MyComm),
35  group_lproc(gt.group_lproc)
36 {
37  gt.groupmaster_lproc.Copy(groupmaster_lproc);
38  gt.lproc_proc.Copy(lproc_proc);
39  gt.group_mgroup.Copy(group_mgroup);
40 }
41 
42 void GroupTopology::ProcToLProc()
43 {
44  int NRanks;
45  MPI_Comm_size(MyComm, &NRanks);
46 
47  map<int, int> proc_lproc;
48 
49  int lproc_counter = 0;
50  for (int i = 0; i < group_lproc.Size_of_connections(); i++)
51  {
52  const pair<const int, int> p(group_lproc.GetJ()[i], lproc_counter);
53  if (proc_lproc.insert(p).second)
54  {
55  lproc_counter++;
56  }
57  }
58  // Note: group_lproc.GetJ()[0] == MyRank --> proc_lproc[MyRank] == 0
59 
60  lproc_proc.SetSize(lproc_counter);
61  for (map<int, int>::iterator it = proc_lproc.begin();
62  it != proc_lproc.end(); ++it)
63  {
64  lproc_proc[it->second] = it->first;
65  }
66 
67  for (int i = 0; i < group_lproc.Size_of_connections(); i++)
68  {
69  group_lproc.GetJ()[i] = proc_lproc[group_lproc.GetJ()[i]];
70  }
71 
72  for (int i = 0; i < NGroups(); i++)
73  {
74  groupmaster_lproc[i] = proc_lproc[groupmaster_lproc[i]];
75  }
76 }
77 
78 void GroupTopology::Create(ListOfIntegerSets &groups, int mpitag)
79 {
80  groups.AsTable(group_lproc); // group_lproc = group_proc
81 
82  Table group_mgroupandproc;
83  group_mgroupandproc.SetDims(NGroups(),
84  group_lproc.Size_of_connections() + NGroups());
85  for (int i = 0; i < NGroups(); i++)
86  {
87  int j = group_mgroupandproc.GetI()[i];
88  group_mgroupandproc.GetI()[i+1] = j + group_lproc.RowSize(i) + 1;
89  group_mgroupandproc.GetJ()[j] = i;
90  j++;
91  for (int k = group_lproc.GetI()[i];
92  j < group_mgroupandproc.GetI()[i+1]; j++, k++)
93  {
94  group_mgroupandproc.GetJ()[j] = group_lproc.GetJ()[k];
95  }
96  }
97 
98  // build groupmaster_lproc with lproc = proc
99  groupmaster_lproc.SetSize(NGroups());
100 
101  // simplest choice of the group owner
102  for (int i = 0; i < NGroups(); i++)
103  {
104  groupmaster_lproc[i] = groups.PickElementInSet(i);
105  }
106 
107  // load-balanced choice of the group owner, which however can lead to
108  // isolated dofs
109  // for (i = 0; i < NGroups(); i++)
110  // groupmaster_lproc[i] = groups.PickRandomElementInSet(i);
111 
112  ProcToLProc();
113 
114  // build group_mgroup
115  group_mgroup.SetSize(NGroups());
116 
117  int send_counter = 0;
118  int recv_counter = 0;
119  for (int i = 1; i < NGroups(); i++)
120  if (groupmaster_lproc[i] != 0) // we are not the master
121  {
122  recv_counter++;
123  }
124  else
125  {
126  send_counter += group_lproc.RowSize(i)-1;
127  }
128 
129  MPI_Request *requests = new MPI_Request[send_counter];
130  MPI_Status *statuses = new MPI_Status[send_counter];
131 
132  int max_recv_size = 0;
133  send_counter = 0;
134  for (int i = 1; i < NGroups(); i++)
135  {
136  if (groupmaster_lproc[i] == 0) // we are the master
137  {
138  group_mgroup[i] = i;
139 
140  for (int j = group_lproc.GetI()[i];
141  j < group_lproc.GetI()[i+1]; j++)
142  {
143  if (group_lproc.GetJ()[j] != 0)
144  {
145  MPI_Isend(group_mgroupandproc.GetRow(i),
146  group_mgroupandproc.RowSize(i),
147  MPI_INT,
148  lproc_proc[group_lproc.GetJ()[j]],
149  mpitag,
150  MyComm,
151  &requests[send_counter]);
152  send_counter++;
153  }
154  }
155  }
156  else // we are not the master
157  if (max_recv_size < group_lproc.RowSize(i))
158  {
159  max_recv_size = group_lproc.RowSize(i);
160  }
161  }
162  max_recv_size++;
163 
164  IntegerSet group;
165  if (recv_counter > 0)
166  {
167  int count;
168  MPI_Status status;
169  int *recv_buf = new int[max_recv_size];
170  for ( ; recv_counter > 0; recv_counter--)
171  {
172  MPI_Recv(recv_buf, max_recv_size, MPI_INT,
173  MPI_ANY_SOURCE, mpitag, MyComm, &status);
174 
175  MPI_Get_count(&status, MPI_INT, &count);
176 
177  group.Recreate(count-1, recv_buf+1);
178  int g = groups.Lookup(group);
179  group_mgroup[g] = recv_buf[0];
180 
181  if (lproc_proc[groupmaster_lproc[g]] != status.MPI_SOURCE)
182  {
183  cerr << "\n\n\nGroupTopology::GroupTopology: "
184  << MyRank() << ": ERROR\n\n\n" << endl;
185  mfem_error();
186  }
187  }
188  delete [] recv_buf;
189  }
190 
191  MPI_Waitall(send_counter, requests, statuses);
192 
193  delete [] statuses;
194  delete [] requests;
195 }
196 
197 
198 // specializations of the MPITypeMap static member
199 template<> const MPI_Datatype MPITypeMap<int>::mpi_type = MPI_INT;
200 template<> const MPI_Datatype MPITypeMap<double>::mpi_type = MPI_DOUBLE;
201 
202 
204  : gtopo(gt)
205 {
206  group_buf_size = 0;
207  requests = NULL;
208  statuses = NULL;
209 }
210 
212 {
213  group_ldof.MakeI(gtopo.NGroups());
214  for (int i = 0; i < ldof_group.Size(); i++)
215  {
216  int group = ldof_group[i];
217  if (group != 0)
218  {
219  group_ldof.AddAColumnInRow(group);
220  }
221  }
222  group_ldof.MakeJ();
223 
224  for (int i = 0; i < ldof_group.Size(); i++)
225  {
226  int group = ldof_group[i];
227  if (group != 0)
228  {
229  group_ldof.AddConnection(group, i);
230  }
231  }
232  group_ldof.ShiftUpI();
233 
234  Finalize();
235 }
236 
238 {
239  int request_counter = 0;
240 
241  for (int gr = 1; gr < group_ldof.Size(); gr++)
242  if (group_ldof.RowSize(gr) != 0)
243  {
244  int gr_requests;
245  if (!gtopo.IAmMaster(gr)) // we are not the master
246  {
247  gr_requests = 1;
248  }
249  else
250  {
251  gr_requests = gtopo.GetGroupSize(gr)-1;
252  }
253 
254  request_counter += gr_requests;
255  group_buf_size += gr_requests * group_ldof.RowSize(gr);
256  }
257 
258  requests = new MPI_Request[request_counter];
259  statuses = new MPI_Status[request_counter];
260 }
261 
262 template <class T>
264 {
265  if (group_buf_size == 0)
266  {
267  return;
268  }
269 
270  group_buf.SetSize(group_buf_size*sizeof(T));
271  T *buf = (T *)group_buf.GetData();
272 
273  int i, gr, request_counter = 0;
274 
275  for (gr = 1; gr < group_ldof.Size(); gr++)
276  {
277  const int nldofs = group_ldof.RowSize(gr);
278 
279  // ignore groups without dofs
280  if (nldofs == 0)
281  {
282  continue;
283  }
284 
285  if (!gtopo.IAmMaster(gr)) // we are not the master
286  {
287  MPI_Irecv(buf,
288  nldofs,
290  gtopo.GetGroupMasterRank(gr),
291  40822 + gtopo.GetGroupMasterGroup(gr),
292  gtopo.GetComm(),
293  &requests[request_counter]);
294  request_counter++;
295  }
296  else // we are the master
297  {
298  // fill send buffer
299  const int *ldofs = group_ldof.GetRow(gr);
300  for (i = 0; i < nldofs; i++)
301  {
302  buf[i] = ldata[ldofs[i]];
303  }
304 
305  const int gs = gtopo.GetGroupSize(gr);
306  const int *nbs = gtopo.GetGroup(gr);
307  for (i = 0; i < gs; i++)
308  {
309  if (nbs[i] != 0)
310  {
311  MPI_Isend(buf,
312  nldofs,
314  gtopo.GetNeighborRank(nbs[i]),
315  40822 + gtopo.GetGroupMasterGroup(gr),
316  gtopo.GetComm(),
317  &requests[request_counter]);
318  request_counter++;
319  }
320  }
321  }
322  buf += nldofs;
323  }
324 
325  MPI_Waitall(request_counter, requests, statuses);
326 
327  // copy the received data from the buffer to ldata
328  buf = (T *)group_buf.GetData();
329  for (gr = 1; gr < group_ldof.Size(); gr++)
330  {
331  const int nldofs = group_ldof.RowSize(gr);
332 
333  // ignore groups without dofs
334  if (nldofs == 0)
335  {
336  continue;
337  }
338 
339  if (!gtopo.IAmMaster(gr)) // we are not the master
340  {
341  const int *ldofs = group_ldof.GetRow(gr);
342  for (i = 0; i < nldofs; i++)
343  {
344  ldata[ldofs[i]] = buf[i];
345  }
346  }
347  buf += nldofs;
348  }
349 }
350 
351 template <class T>
352 void GroupCommunicator::Reduce(T *ldata, void (*Op)(OpData<T>))
353 {
354  if (group_buf_size == 0)
355  {
356  return;
357  }
358 
359  int i, gr, request_counter = 0;
360  OpData<T> opd;
361 
362  group_buf.SetSize(group_buf_size*sizeof(T));
363  opd.ldata = ldata;
364  opd.buf = (T *)group_buf.GetData();
365  for (gr = 1; gr < group_ldof.Size(); gr++)
366  {
367  opd.nldofs = group_ldof.RowSize(gr);
368 
369  // ignore groups without dofs
370  if (opd.nldofs == 0)
371  {
372  continue;
373  }
374 
375  opd.ldofs = group_ldof.GetRow(gr);
376 
377  if (!gtopo.IAmMaster(gr)) // we are not the master
378  {
379  for (i = 0; i < opd.nldofs; i++)
380  {
381  opd.buf[i] = ldata[opd.ldofs[i]];
382  }
383 
384  MPI_Isend(opd.buf,
385  opd.nldofs,
387  gtopo.GetGroupMasterRank(gr),
388  43822 + gtopo.GetGroupMasterGroup(gr),
389  gtopo.GetComm(),
390  &requests[request_counter]);
391  request_counter++;
392  opd.buf += opd.nldofs;
393  }
394  else // we are the master
395  {
396  const int gs = gtopo.GetGroupSize(gr);
397  const int *nbs = gtopo.GetGroup(gr);
398  for (i = 0; i < gs; i++)
399  {
400  if (nbs[i] != 0)
401  {
402  MPI_Irecv(opd.buf,
403  opd.nldofs,
405  gtopo.GetNeighborRank(nbs[i]),
406  43822 + gtopo.GetGroupMasterGroup(gr),
407  gtopo.GetComm(),
408  &requests[request_counter]);
409  request_counter++;
410  opd.buf += opd.nldofs;
411  }
412  }
413  }
414  }
415 
416  MPI_Waitall(request_counter, requests, statuses);
417 
418  // perform the reduce operation
419  opd.buf = (T *)group_buf.GetData();
420  for (gr = 1; gr < group_ldof.Size(); gr++)
421  {
422  opd.nldofs = group_ldof.RowSize(gr);
423 
424  // ignore groups without dofs
425  if (opd.nldofs == 0)
426  {
427  continue;
428  }
429 
430  if (!gtopo.IAmMaster(gr)) // we are not the master
431  {
432  opd.buf += opd.nldofs;
433  }
434  else // we are the master
435  {
436  opd.ldofs = group_ldof.GetRow(gr);
437  opd.nb = gtopo.GetGroupSize(gr)-1;
438  Op(opd);
439  opd.buf += opd.nb * opd.nldofs;
440  }
441  }
442 }
443 
444 template <class T>
446 {
447  for (int i = 0; i < opd.nldofs; i++)
448  {
449  T data = opd.ldata[opd.ldofs[i]];
450  for (int j = 0; j < opd.nb; j++)
451  {
452  data += opd.buf[j*opd.nldofs+i];
453  }
454  opd.ldata[opd.ldofs[i]] = data;
455  }
456 }
457 
458 template <class T>
460 {
461  for (int i = 0; i < opd.nldofs; i++)
462  {
463  T data = opd.ldata[opd.ldofs[i]];
464  for (int j = 0; j < opd.nb; j++)
465  {
466  T b = opd.buf[j*opd.nldofs+i];
467  if (data > b)
468  {
469  data = b;
470  }
471  }
472  opd.ldata[opd.ldofs[i]] = data;
473  }
474 }
475 
476 template <class T>
478 {
479  for (int i = 0; i < opd.nldofs; i++)
480  {
481  T data = opd.ldata[opd.ldofs[i]];
482  for (int j = 0; j < opd.nb; j++)
483  {
484  T b = opd.buf[j*opd.nldofs+i];
485  if (data < b)
486  {
487  data = b;
488  }
489  }
490  opd.ldata[opd.ldofs[i]] = data;
491  }
492 }
493 
494 template <class T>
496 {
497  for (int i = 0; i < opd.nldofs; i++)
498  {
499  T data = opd.ldata[opd.ldofs[i]];
500  for (int j = 0; j < opd.nb; j++)
501  {
502  data |= opd.buf[j*opd.nldofs+i];
503  }
504  opd.ldata[opd.ldofs[i]] = data;
505  }
506 }
507 
509 {
510  delete [] statuses;
511  delete [] requests;
512 }
513 
514 // @cond DOXYGEN_SKIP
515 
516 // instantiate GroupCommunicator::Bcast and Reduce for int and double
517 template void GroupCommunicator::Bcast<int>(int *);
518 template void GroupCommunicator::Reduce<int>(int *, void (*)(OpData<int>));
519 
520 template void GroupCommunicator::Bcast<double>(double *);
521 template void GroupCommunicator::Reduce<double>(
522  double *, void (*)(OpData<double>));
523 
524 // @endcond
525 
526 // instantiate reduce operators for int and double
527 template void GroupCommunicator::Sum<int>(OpData<int>);
528 template void GroupCommunicator::Min<int>(OpData<int>);
529 template void GroupCommunicator::Max<int>(OpData<int>);
530 template void GroupCommunicator::BitOR<int>(OpData<int>);
531 
532 template void GroupCommunicator::Sum<double>(OpData<double>);
533 template void GroupCommunicator::Min<double>(OpData<double>);
534 template void GroupCommunicator::Max<double>(OpData<double>);
535 
536 
537 #ifdef __bgq__
538 static void DebugRankCoords(int** coords, int dim, int size)
539 {
540  for (int i = 0; i < size; i++)
541  {
542  std::cout << "Rank " << i << " coords: ";
543  for (int j = 0; j < dim; j++)
544  {
545  std::cout << coords[i][j] << " ";
546  }
547  std::cout << endl;
548  }
549 }
550 
551 struct CompareCoords
552 {
553  CompareCoords(int coord) : coord(coord) {}
554  int coord;
555 
556  bool operator()(int* const &a, int* const &b) const
557  { return a[coord] < b[coord]; }
558 };
559 
560 void KdTreeSort(int** coords, int d, int dim, int size)
561 {
562  if (size > 1)
563  {
564  bool all_same = true;
565  for (int i = 1; i < size && all_same; i++)
566  {
567  for (int j = 0; j < dim; j++)
568  {
569  if (coords[i][j] != coords[0][j]) { all_same = false; break; }
570  }
571  }
572  if (all_same) { return; }
573 
574  // sort by coordinate 'd'
575  std::sort(coords, coords + size, CompareCoords(d));
576  int next = (d + 1) % dim;
577 
578  if (coords[0][d] < coords[size-1][d])
579  {
580  KdTreeSort(coords, next, dim, size/2);
581  KdTreeSort(coords + size/2, next, dim, size - size/2);
582  }
583  else
584  {
585  // skip constant dimension
586  KdTreeSort(coords, next, dim, size);
587  }
588  }
589 }
590 
591 MPI_Comm ReorderRanksZCurve(MPI_Comm comm)
592 {
593  MPI_Status status;
594 
595  int rank, size;
596  MPI_Comm_rank(comm, &rank);
597  MPI_Comm_size(comm, &size);
598 
599  int dim;
600  MPIX_Torus_ndims(&dim);
601 
602  int* mycoords = new int[dim + 1];
603  MPIX_Rank2torus(rank, mycoords);
604 
605  MPI_Send(mycoords, dim, MPI_INT, 0, 111, comm);
606  delete [] mycoords;
607 
608  if (rank == 0)
609  {
610  int** coords = new int*[size];
611  for (int i = 0; i < size; i++)
612  {
613  coords[i] = new int[dim + 1];
614  coords[i][dim] = i;
615  MPI_Recv(coords[i], dim, MPI_INT, i, 111, comm, &status);
616  }
617 
618  KdTreeSort(coords, 0, dim, size);
619 
620  //DebugRankCoords(coords, dim, size);
621 
622  for (int i = 0; i < size; i++)
623  {
624  MPI_Send(&coords[i][dim], 1, MPI_INT, i, 112, comm);
625  delete [] coords[i];
626  }
627  delete [] coords;
628  }
629 
630  int new_rank;
631  MPI_Recv(&new_rank, 1, MPI_INT, 0, 112, comm, &status);
632 
633  MPI_Comm new_comm;
634  MPI_Comm_split(comm, 0, new_rank, &new_comm);
635  return new_comm;
636 }
637 
638 #else // __bgq__
639 
640 MPI_Comm ReorderRanksZCurve(MPI_Comm comm)
641 {
642  // pass
643  return comm;
644 }
645 #endif // __bgq__
646 
647 
648 } // namespace mfem
649 
650 #endif
int Lookup(IntegerSet &s)
Definition: sets.cpp:95
int GetGroupMasterRank(int g) const
void Create(ListOfIntegerSets &groups, int mpitag)
int Size() const
Logical size of the array.
Definition: array.hpp:109
void Recreate(const int n, const int *p)
Definition: sets.cpp:59
int GetGroupMasterGroup(int g) const
int * GetJ()
Definition: table.hpp:108
GroupCommunicator(GroupTopology &gt)
void AsTable(Table &t)
Definition: sets.cpp:107
Helper struct to convert a C++ type to an MPI type.
MPI_Comm ReorderRanksZCurve(MPI_Comm comm)
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
Definition: table.cpp:74
int GetGroupSize(int g) const
void SetDims(int rows, int nnz)
Definition: table.cpp:132
const int * GetGroup(int g) const
void Copy(Array &copy) const
Create a copy of the current array.
Definition: array.hpp:160
T * GetData()
Returns the data.
Definition: array.hpp:91
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:179
int Size_of_connections() const
Definition: table.hpp:92
int PickElementInSet(int i)
Definition: sets.hpp:58
bool IAmMaster(int g) const
int dim
Definition: ex3.cpp:47
static void Sum(OpData< T >)
Reduce operation Sum, instantiated for int and double.
void Finalize()
Allocate internal buffers after the GroupLDofTable is defined.
void AddConnection(int r, int c)
Definition: table.hpp:74
static void Min(OpData< T >)
Reduce operation Min, instantiated for int and double.
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:86
int GetNeighborRank(int i) const
static void Max(OpData< T >)
Reduce operation Max, instantiated for int and double.
void AddAColumnInRow(int r)
Definition: table.hpp:71
void mfem_error(const char *msg)
Definition: error.cpp:23
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:331
void ShiftUpI()
Definition: table.cpp:107
A set of integers.
Definition: sets.hpp:23
void MakeJ()
Definition: table.cpp:84
void Reduce(T *ldata, void(*Op)(OpData< T >))
void KdTreeSort(int **coords, int d, int dim, int size)
int NGroups() const
void Create(Array< int > &ldof_group)
int * GetI()
Definition: table.hpp:107
int RowSize(int i) const
Definition: table.hpp:102
List of integer sets.
Definition: sets.hpp:49
static void BitOR(OpData< T >)
Reduce operation bitwise OR, instantiated for int only.