dune-grid  2.4.1-rc2
torus.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 #ifndef DUNE_GRID_YASPGRID_TORUS_HH
4 #define DUNE_GRID_YASPGRID_TORUS_HH
5 
6 #include <bitset>
7 #include <cmath>
8 #include <deque>
9 #include <vector>
10 
11 #if HAVE_MPI
12 #include <mpi.h>
13 #endif
14 
15 #include <dune/common/array.hh>
16 #include <dune/common/binaryfunctions.hh>
18 
19 #include "partitioning.hh"
20 
25 namespace Dune
26 {
27 
41  template<class CollectiveCommunication, int d>
42  class Torus {
43  public:
45  typedef Dune::array<int, d> iTupel;
46 
47 
48  private:
49  struct CommPartner {
50  int rank;
51  iTupel delta;
52  int index;
53  };
54 
55  struct CommTask {
56  int rank; // process to send to / receive from
57  void *buffer; // buffer to send / receive
58  int size; // size of buffer
59 #if HAVE_MPI
60  MPI_Request request; // used by MPI to handle request
61 #else
62  int request;
63 #endif
64  int flag; // used by MPI
65  };
66 
67  public:
69  Torus ()
70  {}
71 
73  Torus (CollectiveCommunication comm, int tag, iTupel size, const YLoadBalance<d>* lb)
74  : _comm(comm), _tag(tag)
75  {
76  // determine dimensions
77  lb->loadbalance(size, _comm.size(), _dims);
78 
79  // compute increments for lexicographic ordering
80  int inc = 1;
81  for (int i=0; i<d; i++)
82  {
83  _increment[i] = inc;
84  inc *= _dims[i];
85  }
86 
87  // check whether the load balancing matches the size of the communicator
88  if (inc != _comm.size())
89  DUNE_THROW(Dune::Exception, "Communicator size and result of the given load balancer do not match!");
90 
91  // make full schedule
92  proclists();
93  }
94 
96  int rank () const
97  {
98  return _comm.rank();
99  }
100 
102  iTupel coord () const
103  {
104  return rank_to_coord(_comm.rank());
105  }
106 
108  int procs () const
109  {
110  return _comm.size();
111  }
112 
114  const iTupel & dims () const
115  {
116  return _dims;
117  }
118 
120  int dims (int i) const
121  {
122  return _dims[i];
123  }
124 
126  CollectiveCommunication comm () const
127  {
128  return _comm;
129  }
130 
132  int tag () const
133  {
134  return _tag;
135  }
136 
138  bool inside (iTupel c) const
139  {
140  for (int i=d-1; i>=0; i--)
141  if (c[i]<0 || c[i]>=_dims[i]) return false;
142  return true;
143  }
144 
146  iTupel rank_to_coord (int rank) const
147  {
148  iTupel coord;
149  rank = rank%_comm.size();
150  for (int i=d-1; i>=0; i--)
151  {
152  coord[i] = rank/_increment[i];
153  rank = rank%_increment[i];
154  }
155  return coord;
156  }
157 
159  int coord_to_rank (iTupel coord) const
160  {
161  for (int i=0; i<d; i++) coord[i] = coord[i]%_dims[i];
162  int rank = 0;
163  for (int i=0; i<d; i++) rank += coord[i]*_increment[i];
164  return rank;
165  }
166 
168  int rank_relative (int rank, int dir, int cnt) const
169  {
170  iTupel coord = rank_to_coord(rank);
171  coord[dir] = (coord[dir]+_dims[dir]+cnt)%_dims[dir];
172  return coord_to_rank(coord);
173  }
174 
176  int color (const iTupel & coord) const
177  {
178  int c = 0;
179  int power = 1;
180 
181  // interior coloring
182  for (int i=0; i<d; i++)
183  {
184  if (coord[i]%2==1) c += power;
185  power *= 2;
186  }
187 
188  // extra colors for boundary processes
189  for (int i=0; i<d; i++)
190  {
191  if (_dims[i]>1 && coord[i]==_dims[i]-1) c += power;
192  power *= 2;
193  }
194 
195  return c;
196  }
197 
199  int color (int rank) const
200  {
201  return color(rank_to_coord(rank));
202  }
203 
205  int neighbors () const
206  {
207  int n=1;
208  for (int i=0; i<d; ++i)
209  n *= 3;
210  return n-1;
211  }
212 
214  bool is_neighbor (iTupel delta, std::bitset<d> periodic) const
215  {
216  iTupel coord = rank_to_coord(_comm.rank()); // my own coordinate with 0 <= c_i < dims_i
217 
218  for (int i=0; i<d; i++)
219  {
220  if (delta[i]<0)
221  {
222  // if I am on the boundary and domain is not periodic => no neighbor
223  if (coord[i]==0 && periodic[i]==false) return false;
224  }
225  if (delta[i]>0)
226  {
227  // if I am on the boundary and domain is not periodic => no neighbor
228  if (coord[i]==_dims[i]-1 && periodic[i]==false) return false;
229  }
230  }
231  return true;
232  }
233 
241  double partition (int rank, iTupel origin_in, iTupel size_in, iTupel& origin_out, iTupel& size_out) const
242  {
243  iTupel coord = rank_to_coord(rank);
244  double maxsize = 1;
245  double sz = 1;
246 
247  // make a tensor product partition
248  for (int i=0; i<d; i++)
249  {
250  // determine
251  int m = size_in[i]/_dims[i];
252  int r = size_in[i]%_dims[i];
253 
254  sz *= size_in[i];
255 
256  if (coord[i]<_dims[i]-r)
257  {
258  origin_out[i] = origin_in[i] + coord[i]*m;
259  size_out[i] = m;
260  maxsize *= m;
261  }
262  else
263  {
264  origin_out[i] = origin_in[i] + (_dims[i]-r)*m + (coord[i]-(_dims[i]-r))*(m+1);
265  size_out[i] = m+1;
266  maxsize *= m+1;
267  }
268  }
269  return maxsize/(sz/_comm.size());
270  }
271 
279  public:
281  ProcListIterator (typename std::deque<CommPartner>::const_iterator iter)
282  {
283  i = iter;
284  }
285 
287  int rank () const
288  {
289  return i->rank;
290  }
291 
293  iTupel delta () const
294  {
295  return i->delta;
296  }
297 
299  int index () const
300  {
301  return i->index;
302  }
303 
305  int distance () const
306  {
307  int dist = 0;
308  iTupel delta=i->delta;
309  for (int j=0; j<d; ++j)
310  dist += std::abs(delta[j]);
311  return dist;
312  }
313 
315  bool operator== (const ProcListIterator& iter)
316  {
317  return i == iter.i;
318  }
319 
320 
322  bool operator!= (const ProcListIterator& iter)
323  {
324  return i != iter.i;
325  }
326 
328  ProcListIterator& operator++ ()
329  {
330  ++i;
331  return *this;
332  }
333 
334  private:
335  typename std::deque<CommPartner>::const_iterator i;
336  };
337 
339  ProcListIterator sendbegin () const
340  {
341  return ProcListIterator(_sendlist.begin());
342  }
343 
345  ProcListIterator sendend () const
346  {
347  return ProcListIterator(_sendlist.end());
348  }
349 
351  ProcListIterator recvbegin () const
352  {
353  return ProcListIterator(_recvlist.begin());
354  }
355 
357  ProcListIterator recvend () const
358  {
359  return ProcListIterator(_recvlist.end());
360  }
361 
363  void send (int rank, void* buffer, int size) const
364  {
365  CommTask task;
366  task.rank = rank;
367  task.buffer = buffer;
368  task.size = size;
369  if (rank!=_comm.rank())
370  _sendrequests.push_back(task);
371  else
372  _localsendrequests.push_back(task);
373  }
374 
376  void recv (int rank, void* buffer, int size) const
377  {
378  CommTask task;
379  task.rank = rank;
380  task.buffer = buffer;
381  task.size = size;
382  if (rank!=_comm.rank())
383  _recvrequests.push_back(task);
384  else
385  _localrecvrequests.push_back(task);
386  }
387 
389  void exchange () const
390  {
391  // handle local requests first
392  if (_localsendrequests.size()!=_localrecvrequests.size())
393  {
394  std::cout << "[" << rank() << "]: ERROR: local sends/receives do not match in exchange!" << std::endl;
395  return;
396  }
397  for (unsigned int i=0; i<_localsendrequests.size(); i++)
398  {
399  if (_localsendrequests[i].size!=_localrecvrequests[i].size)
400  {
401  std::cout << "[" << rank() << "]: ERROR: size in local sends/receive does not match in exchange!" << std::endl;
402  return;
403  }
404  memcpy(_localrecvrequests[i].buffer,_localsendrequests[i].buffer,_localsendrequests[i].size);
405  }
406  _localsendrequests.clear();
407  _localrecvrequests.clear();
408 
409 #if HAVE_MPI
410  // handle foreign requests
411  int sends=0;
412  int recvs=0;
413 
414  // issue sends to foreign processes
415  for (unsigned int i=0; i<_sendrequests.size(); i++)
416  if (_sendrequests[i].rank!=rank())
417  {
418  // std::cout << "[" << rank() << "]" << " send " << _sendrequests[i].size << " bytes "
419  // << "to " << _sendrequests[i].rank << " p=" << _sendrequests[i].buffer << std::endl;
420  MPI_Isend(_sendrequests[i].buffer, _sendrequests[i].size, MPI_BYTE,
421  _sendrequests[i].rank, _tag, _comm, &(_sendrequests[i].request));
422  _sendrequests[i].flag = false;
423  sends++;
424  }
425 
426  // issue receives from foreign processes
427  for (unsigned int i=0; i<_recvrequests.size(); i++)
428  if (_recvrequests[i].rank!=rank())
429  {
430  // std::cout << "[" << rank() << "]" << " recv " << _recvrequests[i].size << " bytes "
431  // << "fm " << _recvrequests[i].rank << " p=" << _recvrequests[i].buffer << std::endl;
432  MPI_Irecv(_recvrequests[i].buffer, _recvrequests[i].size, MPI_BYTE,
433  _recvrequests[i].rank, _tag, _comm, &(_recvrequests[i].request));
434  _recvrequests[i].flag = false;
435  recvs++;
436  }
437 
438  // poll sends
439  while (sends>0)
440  {
441  for (unsigned int i=0; i<_sendrequests.size(); i++)
442  if (!_sendrequests[i].flag)
443  {
444  MPI_Status status;
445  MPI_Test( &(_sendrequests[i].request), &(_sendrequests[i].flag), &status);
446  if (_sendrequests[i].flag)
447  {
448  sends--;
449  // std::cout << "[" << rank() << "]" << " send to " << _sendrequests[i].rank << " OK" << std::endl;
450  }
451  }
452  }
453 
454  // poll receives
455  while (recvs>0)
456  {
457  for (unsigned int i=0; i<_recvrequests.size(); i++)
458  if (!_recvrequests[i].flag)
459  {
460  MPI_Status status;
461  MPI_Test( &(_recvrequests[i].request), &(_recvrequests[i].flag), &status);
462  if (_recvrequests[i].flag)
463  {
464  recvs--;
465  // std::cout << "[" << rank() << "]" << " recv fm " << _recvrequests[i].rank << " OK" << std::endl;
466  }
467 
468  }
469  }
470 
471  // clear request buffers
472  _sendrequests.clear();
473  _recvrequests.clear();
474 #endif
475  }
476 
478  double global_max (double x) const
479  {
480  double res = 0.0;
481  _comm.template allreduce<Dune::Max<double>,double>(&x, &res, 1);
482  return res;
483  }
484 
486  void print (std::ostream& s) const
487  {
488  s << "[" << rank() << "]: Torus " << procs() << " processor(s) arranged as " << dims() << std::endl;
489  for (ProcListIterator i=sendbegin(); i!=sendend(); ++i)
490  {
491  s << "[" << rank() << "]: send to "
492  << "rank=" << i.rank()
493  << " index=" << i.index()
494  << " delta=" << i.delta() << " dist=" << i.distance() << std::endl;
495  }
496  for (ProcListIterator i=recvbegin(); i!=recvend(); ++i)
497  {
498  s << "[" << rank() << "]: recv from "
499  << "rank=" << i.rank()
500  << " index=" << i.index()
501  << " delta=" << i.delta() << " dist=" << i.distance() << std::endl;
502  }
503  }
504 
505  private:
506 
507  void proclists ()
508  {
509  // compile the full neighbor list
510  CommPartner cp;
511  iTupel delta;
512 
513  std::fill(delta.begin(), delta.end(), -1);
514  bool ready = false;
515  iTupel me, nb;
516  me = rank_to_coord(_comm.rank());
517  int index = 0;
518  int last = neighbors()-1;
519  while (!ready)
520  {
521  // find neighbors coordinates
522  for (int i=0; i<d; i++)
523  nb[i] = ( me[i]+_dims[i]+delta[i] ) % _dims[i];
524 
525  // find neighbors rank
526  int nbrank = coord_to_rank(nb);
527 
528  // check if delta is not zero
529  for (int i=0; i<d; i++)
530  if (delta[i]!=0)
531  {
532  cp.rank = nbrank;
533  cp.delta = delta;
534  cp.index = index;
535  _recvlist.push_back(cp);
536  cp.index = last-index;
537  _sendlist.push_front(cp);
538  index++;
539  break;
540  }
541 
542  // next neighbor
543  ready = true;
544  for (int i=0; i<d; i++)
545  if (delta[i]<1)
546  {
547  (delta[i])++;
548  ready=false;
549  break;
550  }
551  else
552  {
553  delta[i] = -1;
554  }
555  }
556 
557  }
558 
559  CollectiveCommunication _comm;
560 
561  iTupel _dims;
562  iTupel _increment;
563  int _tag;
564  std::deque<CommPartner> _sendlist;
565  std::deque<CommPartner> _recvlist;
566 
567  mutable std::vector<CommTask> _sendrequests;
568  mutable std::vector<CommTask> _recvrequests;
569  mutable std::vector<CommTask> _localsendrequests;
570  mutable std::vector<CommTask> _localrecvrequests;
571 
572  };
573 
575  template <class CollectiveCommunication, int d>
576  inline std::ostream& operator<< (std::ostream& s, const Torus<CollectiveCommunication, d> & t)
577  {
578  t.print(s);
579  return s;
580  }
581 }
582 
583 #endif
const iTupel & dims() const
return dimensions of torus
Definition: torus.hh:114
Dune::array< int, d > iTupel
type used to pass tupels in and out
Definition: torus.hh:45
ProcListIterator recvbegin() const
first process in receive list
Definition: torus.hh:351
iTupel delta() const
return distance vector
Definition: torus.hh:293
Torus(CollectiveCommunication comm, int tag, iTupel size, const YLoadBalance< d > *lb)
make partitioner from communicator and coarse mesh size
Definition: torus.hh:73
static double dist(const double *x, const double *y)
Definition: ghmesh.cc:149
virtual void loadbalance(const iTupel &, int, iTupel &) const =0
iTupel coord() const
return own coordinates
Definition: torus.hh:102
double partition(int rank, iTupel origin_in, iTupel size_in, iTupel &origin_out, iTupel &size_out) const
partition the given grid onto the torus and return the piece of the process with given rank; returns ...
Definition: torus.hh:241
int procs() const
return number of processes
Definition: torus.hh:108
int rank_relative(int rank, int dir, int cnt) const
return rank of process where its coordinate in direction dir has offset cnt (handles periodic case) ...
Definition: torus.hh:168
int coord_to_rank(iTupel coord) const
map coordinate in torus to rank using lexicographic ordering
Definition: torus.hh:159
a base class for the yaspgrid partitioning strategy The name might be irritating. It will probably ch...
Definition: partitioning.hh:23
double global_max(double x) const
global max
Definition: torus.hh:478
This file provides tools to partition YaspGrids. If you want to write your own partitioner, inherit from YLoadBalance and implement the loadbalance() method. You can also browse this file for already available useful partitioners, like YaspFixedSizePartitioner.
Definition: torus.hh:42
ProcListIterator recvend() const
last process in receive list
Definition: torus.hh:357
int rank() const
return rank of neighboring process
Definition: torus.hh:287
void print(std::ostream &s) const
print contents of torus object
Definition: torus.hh:486
int rank() const
return own rank
Definition: torus.hh:96
int neighbors() const
return the number of neighbors, which is
Definition: torus.hh:205
int distance() const
return 1-norm of distance vector
Definition: torus.hh:305
ProcListIterator(typename std::deque< CommPartner >::const_iterator iter)
make an iterator
Definition: torus.hh:281
CollectiveCommunication comm() const
return communicator
Definition: torus.hh:126
void exchange() const
exchange messages stored in request buffers; clear request buffers afterwards
Definition: torus.hh:389
bool inside(iTupel c) const
return true if coordinate is inside torus
Definition: torus.hh:138
int color(int rank) const
assign color to given rank
Definition: torus.hh:199
void send(int rank, void *buffer, int size) const
store a send request; buffers are sent in order; handles also local requests with memcpy ...
Definition: torus.hh:363
int color(const iTupel &coord) const
assign color to given coordinate
Definition: torus.hh:176
Definition: torus.hh:278
iTupel rank_to_coord(int rank) const
map rank to coordinate in torus using lexicographic ordering
Definition: torus.hh:146
Include standard header files.
Definition: agrid.hh:59
ProcListIterator sendend() const
end of send list
Definition: torus.hh:345
ProcListIterator sendbegin() const
first process in send list
Definition: torus.hh:339
int index() const
return index in proclist
Definition: torus.hh:299
int dims(int i) const
return dimensions of torus in direction i
Definition: torus.hh:120
bool is_neighbor(iTupel delta, std::bitset< d > periodic) const
return true if neighbor with given delta is a neighbor under the given periodicity ...
Definition: torus.hh:214
void abs(const DofVectorPointer< int > &dofVector)
Definition: dofvector.hh:326
void recv(int rank, void *buffer, int size) const
store a receive request; buffers are received in order; handles also local requests with memcpy ...
Definition: torus.hh:376
Torus()
constructor making uninitialized object
Definition: torus.hh:69
int tag() const
return tag used by torus
Definition: torus.hh:132