dune-grid-glue  2.8.0
ringcomm.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 /* IMPLEMENTATION OF CLASS G R I D G L U E */
4 
6 #define CheckMPIStatus(A,B) {}
7 
8 #include <mpi.h>
9 #include <functional>
10 #include <utility>
11 
12 #include <dune/common/fvector.hh>
13 #include <dune/common/hybridutilities.hh>
14 
15 #include <dune/geometry/type.hh>
16 
17 namespace Dune {
18 namespace Parallel {
19 
20  namespace Impl {
21 
23  template<typename T>
24  struct MPITypeInfo {};
25 
26  template<>
27  struct MPITypeInfo< int >
28  {
29  static const unsigned int size = 1;
30  static inline MPI_Datatype getType()
31  {
32  return MPI_INT;
33  }
34  };
35 
36  template<typename K, int N>
37  struct MPITypeInfo< Dune::FieldVector<K,N> >
38  {
39  static const unsigned int size = N;
40  static inline MPI_Datatype getType()
41  {
42  return Dune::MPITraits<K>::getType();
43  }
44  };
45 
46  template<>
47  struct MPITypeInfo< unsigned int >
48  {
49  static const unsigned int size = 1;
50  static inline MPI_Datatype getType()
51  {
52  return MPI_UNSIGNED;
53  }
54  };
55 
56  template<>
57  struct MPITypeInfo< Dune::GeometryType >
58  {
59  static const unsigned int size = 1;
60  static inline MPI_Datatype getType()
61  {
62  return Dune::MPITraits< Dune::GeometryType >::getType();
63  }
64  };
65 
66  template<typename T>
67  void MPI_SetVectorSize(
68  std::vector<T> & data,
69  MPI_Status & status)
70  {
71  typedef MPITypeInfo<T> Info;
72  int sz;
73  MPI_Get_count(&status, Info::getType(), &sz);
74  assert(sz%Info::size == 0);
75  data.resize(sz/Info::size);
76  }
77 
87  template<typename T>
88  void MPI_SendVectorInRing(
89  std::vector<T> & data,
90  std::vector<T> & next,
91  int tag,
92  int rightrank,
93  int leftrank,
94  MPI_Comm comm,
95  MPI_Request& r_send,
96  MPI_Request& r_recv
97  )
98  {
99  // mpi status stuff
100  [[maybe_unused]] int result = 0;
101  typedef MPITypeInfo<T> Info;
102  // resize next buffer to maximum size
103  next.resize(next.capacity());
104  // send data (explicitly send data.size elements)
105  result =
106  MPI_Isend(
107  &(data[0]), Info::size*data.size(), Info::getType(), rightrank, tag,
108  comm, &r_send);
109  // receive up to maximum size. The acutal size is stored in the status
110  result =
111  MPI_Irecv(
112  &(next[0]), Info::size*next.size(), Info::getType(), leftrank, tag,
113  comm, &r_recv);
114  // // check result
115  // MPI_Status status;
116  // CheckMPIStatus(result, status);
117  }
118 
119  template<typename T>
120  using ptr_t = T*;
121 
122  /* these helper structs are needed as long as we still support
123  C++11, as we can't use variadic lambdas */
124  template<typename... Args>
125  struct call_MPI_SendVectorInRing
126  {
127  std::tuple<Args...> & remotedata;
128  std::tuple<Args...> & nextdata;
129  int & tag;
130  int & rightrank;
131  int & leftrank;
132  MPI_Comm & mpicomm;
133  std::array<MPI_Request,sizeof...(Args)> & requests_send;
134  std::array<MPI_Request,sizeof...(Args)> & requests_recv;
135 
136  template<typename I>
137  void operator()(I i)
138  {
139  MPI_SendVectorInRing(
140  std::get<i>(remotedata),
141  std::get<i>(nextdata),
142  tag+i,
143  rightrank, leftrank, mpicomm,
144  requests_send[i],
145  requests_recv[i]);
146  }
147  };
148  template<typename... Args>
149  struct call_MPI_SetVectorSize
150  {
151  std::tuple<Args...> & nextdata;
152  std::array<MPI_Status,sizeof...(Args)> & status_recv;
153 
154  template<typename I>
155  void operator()(I i)
156  {
157  MPI_SetVectorSize(std::get<i>(nextdata),status_recv[i]);
158  }
159  };
160 
161  template<typename OP, std::size_t... Indices, typename... Args>
162  void MPI_AllApply_impl(MPI_Comm mpicomm,
163  OP && op,
164  std::index_sequence<Indices...> indices,
165  const Args&... data)
166  {
167  constexpr std::size_t N = sizeof...(Args);
168  int myrank = 0;
169  int commsize = 0;
170 #if HAVE_MPI
171  MPI_Comm_rank(mpicomm, &myrank);
172  MPI_Comm_size(mpicomm, &commsize);
173 #endif // HAVE_MPI
174 
175  if (commsize > 1)
176  {
177 #ifdef DEBUG_GRIDGLUE_PARALLELMERGE
178  std::cout << myrank << " Start Communication, size " << commsize << std::endl;
179 #endif
180 
181  // get data sizes
182  std::array<unsigned int, N> size({ ((unsigned int)data.size())... });
183 
184  // communicate max data size
185  std::array<unsigned int, N> maxSize;
186  MPI_Allreduce(&size, &maxSize,
187  size.size(), MPI_UNSIGNED, MPI_MAX, mpicomm);
188 #ifdef DEBUG_GRIDGLUE_PARALLELMERGE
189  std::cout << myrank << " maxSize " << "done... " << std::endl;
190 #endif
191 
192  // allocate receiving buffers with maxsize to ensure sufficient buffer size for communication
193  std::tuple<Args...> remotedata { Args(maxSize[Indices])... };
194 
195  // copy local data to receiving buffer
196  remotedata = std::tie(data...);
197 
198  // allocate second set of receiving buffers necessary for async communication
199  std::tuple<Args...> nextdata { Args(maxSize[Indices])... };
200 
201  // communicate data in the ring
202  int rightrank = (myrank + 1 + commsize) % commsize;
203  int leftrank = (myrank - 1 + commsize) % commsize;
204 
205  std::cout << myrank << ": size = " << commsize << std::endl;
206  std::cout << myrank << ": left = " << leftrank
207  << " right = " << rightrank << std::endl;
208 
209  // currently the remote data is our own data
210  int remoterank = myrank;
211 
212  for (int i=1; i<commsize; i++)
213  {
214  // in this iteration we will receive data from nextrank
215  int nextrank = (myrank - i + commsize) % commsize;
216 
217  std::cout << myrank << ": next = " << nextrank << std::endl;
218 
219  // send remote data to right neighbor and receive from left neighbor
220  std::array<MPI_Request,N> requests_send;
221  std::array<MPI_Request,N> requests_recv;
222 
223  int tag = 0;
224  Dune::Hybrid::forEach(indices,
225  // [&](auto i){
226  // MPI_SendVectorInRing(
227  // std::get<i>(remotedata),
228  // std::get<i>(nextdata),
229  // tag+i,
230  // rightrank, leftrank, mpicomm,
231  // requests_send[i],
232  // requests_recv[i]);
233  // });
234  call_MPI_SendVectorInRing<Args...>({
235  remotedata,
236  nextdata,
237  tag,
238  rightrank, leftrank, mpicomm,
239  requests_send,
240  requests_recv
241  }));
242 
243  // apply operator
244  op(remoterank,std::get<Indices>(remotedata)...);
245 
246  // wait for communication to finalize
247  std::array<MPI_Status,N> status_send;
248  std::array<MPI_Status,N> status_recv;
249  MPI_Waitall(N,&requests_recv[0],&status_recv[0]);
250 
251  // we finished receiving from nextrank and thus remoterank = nextrank
252  remoterank = nextrank;
253 
254  // get current data sizes
255  // and resize vectors
256  Dune::Hybrid::forEach(indices,
257  // [&](auto i){
258  // MPI_SetVectorSize(std::get<i>(nextdata),status_recv[i]);
259  // });
260  call_MPI_SetVectorSize<Args...>({
261  nextdata, status_recv
262  }));
263 
264  MPI_Waitall(N,&requests_send[0],&status_send[0]);
265 
266  // swap the communication buffers
267  std::swap(remotedata,nextdata);
268  }
269 
270  // last apply (or the only one in the case of sequential application)
271  op(remoterank,std::get<Indices>(remotedata)...);
272  }
273  else // sequential
274  {
275  op(myrank,data...);
276  }
277  }
278 
279  } // end namespace Impl
280 
294 template<typename OP, typename... Args>
295 void MPI_AllApply(MPI_Comm mpicomm,
296  OP && op,
297  const Args& ... data)
298 {
299  Impl::MPI_AllApply_impl(
300  mpicomm,
301  std::forward<OP>(op),
302  std::make_index_sequence<sizeof...(Args)>(),
303  data...
304  );
305 }
306 
307 } // end namespace Parallel
308 } // end namespace Dune
Definition: gridglue.hh:35
void MPI_AllApply(MPI_Comm mpicomm, OP &&op, const Args &... data)
apply an operator locally to a difstributed data set
Definition: ringcomm.hh:295