MFEM  v3.3
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 "text.hpp"
26 
27 #include <iostream>
28 #include <map>
29 
30 using namespace std;
31 
32 namespace mfem
33 {
34 
35 GroupTopology::GroupTopology(const GroupTopology &gt)
36  : MyComm(gt.MyComm),
37  group_lproc(gt.group_lproc)
38 {
39  gt.groupmaster_lproc.Copy(groupmaster_lproc);
40  gt.lproc_proc.Copy(lproc_proc);
41  gt.group_mgroup.Copy(group_mgroup);
42 }
43 
44 void GroupTopology::ProcToLProc()
45 {
46  int NRanks;
47  MPI_Comm_size(MyComm, &NRanks);
48 
49  map<int, int> proc_lproc;
50 
51  int lproc_counter = 0;
52  for (int i = 0; i < group_lproc.Size_of_connections(); i++)
53  {
54  const pair<const int, int> p(group_lproc.GetJ()[i], lproc_counter);
55  if (proc_lproc.insert(p).second)
56  {
57  lproc_counter++;
58  }
59  }
60  // Note: group_lproc.GetJ()[0] == MyRank --> proc_lproc[MyRank] == 0
61 
62  lproc_proc.SetSize(lproc_counter);
63  for (map<int, int>::iterator it = proc_lproc.begin();
64  it != proc_lproc.end(); ++it)
65  {
66  lproc_proc[it->second] = it->first;
67  }
68 
69  for (int i = 0; i < group_lproc.Size_of_connections(); i++)
70  {
71  group_lproc.GetJ()[i] = proc_lproc[group_lproc.GetJ()[i]];
72  }
73 
74  for (int i = 0; i < NGroups(); i++)
75  {
76  groupmaster_lproc[i] = proc_lproc[groupmaster_lproc[i]];
77  }
78 }
79 
80 void GroupTopology::Create(ListOfIntegerSets &groups, int mpitag)
81 {
82  groups.AsTable(group_lproc); // group_lproc = group_proc
83 
84  Table group_mgroupandproc;
85  group_mgroupandproc.SetDims(NGroups(),
86  group_lproc.Size_of_connections() + NGroups());
87  for (int i = 0; i < NGroups(); i++)
88  {
89  int j = group_mgroupandproc.GetI()[i];
90  group_mgroupandproc.GetI()[i+1] = j + group_lproc.RowSize(i) + 1;
91  group_mgroupandproc.GetJ()[j] = i;
92  j++;
93  for (int k = group_lproc.GetI()[i];
94  j < group_mgroupandproc.GetI()[i+1]; j++, k++)
95  {
96  group_mgroupandproc.GetJ()[j] = group_lproc.GetJ()[k];
97  }
98  }
99 
100  // build groupmaster_lproc with lproc = proc
101  groupmaster_lproc.SetSize(NGroups());
102 
103  // simplest choice of the group owner
104  for (int i = 0; i < NGroups(); i++)
105  {
106  groupmaster_lproc[i] = groups.PickElementInSet(i);
107  }
108 
109  // load-balanced choice of the group owner, which however can lead to
110  // isolated dofs
111  // for (i = 0; i < NGroups(); i++)
112  // groupmaster_lproc[i] = groups.PickRandomElementInSet(i);
113 
114  ProcToLProc();
115 
116  // build group_mgroup
117  group_mgroup.SetSize(NGroups());
118  group_mgroup[0] = 0; // the local group
119 
120  int send_counter = 0;
121  int recv_counter = 0;
122  for (int i = 1; i < NGroups(); i++)
123  if (groupmaster_lproc[i] != 0) // we are not the master
124  {
125  recv_counter++;
126  }
127  else
128  {
129  send_counter += group_lproc.RowSize(i)-1;
130  }
131 
132  MPI_Request *requests = new MPI_Request[send_counter];
133  MPI_Status *statuses = new MPI_Status[send_counter];
134 
135  int max_recv_size = 0;
136  send_counter = 0;
137  for (int i = 1; i < NGroups(); i++)
138  {
139  if (groupmaster_lproc[i] == 0) // we are the master
140  {
141  group_mgroup[i] = i;
142 
143  for (int j = group_lproc.GetI()[i];
144  j < group_lproc.GetI()[i+1]; j++)
145  {
146  if (group_lproc.GetJ()[j] != 0)
147  {
148  MPI_Isend(group_mgroupandproc.GetRow(i),
149  group_mgroupandproc.RowSize(i),
150  MPI_INT,
151  lproc_proc[group_lproc.GetJ()[j]],
152  mpitag,
153  MyComm,
154  &requests[send_counter]);
155  send_counter++;
156  }
157  }
158  }
159  else // we are not the master
160  if (max_recv_size < group_lproc.RowSize(i))
161  {
162  max_recv_size = group_lproc.RowSize(i);
163  }
164  }
165  max_recv_size++;
166 
167  IntegerSet group;
168  if (recv_counter > 0)
169  {
170  int count;
171  MPI_Status status;
172  int *recv_buf = new int[max_recv_size];
173  for ( ; recv_counter > 0; recv_counter--)
174  {
175  MPI_Recv(recv_buf, max_recv_size, MPI_INT,
176  MPI_ANY_SOURCE, mpitag, MyComm, &status);
177 
178  MPI_Get_count(&status, MPI_INT, &count);
179 
180  group.Recreate(count-1, recv_buf+1);
181  int g = groups.Lookup(group);
182  group_mgroup[g] = recv_buf[0];
183 
184  if (lproc_proc[groupmaster_lproc[g]] != status.MPI_SOURCE)
185  {
186  cerr << "\n\n\nGroupTopology::GroupTopology: "
187  << MyRank() << ": ERROR\n\n\n" << endl;
188  mfem_error();
189  }
190  }
191  delete [] recv_buf;
192  }
193 
194  MPI_Waitall(send_counter, requests, statuses);
195 
196  delete [] statuses;
197  delete [] requests;
198 }
199 
200 void GroupTopology::Save(ostream &out) const
201 {
202  out << "\ncommunication_groups\n";
203  out << "number_of_groups " << NGroups() << "\n\n";
204 
205  out << "# number of entities in each group, followed by group ids in group\n";
206  for (int group_id = 0; group_id < NGroups(); ++group_id)
207  {
208  int group_size = GetGroupSize(group_id);
209  const int * group_ptr = GetGroup(group_id);
210  out << group_size;
211  for ( int group_member_index = 0; group_member_index < group_size;
212  ++group_member_index)
213  {
214  out << " " << GetNeighborRank( group_ptr[group_member_index] );
215  }
216  out << "\n";
217  }
218 
219  // For future use, optional ownership strategy.
220  // out << "# ownership";
221 }
222 
223 void GroupTopology::Load(istream &in)
224 {
225  // Load in group topology and create list of integer sets. Use constructor
226  // that uses list of integer sets.
227  std::string ident;
228 
229  // Read in number of groups
230  int number_of_groups = -1;
231  in >> ident;
232  MFEM_VERIFY(ident == "number_of_groups",
233  "GroupTopology::Load - expected 'number_of_groups' entry.");
234  in >> number_of_groups;
235 
236  // Skip number of entries in each group comment.
237  skip_comment_lines(in, '#');
238 
239  ListOfIntegerSets integer_sets;
240  for (int group_id = 0; group_id < number_of_groups; ++group_id)
241  {
242  IntegerSet integer_set;
243  Array<int>& array = integer_set;
244  int group_size;
245  in >> group_size;
246  array.Reserve(group_size);
247  for ( int index = 0; index < group_size; ++index )
248  {
249  int value;
250  in >> value;
251  array.Append(value);
252  }
253  integer_sets.Insert(integer_set);
254  }
255 
256  Create(integer_sets, 823);
257 }
258 
259 
260 // specializations of the MPITypeMap static member
261 template<> const MPI_Datatype MPITypeMap<int>::mpi_type = MPI_INT;
262 template<> const MPI_Datatype MPITypeMap<double>::mpi_type = MPI_DOUBLE;
263 
264 
266  : gtopo(gt)
267 {
268  group_buf_size = 0;
269  requests = NULL;
270  statuses = NULL;
271 }
272 
274 {
275  group_ldof.MakeI(gtopo.NGroups());
276  for (int i = 0; i < ldof_group.Size(); i++)
277  {
278  int group = ldof_group[i];
279  if (group != 0)
280  {
281  group_ldof.AddAColumnInRow(group);
282  }
283  }
284  group_ldof.MakeJ();
285 
286  for (int i = 0; i < ldof_group.Size(); i++)
287  {
288  int group = ldof_group[i];
289  if (group != 0)
290  {
291  group_ldof.AddConnection(group, i);
292  }
293  }
294  group_ldof.ShiftUpI();
295 
296  Finalize();
297 }
298 
300 {
301  int request_counter = 0;
302 
303  for (int gr = 1; gr < group_ldof.Size(); gr++)
304  if (group_ldof.RowSize(gr) != 0)
305  {
306  int gr_requests;
307  if (!gtopo.IAmMaster(gr)) // we are not the master
308  {
309  gr_requests = 1;
310  }
311  else
312  {
313  gr_requests = gtopo.GetGroupSize(gr)-1;
314  }
315 
316  request_counter += gr_requests;
317  group_buf_size += gr_requests * group_ldof.RowSize(gr);
318  }
319 
320  requests = new MPI_Request[request_counter];
321  statuses = new MPI_Status[request_counter];
322 }
323 
324 template <class T>
325 void GroupCommunicator::Bcast(T *ldata, int layout)
326 {
327  if (group_buf_size == 0) { return; }
328 
329  T *buf;
330  if (layout == 0)
331  {
332  group_buf.SetSize(group_buf_size*sizeof(T));
333  buf = (T *)group_buf.GetData();
334  }
335  else
336  {
337  buf = ldata;
338  }
339 
340  int i, gr, request_counter = 0;
341 
342  for (gr = 1; gr < group_ldof.Size(); gr++)
343  {
344  const int nldofs = group_ldof.RowSize(gr);
345 
346  // ignore groups without dofs
347  if (nldofs == 0) { continue; }
348 
349  if (!gtopo.IAmMaster(gr)) // we are not the master
350  {
351  MPI_Irecv(buf,
352  nldofs,
354  gtopo.GetGroupMasterRank(gr),
355  40822 + gtopo.GetGroupMasterGroup(gr),
356  gtopo.GetComm(),
357  &requests[request_counter]);
358  request_counter++;
359  }
360  else // we are the master
361  {
362  if (layout == 0)
363  {
364  // fill send buffer
365  const int *ldofs = group_ldof.GetRow(gr);
366  for (i = 0; i < nldofs; i++)
367  {
368  buf[i] = ldata[ldofs[i]];
369  }
370  }
371 
372  const int gs = gtopo.GetGroupSize(gr);
373  const int *nbs = gtopo.GetGroup(gr);
374  for (i = 0; i < gs; i++)
375  {
376  if (nbs[i] != 0)
377  {
378  MPI_Isend(buf,
379  nldofs,
381  gtopo.GetNeighborRank(nbs[i]),
382  40822 + gtopo.GetGroupMasterGroup(gr),
383  gtopo.GetComm(),
384  &requests[request_counter]);
385  request_counter++;
386  }
387  }
388  }
389  buf += nldofs;
390  }
391 
392  MPI_Waitall(request_counter, requests, statuses);
393 
394  if (layout == 0)
395  {
396  // copy the received data from the buffer to ldata
397  buf = (T *)group_buf.GetData();
398  for (gr = 1; gr < group_ldof.Size(); gr++)
399  {
400  const int nldofs = group_ldof.RowSize(gr);
401 
402  // ignore groups without dofs
403  if (nldofs == 0) { continue; }
404 
405  if (!gtopo.IAmMaster(gr)) // we are not the master
406  {
407  const int *ldofs = group_ldof.GetRow(gr);
408  for (i = 0; i < nldofs; i++)
409  {
410  ldata[ldofs[i]] = buf[i];
411  }
412  }
413  buf += nldofs;
414  }
415  }
416 }
417 
418 template <class T>
419 void GroupCommunicator::Reduce(T *ldata, void (*Op)(OpData<T>))
420 {
421  if (group_buf_size == 0) { return; }
422 
423  int i, gr, request_counter = 0;
424  OpData<T> opd;
425 
426  group_buf.SetSize(group_buf_size*sizeof(T));
427  opd.ldata = ldata;
428  opd.buf = (T *)group_buf.GetData();
429  for (gr = 1; gr < group_ldof.Size(); gr++)
430  {
431  opd.nldofs = group_ldof.RowSize(gr);
432 
433  // ignore groups without dofs
434  if (opd.nldofs == 0) { continue; }
435 
436  opd.ldofs = group_ldof.GetRow(gr);
437 
438  if (!gtopo.IAmMaster(gr)) // we are not the master
439  {
440  for (i = 0; i < opd.nldofs; i++)
441  {
442  opd.buf[i] = ldata[opd.ldofs[i]];
443  }
444 
445  MPI_Isend(opd.buf,
446  opd.nldofs,
448  gtopo.GetGroupMasterRank(gr),
449  43822 + gtopo.GetGroupMasterGroup(gr),
450  gtopo.GetComm(),
451  &requests[request_counter]);
452  request_counter++;
453  opd.buf += opd.nldofs;
454  }
455  else // we are the master
456  {
457  const int gs = gtopo.GetGroupSize(gr);
458  const int *nbs = gtopo.GetGroup(gr);
459  for (i = 0; i < gs; i++)
460  {
461  if (nbs[i] != 0)
462  {
463  MPI_Irecv(opd.buf,
464  opd.nldofs,
466  gtopo.GetNeighborRank(nbs[i]),
467  43822 + gtopo.GetGroupMasterGroup(gr),
468  gtopo.GetComm(),
469  &requests[request_counter]);
470  request_counter++;
471  opd.buf += opd.nldofs;
472  }
473  }
474  }
475  }
476 
477  MPI_Waitall(request_counter, requests, statuses);
478 
479  // perform the reduce operation
480  opd.buf = (T *)group_buf.GetData();
481  for (gr = 1; gr < group_ldof.Size(); gr++)
482  {
483  opd.nldofs = group_ldof.RowSize(gr);
484 
485  // ignore groups without dofs
486  if (opd.nldofs == 0) { continue; }
487 
488  if (!gtopo.IAmMaster(gr)) // we are not the master
489  {
490  opd.buf += opd.nldofs;
491  }
492  else // we are the master
493  {
494  opd.ldofs = group_ldof.GetRow(gr);
495  opd.nb = gtopo.GetGroupSize(gr)-1;
496  Op(opd);
497  opd.buf += opd.nb * opd.nldofs;
498  }
499  }
500 }
501 
502 template <class T>
504 {
505  for (int i = 0; i < opd.nldofs; i++)
506  {
507  T data = opd.ldata[opd.ldofs[i]];
508  for (int j = 0; j < opd.nb; j++)
509  {
510  data += opd.buf[j*opd.nldofs+i];
511  }
512  opd.ldata[opd.ldofs[i]] = data;
513  }
514 }
515 
516 template <class T>
518 {
519  for (int i = 0; i < opd.nldofs; i++)
520  {
521  T data = opd.ldata[opd.ldofs[i]];
522  for (int j = 0; j < opd.nb; j++)
523  {
524  T b = opd.buf[j*opd.nldofs+i];
525  if (data > b)
526  {
527  data = b;
528  }
529  }
530  opd.ldata[opd.ldofs[i]] = data;
531  }
532 }
533 
534 template <class T>
536 {
537  for (int i = 0; i < opd.nldofs; i++)
538  {
539  T data = opd.ldata[opd.ldofs[i]];
540  for (int j = 0; j < opd.nb; j++)
541  {
542  T b = opd.buf[j*opd.nldofs+i];
543  if (data < b)
544  {
545  data = b;
546  }
547  }
548  opd.ldata[opd.ldofs[i]] = data;
549  }
550 }
551 
552 template <class T>
554 {
555  for (int i = 0; i < opd.nldofs; i++)
556  {
557  T data = opd.ldata[opd.ldofs[i]];
558  for (int j = 0; j < opd.nb; j++)
559  {
560  data |= opd.buf[j*opd.nldofs+i];
561  }
562  opd.ldata[opd.ldofs[i]] = data;
563  }
564 }
565 
567 {
568  delete [] statuses;
569  delete [] requests;
570 }
571 
572 // @cond DOXYGEN_SKIP
573 
574 // instantiate GroupCommunicator::Bcast and Reduce for int and double
575 template void GroupCommunicator::Bcast<int>(int *);
576 template void GroupCommunicator::Bcast<int>(int *, int);
577 template void GroupCommunicator::Reduce<int>(int *, void (*)(OpData<int>));
578 
579 template void GroupCommunicator::Bcast<double>(double *);
580 template void GroupCommunicator::Bcast<double>(double *, int);
581 template void GroupCommunicator::Reduce<double>(
582  double *, void (*)(OpData<double>));
583 
584 // @endcond
585 
586 // instantiate reduce operators for int and double
587 template void GroupCommunicator::Sum<int>(OpData<int>);
588 template void GroupCommunicator::Min<int>(OpData<int>);
589 template void GroupCommunicator::Max<int>(OpData<int>);
590 template void GroupCommunicator::BitOR<int>(OpData<int>);
591 
592 template void GroupCommunicator::Sum<double>(OpData<double>);
593 template void GroupCommunicator::Min<double>(OpData<double>);
594 template void GroupCommunicator::Max<double>(OpData<double>);
595 
596 
597 #ifdef __bgq__
598 static void DebugRankCoords(int** coords, int dim, int size)
599 {
600  for (int i = 0; i < size; i++)
601  {
602  std::cout << "Rank " << i << " coords: ";
603  for (int j = 0; j < dim; j++)
604  {
605  std::cout << coords[i][j] << " ";
606  }
607  std::cout << endl;
608  }
609 }
610 
611 struct CompareCoords
612 {
613  CompareCoords(int coord) : coord(coord) {}
614  int coord;
615 
616  bool operator()(int* const &a, int* const &b) const
617  { return a[coord] < b[coord]; }
618 };
619 
620 void KdTreeSort(int** coords, int d, int dim, int size)
621 {
622  if (size > 1)
623  {
624  bool all_same = true;
625  for (int i = 1; i < size && all_same; i++)
626  {
627  for (int j = 0; j < dim; j++)
628  {
629  if (coords[i][j] != coords[0][j]) { all_same = false; break; }
630  }
631  }
632  if (all_same) { return; }
633 
634  // sort by coordinate 'd'
635  std::sort(coords, coords + size, CompareCoords(d));
636  int next = (d + 1) % dim;
637 
638  if (coords[0][d] < coords[size-1][d])
639  {
640  KdTreeSort(coords, next, dim, size/2);
641  KdTreeSort(coords + size/2, next, dim, size - size/2);
642  }
643  else
644  {
645  // skip constant dimension
646  KdTreeSort(coords, next, dim, size);
647  }
648  }
649 }
650 
651 MPI_Comm ReorderRanksZCurve(MPI_Comm comm)
652 {
653  MPI_Status status;
654 
655  int rank, size;
656  MPI_Comm_rank(comm, &rank);
657  MPI_Comm_size(comm, &size);
658 
659  int dim;
660  MPIX_Torus_ndims(&dim);
661 
662  int* mycoords = new int[dim + 1];
663  MPIX_Rank2torus(rank, mycoords);
664 
665  MPI_Send(mycoords, dim, MPI_INT, 0, 111, comm);
666  delete [] mycoords;
667 
668  if (rank == 0)
669  {
670  int** coords = new int*[size];
671  for (int i = 0; i < size; i++)
672  {
673  coords[i] = new int[dim + 1];
674  coords[i][dim] = i;
675  MPI_Recv(coords[i], dim, MPI_INT, i, 111, comm, &status);
676  }
677 
678  KdTreeSort(coords, 0, dim, size);
679 
680  //DebugRankCoords(coords, dim, size);
681 
682  for (int i = 0; i < size; i++)
683  {
684  MPI_Send(&coords[i][dim], 1, MPI_INT, i, 112, comm);
685  delete [] coords[i];
686  }
687  delete [] coords;
688  }
689 
690  int new_rank;
691  MPI_Recv(&new_rank, 1, MPI_INT, 0, 112, comm, &status);
692 
693  MPI_Comm new_comm;
694  MPI_Comm_split(comm, 0, new_rank, &new_comm);
695  return new_comm;
696 }
697 
698 #else // __bgq__
699 
700 MPI_Comm ReorderRanksZCurve(MPI_Comm comm)
701 {
702  // pass
703  return comm;
704 }
705 #endif // __bgq__
706 
707 
708 } // namespace mfem
709 
710 #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
void skip_comment_lines(std::istream &is, const char comment_char)
Definition: text.hpp:25
int PickElementInSet(int i)
Definition: sets.hpp:63
bool IAmMaster(int g) const
int dim
Definition: ex3.cpp:47
static void Sum(OpData< T >)
Reduce operation Sum, instantiated for int and double.
int Append(const T &el)
Append element to array, resize if necessary.
Definition: array.hpp:394
void Finalize()
Allocate internal buffers after the GroupLDofTable is defined.
void AddConnection(int r, int c)
Definition: table.hpp:74
int Insert(IntegerSet &s)
Definition: sets.cpp:82
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:122
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
void Save(std::ostream &out) const
Save the data in a stream.
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:106
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:349
void ShiftUpI()
Definition: table.cpp:107
Data structure on which we define reduce operations.
A set of integers.
Definition: sets.hpp:23
void MakeJ()
Definition: table.cpp:84
void Reduce(T *ldata, void(*Op)(OpData< T >))
Reduce within each group where the master is the root.
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:54
void Bcast(T *data, int layout)
Broadcast within each group where the master is the root. The data layout can be: ...
void Load(std::istream &in)
Load the data from a stream.
static void BitOR(OpData< T >)
Reduce operation bitwise OR, instantiated for int only.