OPAL (Object Oriented Parallel Accelerator Library)  2.2.0
OPAL
FieldLayout.hpp
Go to the documentation of this file.
1 // -*- C++ -*-
2 /***************************************************************************
3  *
4  * The IPPL Framework
5  *
6  * This program was prepared by PSI.
7  * All rights in the program are reserved by PSI.
8  * Neither PSI nor the author(s)
9  * makes any warranty, express or implied, or assumes any liability or
10  * responsibility for the use of this software
11  *
12  * Visit www.amas.web.psi for more details
13  *
14  ***************************************************************************/
15 
16 // -*- C++ -*-
17 /***************************************************************************
18  *
19  * The IPPL Framework
20  *
21  *
22  * Visit http://people.web.psi.ch/adelmann/ for more details
23  *
24  ***************************************************************************/
25 
26 // include files
27 #include "FieldLayout/FieldLayout.h"
28 #include "FieldLayout/VRB.h"
29 #include "Message/Communicate.h"
30 #include "Message/Message.h"
31 #include "Utility/DiscMeta.h"
32 #include "Utility/IpplInfo.h"
33 #include "Utility/IpplStats.h"
34 #include "Utility/PAssert.h"
35 
36 
37 #include <cstdlib>
38 #include <limits>
39 
41 // Default constructor, which should only be used if you are going to
42 // call 'initialize' soon after (before using in any context)
43 template<unsigned Dim>
44 FieldLayout<Dim>::FieldLayout() {
45 
46 
47 
48  // we have one more FieldLayout, indicate this
49  //INCIPPLSTAT(incFieldLayouts);
50 
51  // just initialize basic things
52  vnodesPerDirection_m = 0;
53 
54  // for this kind of construction, we just take it that the user is
55  // requesting all parallel axes
56  for (unsigned int dl=0; dl < Dim; ++dl) RequestedLayout[dl] = PARALLEL;
57 }
58 
59 
61 // Constructor which reads in FieldLayout data from a file. If the
62 // file contains data for an equal number of nodes as we are running on,
63 // then that vnode -> pnode mapping will be used. If the file does not
64 // contain info for the same number of pnodes, the vnodes will be
65 // distributed in some other manner.
66 template<unsigned Dim>
67 FieldLayout<Dim>::FieldLayout(const char *filename) {
68 
69 
70 
71  // we have one more FieldLayout, indicate this
72  //INCIPPLSTAT(incFieldLayouts);
73 
74  // try to initialize ourselves, by reading the info from the file.
75  vnodesPerDirection_m = 0;
76 
77  // for this kind of construction, we just take it that the user is
78  // requesting all parallel axes
79  for (unsigned int dl=0; dl < Dim; ++dl) RequestedLayout[dl] = PARALLEL;
80 
81  // read in data from file
82  read(filename);
83 }
84 
85 
87 // Destructor: Everything deletes itself automatically ... the base
88 // class destructors inform all the FieldLayoutUser's we're going away.
89 template<unsigned Dim>
90 FieldLayout<Dim>::~FieldLayout() {
91  if (vnodesPerDirection_m != 0)
92  delete [] vnodesPerDirection_m;
93 }
94 
95 
97 
98 // Initialization functions, only to be called by the user of FieldLayout
99 // objects when the FieldLayout was created using the default constructor;
100 // otherwise these are only called internally by the various non-default
101 // FieldLayout constructors:
102 
103 //-----------------------------------------------------------------------------
104 // These specify only a total number of vnodes, allowing the constructor
105 // complete control on how to do the vnode partitioning of the index space:
106 
107 template<unsigned Dim>
108 void
109 FieldLayout<Dim>::initialize(const NDIndex<Dim>& domain,
110  e_dim_tag *p, int vnodes) {
111 
112 
113  setup(domain, p, vnodes);
114 }
115 
116 
117 template<unsigned Dim>
118 void
119 FieldLayout<Dim>::initialize(const Index& i1, e_dim_tag p1, int vnodes) {
120 
121 
122 
123  PInsist(Dim==1,
124  "Number of arguments does not match dimension of FieldLayout!!");
125  NDIndex<Dim> ndi(i1);
126  setup(ndi,&p1,vnodes);
127 }
128 
129 template<unsigned Dim>
130 void
131 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2,
132  e_dim_tag p1, e_dim_tag p2, int vnodes) {
133 
134 
135 
136  PInsist(Dim==2,
137  "Number of arguments does not match dimension of FieldLayout!!");
138  e_dim_tag par[Dim];
139  par[0] = p1;
140  par[1] = p2;
141  NDIndex<Dim> ndi;
142  ndi[0] = i1;
143  ndi[1] = i2;
144  setup(ndi,par,vnodes);
145 }
146 template<unsigned Dim>
147 void
148 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2, const Index& i3,
149  e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
150  int vnodes) {
151 
152  PInsist(Dim==3,
153  "Number of arguments does not match dimension of FieldLayout!!");
154  e_dim_tag par[Dim];
155  par[0] = p1;
156  par[1] = p2;
157  par[2] = p3;
158  NDIndex<Dim> ndi;
159  ndi[0] = i1;
160  ndi[1] = i2;
161  ndi[2] = i3;
162  setup(ndi,par,vnodes);
163 }
164 template<unsigned Dim>
165 void
166 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2, const Index& i3,
167  const Index& i4,
168  e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
169  e_dim_tag p4,
170  int vnodes) {
171 
172 
173  PInsist(Dim==4,
174  "Number of arguments does not match dimension of FieldLayout!!");
175  e_dim_tag par[Dim];
176  par[0] = p1;
177  par[1] = p2;
178  par[2] = p3;
179  par[3] = p4;
180  NDIndex<Dim> ndi;
181  ndi[0] = i1;
182  ndi[1] = i2;
183  ndi[2] = i3;
184  ndi[3] = i4;
185  setup(ndi,par,vnodes);
186 }
187 template<unsigned Dim>
188 void
189 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2, const Index& i3,
190  const Index& i4, const Index& i5,
191  e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
192  e_dim_tag p4, e_dim_tag p5,
193  int vnodes) {
194 
195  PInsist(Dim==5,
196  "Number of arguments does not match dimension of FieldLayout!!");
197  e_dim_tag par[Dim];
198  par[0] = p1;
199  par[1] = p2;
200  par[2] = p3;
201  par[3] = p4;
202  par[4] = p5;
203  NDIndex<Dim> ndi;
204  ndi[0] = i1;
205  ndi[1] = i2;
206  ndi[2] = i3;
207  ndi[3] = i4;
208  ndi[4] = i5;
209  setup(ndi,par,vnodes);
210 }
211 template<unsigned Dim>
212 void
213 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2, const Index& i3,
214  const Index& i4, const Index& i5, const Index& i6,
215  e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
216  e_dim_tag p4, e_dim_tag p5, e_dim_tag p6,
217  int vnodes) {
218 
219  PInsist(Dim==6,
220  "Number of arguments does not match dimension of FieldLayout!!");
221  e_dim_tag par[Dim];
222  par[0] = p1;
223  par[1] = p2;
224  par[2] = p3;
225  par[3] = p4;
226  par[4] = p5;
227  par[5] = p6;
228  NDIndex<Dim> ndi;
229  ndi[0] = i1;
230  ndi[1] = i2;
231  ndi[2] = i3;
232  ndi[3] = i4;
233  ndi[4] = i5;
234  ndi[5] = i6;
235  setup(ndi,par,vnodes);
236 }
237 //-----------------------------------------------------------------------------
238 
239 //-----------------------------------------------------------------------------
240 // These specify both the total number of vnodes and the numbers of vnodes
241 // along each dimension for the partitioning of the index space. Obviously
242 // this restricts the number of vnodes to be a product of the numbers along
243 // each dimension (the constructor implementation checks this):
244 
245 template<unsigned Dim>
246 void
247 FieldLayout<Dim>::initialize(const NDIndex<Dim>& domain,
248  e_dim_tag *p, unsigned* vnodesPerDirection,
249  bool recurse, int vnodes) {
250 
251 
252  // Default to correct total vnodes:
253  unsigned vnodesProduct = 1;
254  for (unsigned int d=0; d<Dim; d++) vnodesProduct *= vnodesPerDirection[d];
255  if (vnodes == -1) vnodes = vnodesProduct;
256  // Verify than total vnodes is product of per-dimension vnode counts:
257  if ((unsigned int) vnodes != vnodesProduct) {
258  ERRORMSG("FieldLayout constructor: "
259  << "(vnodes != vnodesPerDirection[0]*vnodesPerDirection[1]*"
260  << "...*vnodesPerDirection[" << Dim-1 << "])"
261  << " ; vnodesPerDirection[0]*vnodesPerDirection[1]*"
262  << "...*vnodesPerDirection[" << Dim-1 << "] = "
263  << vnodesProduct << " ; vnodes = " << vnodes << endl);
264  }
265  setup(domain, p, vnodesPerDirection,recurse,vnodes);
266 }
267 
268 template<unsigned Dim>
269 void
270 FieldLayout<Dim>::initialize(const Index& i1, e_dim_tag p1,
271  unsigned vnodes1, bool recurse, int vnodes) {
272 
273 
274 
275  PInsist(Dim==1,
276  "Number of arguments does not match dimension of FieldLayout!!");
277  NDIndex<Dim> ndi(i1);
278  setup(ndi,&p1,&vnodes1,recurse,vnodes);
279 }
280 
281 template<unsigned Dim>
282 void
283 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2,
284  e_dim_tag p1, e_dim_tag p2,
285  unsigned vnodes1, unsigned vnodes2,
286  bool recurse, int vnodes) {
287 
288  PInsist(Dim==2,
289  "Number of arguments does not match dimension of FieldLayout!!");
290  e_dim_tag par[Dim];
291  par[0] = p1;
292  par[1] = p2;
293  NDIndex<Dim> ndi;
294  ndi[0] = i1;
295  ndi[1] = i2;
296  unsigned vnodesPerDirection[Dim];
297  vnodesPerDirection[0] = vnodes1;
298  vnodesPerDirection[1] = vnodes2;
299  setup(ndi,par,vnodesPerDirection,recurse,vnodes);
300 }
301 template<unsigned Dim>
302 void
303 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2, const Index& i3,
304  e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
305  unsigned vnodes1, unsigned vnodes2,
306  unsigned vnodes3,
307  bool recurse, int vnodes) {
308 
309  PInsist(Dim==3,
310  "Number of arguments does not match dimension of FieldLayout!!");
311  e_dim_tag par[Dim];
312  par[0] = p1;
313  par[1] = p2;
314  par[2] = p3;
315  NDIndex<Dim> ndi;
316  ndi[0] = i1;
317  ndi[1] = i2;
318  ndi[2] = i3;
319  unsigned vnodesPerDirection[Dim];
320  vnodesPerDirection[0] = vnodes1;
321  vnodesPerDirection[1] = vnodes2;
322  vnodesPerDirection[2] = vnodes3;
323  setup(ndi,par,vnodesPerDirection,recurse,vnodes);
324 }
325 template<unsigned Dim>
326 void
327 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2, const Index& i3,
328  const Index& i4,
329  e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
330  e_dim_tag p4,
331  unsigned vnodes1, unsigned vnodes2,
332  unsigned vnodes3, unsigned vnodes4,
333  bool recurse, int vnodes) {
334 
335  PInsist(Dim==4,
336  "Number of arguments does not match dimension of FieldLayout!!");
337  e_dim_tag par[Dim];
338  par[0] = p1;
339  par[1] = p2;
340  par[2] = p3;
341  par[3] = p4;
342  NDIndex<Dim> ndi;
343  ndi[0] = i1;
344  ndi[1] = i2;
345  ndi[2] = i3;
346  ndi[3] = i4;
347  unsigned vnodesPerDirection[Dim];
348  vnodesPerDirection[0] = vnodes1;
349  vnodesPerDirection[1] = vnodes2;
350  vnodesPerDirection[2] = vnodes3;
351  vnodesPerDirection[3] = vnodes4;
352  setup(ndi,par,vnodesPerDirection,recurse,vnodes);
353 }
354 template<unsigned Dim>
355 void
356 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2, const Index& i3,
357  const Index& i4, const Index& i5,
358  e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
359  e_dim_tag p4, e_dim_tag p5,
360  unsigned vnodes1, unsigned vnodes2,
361  unsigned vnodes3, unsigned vnodes4,
362  unsigned vnodes5,
363  bool recurse, int vnodes) {
364 
365  PInsist(Dim==5,
366  "Number of arguments does not match dimension of FieldLayout!!");
367  e_dim_tag par[Dim];
368  par[0] = p1;
369  par[1] = p2;
370  par[2] = p3;
371  par[3] = p4;
372  par[4] = p5;
373  NDIndex<Dim> ndi;
374  ndi[0] = i1;
375  ndi[1] = i2;
376  ndi[2] = i3;
377  ndi[3] = i4;
378  ndi[4] = i5;
379  unsigned vnodesPerDirection[Dim];
380  vnodesPerDirection[0] = vnodes1;
381  vnodesPerDirection[1] = vnodes2;
382  vnodesPerDirection[2] = vnodes3;
383  vnodesPerDirection[3] = vnodes4;
384  vnodesPerDirection[4] = vnodes5;
385  setup(ndi,par,vnodesPerDirection,recurse,vnodes);
386 }
387 template<unsigned Dim>
388 void
389 FieldLayout<Dim>::initialize(const Index& i1, const Index& i2, const Index& i3,
390  const Index& i4, const Index& i5, const Index& i6,
391  e_dim_tag p1, e_dim_tag p2, e_dim_tag p3,
392  e_dim_tag p4, e_dim_tag p5, e_dim_tag p6,
393  unsigned vnodes1, unsigned vnodes2,
394  unsigned vnodes3, unsigned vnodes4,
395  unsigned vnodes5, unsigned vnodes6,
396  bool recurse, int vnodes) {
397 
398  PInsist(Dim==6,
399  "Number of arguments does not match dimension of FieldLayout!!");
400  e_dim_tag par[Dim];
401  par[0] = p1;
402  par[1] = p2;
403  par[2] = p3;
404  par[3] = p4;
405  par[4] = p5;
406  par[5] = p6;
407  NDIndex<Dim> ndi;
408  ndi[0] = i1;
409  ndi[1] = i2;
410  ndi[2] = i3;
411  ndi[3] = i4;
412  ndi[4] = i5;
413  ndi[5] = i6;
414  unsigned vnodesPerDirection[Dim];
415  vnodesPerDirection[0] = vnodes1;
416  vnodesPerDirection[1] = vnodes2;
417  vnodesPerDirection[2] = vnodes3;
418  vnodesPerDirection[3] = vnodes4;
419  vnodesPerDirection[4] = vnodes5;
420  vnodesPerDirection[5] = vnodes6;
421  setup(ndi,par,vnodesPerDirection,recurse,vnodes);
422 }
423 //-----------------------------------------------------------------------------
424 
425 
426 //-----------------------------------------------------------------------------
427 // A version of initialize that takes the total domain, and iterators
428 // over the subdomains the user wants along with iterators over the
429 // node assignments. No communication is done
430 // so these lists must match on all nodes. A bit of error checking
431 // is done for overlapping blocks and illegal nodes, but not exhaustive
432 // error checking.
433 
434 template<unsigned Dim>
435 void
436 FieldLayout<Dim>::initialize(const NDIndex<Dim> &domain,
437  const NDIndex<Dim> *dombegin,
438  const NDIndex<Dim> *domend,
439  const int *nbegin, const int *nend)
440 {
441 
442  // Loop variables
443  int i, j;
444  unsigned int d;
445 
446  // Save the total domain.
447  Domain = domain;
448 
449  // Find the number of vnodes requested
450  int vnodes = (nend - nbegin);
451  PInsist(vnodes > 0,
452  "A user-specified FieldLayout must have at least one vnode.");
453  PInsist(vnodes == (domend - dombegin),
454  "A user-specified FieldLayout must have equal length node and domain lists");
455 
456  // Since we don't know any differently, indicate the requested
457  // layout is all parallel
458  for (d = 0; d < Dim; ++d)
459  RequestedLayout[d] = PARALLEL;
460 
461  // This is not a grid-like layout, so set the pointer for this to 0
462  vnodesPerDirection_m = 0;
463 
464  // Create the empty remote vnode list
465  ac_domain_vnodes *remote_ac = new ac_domain_vnodes(domain);
466  typedef typename ac_gc_domain_vnodes::value_type vntype;
467  Remotes_ac.insert( vntype(gc0(), remote_ac) );
468 
469  // Loop through the vnodes, and add them to our local or remote lists.
470  // Do a sanity check on the vnodes, making sure each one does not
471  // intersect any other. Also, add up the size of each vnode, if it
472  // does not equal the size of the total domain, there are holes
473  // and it is an error.
474  size_t coverage = 0;
475  for (i = 0; i < vnodes; ++i) {
476  // Compare to other vnodes
477  for (j = (i+1); j < vnodes; ++j) {
478  PInsist(! (dombegin[i].touches(dombegin[j])),
479  "A user-specified FieldLayout cannot have overlapping domains.");
480  }
481 
482  // Make sure the processor ID is OK
483  PInsist(nbegin[i] >= 0 && nbegin[i] < Ippl::getNodes(),
484  "A user-specified FieldLayout must have legal node assignments.");
485 
486  // Add in the volume of this domain
487  coverage += dombegin[i].size();
488 
489  // Create a Vnode for this domain
490  Vnode<Dim> *vnode = new Vnode<Dim>(dombegin[i], nbegin[i], i);
491  typedef typename ac_id_vnodes::value_type v1;
492  typedef typename ac_domain_vnodes::value_type v2;
493  bool nosplit = (dombegin[i].size() < 2);
494 
495  // Based on the assigned node, add to our local or remote lists
496  if (nbegin[i] == Ippl::myNode())
497  Local_ac.insert(v1(Unique::get(), vnode));
498  else
499  remote_ac->insert(v2(dombegin[i], vnode), nosplit);
500  }
501 
502  // Check the coverage, make sure it is complete
503  PInsist(coverage == domain.size(),
504  "A user-specified FieldLayout must completely cover the domain.");
505 
506  // At the end, calculate the widthds of all the vnodes
507  calcWidths();
508 }
509 
510 
512 
513 //
514 // Now, the meat of the construction.
515 //
516 
517 //-----------------------------------------------------------------------------
518 //This setup() specifies only a total number of vnodes, taking complete control
519 //on how to do the vnode partitioning of the index space:
520 
521 template<unsigned Dim>
522 void
523 FieldLayout<Dim>::setup(const NDIndex<Dim>& domain,
524  e_dim_tag *userflags, int vnodes)
525 {
526 
527 
528  // we have one more FieldLayout, indicate this
529  //INCIPPLSTAT(incFieldLayouts);
530 
531  // Find the number processors.
532  int nprocs = Ippl::getNodes();
533  int myproc = Ippl::myNode();
534 
535  // If the user didn't specify the number of vnodes, make it equal nprocs
536  if (vnodes <= 0) vnodes = nprocs;
537 
538  Inform dbgmsg("FieldLayout::setup", INFORM_ALL_NODES);
539  // dbgmsg << "*** Domain=" << domain << ", nprocs=" << nprocs;
540  // dbgmsg << ", myproc=" << myproc << ", vnodes=" << vnodes << endl;
541 
542  // If the user did not specify parallel/serial flags then make all parallel.
543  int parallel_count = 0;
544  unsigned int flagdim = 0;
545  long totparelems = 1;
546  for (flagdim=0; flagdim < Dim; ++flagdim) {
547  if (userflags == 0)
548  RequestedLayout[flagdim] = PARALLEL;
549  else
550  RequestedLayout[flagdim] = userflags[flagdim];
551  if (RequestedLayout[flagdim] == PARALLEL) {
552  parallel_count++;
553  totparelems *= domain[flagdim].length();
554  }
555  }
556  // Make sure at least one of the parallel/serial flags is parallel
557  PInsist(parallel_count>0,"At least one dimension of a FieldLayout must be PARALLEL!");
558 
559  // Check to see if we have too few elements to partition. If so, reduced
560  // the number of vnodes (if necessary) to just the number of elements along
561  // parallel dims.
562  if (totparelems < vnodes) {
563  //dbgmsg << "Total parallel lengths = " << totparelems << "; reducing ";
564  //dbgmsg << "vnodes from " << vnodes << " to " << totparelems << endl;
565  vnodes = totparelems;
566  }
567 
568  e_dim_tag *flags = RequestedLayout;
569 
570  // Recursively split the domain until we have generated all the domains.
571  Domain = domain;
572  NDIndex<Dim> *domains_c = new NDIndex<Dim>[vnodes];
573  NDIndex<Dim> *copy_c = new NDIndex<Dim>[vnodes];
574  NDIndex<Dim> leftDomain ;
575 
576  // Start with the whole domain.
577  domains_c[0] = domain;
578  int v;
579  unsigned int d=0;
580 
581  int v1,v2,rm,vtot,vl,vr;
582  double a,lmax,len;
583  vnodesPerDirection_m = 0;
584  for (v=vnodes,rm=0;v>1;v/=2) { rm += (v % 2); }
585  if (rm == 0) {
586 
587  // vnodes is a power of 2
588 
589  // For power-of-two vnodes, allocate storage for vnodesPerDirection_m:
590  //don't--tjw vnodesPerDirection_m = new unsigned[Dim];
591  // Fill in 1 for starters; remains 1 for serial directions:
592  //don't--tjw for (int d2=0; d2<Dim; d2++) vnodesPerDirection_m[d2] = 1;
593 
594  for (v=1; v<vnodes; v*=2) {
595  // Go to the next parallel dimension.
596  while(flags[d] != PARALLEL) if(++d == Dim) d = 0;
597 
598  // Split all the current vnodes.
599  int i,j;
600  for (i=0, j=0; i<v; ++i, j+=2)
601  // Split to the left and to the right, saving both.
602  domains_c[i].split( copy_c[j] , copy_c[j+1] , d );
603  // Copy back.
604  std::copy(copy_c,copy_c+v*2, domains_c);
605 
606  //don't--tjw vnodesPerDirection_m[d] *= 2; // Update to reflect split
607 
608  // On to the next dimension.
609  if (++d == Dim) d = 0;
610  }
611 
612  } else {
613 
614  vtot = 1; // count the number of vnodes to make sure that it worked
615  // vnodes is not a power of 2 so we need to do some fancy splitting
616  // sorry... this would be much cleaner with recursion
617  /*
618  The way this works is to recursively split on the longest dimension.
619  Suppose you request 11 vnodes. It will split the longest dimension
620  in the ratio 5:6 and put the new domains in node 0 and node 5. Then
621  it splits the longest dimension of the 0 domain and puts the results
622  in node 0 and node 2 and then splits the longest dimension of node 5
623  and puts the results in node 5 and node 8. etc.
624  The logic is kind of bizarre, but it works.
625  */
626  for (v=1; v<2*vnodes; ++v) {
627  // kind of reverse the bits of v
628  for (v2=v,v1=1;v2>1;v2/=2) { v1 = 2*v1+(v2%2); }
629  vl = 0; vr = vnodes;
630  while (v1>1) {
631  if ((v1%2)==1) {
632  vl=vl+(vr-vl)/2;
633  } else {
634  vr=vl+(vr-vl)/2;
635  }
636  v1/=2;
637  }
638  v2=vl+(vr-vl)/2;
639 
640  if (v2>vl) {
641  a = v2-vl;
642  a /= vr-vl;
643  vr=v2;
644  leftDomain=domains_c[vl];
645  lmax=0;
646  d=std::numeric_limits<unsigned int>::max();
647  for (unsigned int dd=0;dd<Dim;++dd) {
648  if ( flags[dd] == PARALLEL ) {
649  if ((len = leftDomain[dd].length()) > lmax) {
650  lmax = len;
651  d = dd;
652  }
653  }
654  }
655  domains_c[vl].split( domains_c[vl] , domains_c[vr] , d , a);
656  ++vtot;
657  }
658  }
659  v=vtot;
660  }
661  // Make sure we had a power of two number of vnodes.
662  PAssert_EQ( v, vnodes );
663 
664  // Now make the vnodes, using the domains just generated.
665  // Some of them we store in the local list, others in the remote.
666  ac_domain_vnodes *remote_ac = new ac_domain_vnodes( domain );
667  typedef typename ac_gc_domain_vnodes::value_type vtype;
668  Remotes_ac.insert( vtype(gc0(),remote_ac) );
669  for (v=0; v<vnodes; ++v)
670  {
671  int p = (v*nprocs)/vnodes;
672  bool nosplit = (domains_c[v].size() < 2);
673 
674  // Add v arg to Vnode constructor 3/19/98 --tjw:
675  Vnode<Dim> *vnode = new Vnode<Dim>(domains_c[v], p, v);
676  typedef typename ac_id_vnodes::value_type v1;
677  typedef typename ac_domain_vnodes::value_type v2;
678  if ( p==myproc )
679  Local_ac.insert(v1(Unique::get(),vnode));
680  else
681  remote_ac->insert(v2(domains_c[v], vnode), nosplit);
682  }
683 
684  // Delete the memory we allocated.
685  delete [] domains_c;
686  delete [] copy_c;
687 
688  // Calculate the widths of all the vnodes
689  calcWidths();
690 }
691 
692 //-----------------------------------------------------------------------------
693 
694 //-----------------------------------------------------------------------------
695 //
696 // This setup() specifies both the total number of vnodes and the numbers of
697 // vnodes along each dimension for the partitioning of the index space, where
698 // the user wants this extra control over the partitioning of the index
699 // space. Obviously this restricts the number of vnodes to be a product of the
700 // numbers along each dimension (the constructor implementation checks this).
701 //
702 // The last argument is a bool for the algorithm to use for assigning
703 // vnodes to processors.
704 // If it is false, hand the vnodes to the processors in a very simple
705 // but probably inefficient manner.
706 // If it is true, use a binary recursive algorithm. This will usually be
707 // more efficient because it will generate less communication, but
708 // it will sometimes fail, particularly near the case of one vnode per
709 // processor.
710 //
711 //-----------------------------------------------------------------------------
712 
713 template<unsigned Dim>
714 void
715 FieldLayout<Dim>::setup(const NDIndex<Dim>& domain,
716  e_dim_tag *userflags, unsigned* vnodesPerDirection,
717  bool recurse, int vnodes)
718 {
719 
720  // Find the number processors.
721  int nprocs = Ippl::getNodes();
722  int myproc = Ippl::myNode();
723 
724  // The number of vnodes must have been specified or computed by now:
725  if (vnodes <= 0) ERRORMSG("FieldLayout::setup(): vnodes <= 0 "
726  << "for a vnodes-per-direction product "
727  << "specification; not allowed." << endl);
728 
729  // Loop indices:
730  int v;
731  unsigned int d, d2, vl;
732 
733  int parallel_count = 0;
734  unsigned int flagdim = 0;
735  for (flagdim=0; flagdim < Dim; ++flagdim) {
736  if (userflags == 0)
737  RequestedLayout[flagdim] = PARALLEL;
738  else
739  RequestedLayout[flagdim] = userflags[flagdim];
740 
741  // keep track of the number of parallel dimensions; we need at least one
742  parallel_count += (RequestedLayout[flagdim] == PARALLEL);
743 
744  // make sure any SERIAL dimensions request only one vnode along them:
745  if (RequestedLayout[flagdim] == SERIAL) {
746  bool chk = vnodesPerDirection[flagdim] == 1;
747  PInsist(chk,"SERIAL layout specified, yet vnodesPerDirection is not 1!");
748  }
749  }
750 
751  // Make sure at least one of the parallel/serial flags is parallel
752  PInsist(parallel_count>0,"At least one dimension of a FieldLayout must be PARALLEL!");
753 
754  // Allocate and store vnodesPerDirection_m data-member array:
755  vnodesPerDirection_m = new unsigned[Dim];
756  for (d=0; d<Dim; d++) vnodesPerDirection_m[d] = vnodesPerDirection[d];
757 
758  // The domain of this FieldLayout object under construction:
759  Domain = domain;
760 
761  // Set up a container of NDIndex's to store the index ranges for all vnodes:
762  NDIndex<Dim> *domains_c = new NDIndex<Dim>[vnodes];
763 
764  // Divide the numbers of elements by the numbers of vnodes (along each
765  // dimension). All vnodes will have at least this many elements. Some will
766  // have one extra element:
767  unsigned elementsPerVnode[Dim];
768  unsigned nLargerVnodes[Dim];
769  for (d=0; d<Dim; d++) {
770  elementsPerVnode[d] = domain[d].length()/vnodesPerDirection[d];
771  nLargerVnodes[d] = domain[d].length() % vnodesPerDirection[d];
772  }
773 
774  // Set up the base, bound, and stride for the index range of all
775  // vnodes. Organize by "vnode level." For this kind of partitioning we have a
776  // the equivalent of a Dim-dimensional array of vnodes, with the number of
777  // elements in dimension d being vnodesPerDirection[d]. By "vnode level" we
778  // mean the index along a dimension in the "vnode array" of a given vnode. We
779  // can store the vnode-index bases and bounds 2D arrays of arrays; we only
780  // need a 1D array to store the vnode strides because they are the same
781  // everywhere (and the same as the global strides for each dimension):
782  int stride[Dim];
783  for (d=0; d<Dim; d++) stride[d] = domain[d].stride();
784  int* base[Dim];
785  int* bound[Dim];
786  for (d=0; d<Dim; d++) {
787  base[d] = new int[vnodesPerDirection[d]];
788  bound[d] = new int[vnodesPerDirection[d]];
789  }
790  int length;
791  for (d=0; d<Dim; d++) {
792  // Start things off for the zeroth vnode level using the global index base:
793  length = elementsPerVnode[d];
794  if (nLargerVnodes[d] > 0) length += 1;
795  base[d][0] = domain[d].first();
796  bound[d][0] = base[d][0] + (length-1)*stride[d];
797  // Now go through all the other vnode levels
798  for (vl=1; vl < vnodesPerDirection[d]; vl++) {
799  base[d][vl] = bound[d][vl-1] + stride[d];
800  length = elementsPerVnode[d];
801  if (vl < nLargerVnodes[d]) length += 1;
802  bound[d][vl] = base[d][vl] +(length-1)*stride[d];
803  }
804  }
805 
806  // Now actually initialize the values in the container of NDIndex's which
807  // represent each vnode's subdomain:
808  for (v=0; v<vnodes; v++) {
809  for (d=0; d<Dim; d++) {
810  // Compute this vnode's level in this direction:
811  unsigned denom = 1;
812  for (d2=0; d2<d; d2++) denom *= vnodesPerDirection[d2];
813  int vnodeLevel = (v/denom) % vnodesPerDirection[d];
814  // Now use the precomputed base, bound values to set the Index values:
815  domains_c[v][d] = Index(base[d][vnodeLevel],bound[d][vnodeLevel],
816  stride[d]);
817  denom = 1; // for debugging
818  }
819  }
820 
821  // Now find what processor each vnode will end up on.
822  // This is done with a recursive bisection algorithm.
823  // This produces fairly squarish blocks -- and therefore
824  // less communication at the expense of each processor getting
825  // a less balanced load.
826  int *vnodeProcs = new int[vnodes];
827  int *sizes = new int[Dim];
828  for (v=0; (unsigned int) v<Dim; ++v)
829  sizes[v] = vnodesPerDirection[v];
830 
831  // If we have been instructed to use recursive bisection, do that.
832  // if not, deal them out in a simple manner.
833  if ( recurse )
834  vnodeRecursiveBisection(Dim,sizes,nprocs,vnodeProcs);
835  else
836  for ( v=0; v<vnodes; ++v )
837  vnodeProcs[v] = (v*nprocs)/vnodes;
838 
839 
840  // Now make the vnodes, using the domains just generated.
841  // Some of them we store in the local list, others in the remote.
842  ac_domain_vnodes *remote_ac = new ac_domain_vnodes( domain );
843  typedef typename ac_gc_domain_vnodes::value_type v1;
844  Remotes_ac.insert(v1(gc0(),remote_ac) );
845  for (v=0; v<vnodes; ++v) {
846  int p = vnodeProcs[v];
847  // Add v arg to Vnode constructor 3/19/98 --tjw:
848  Vnode<Dim> *vnode = new Vnode<Dim>(domains_c[v], p, v);
849  typedef typename ac_id_vnodes::value_type v2;
850  if ( p==myproc )
851  Local_ac.insert(v2(Unique::get(),vnode) );
852  else
853  // For domains of extent 1 in all directions, must call
854  // DomainMap::insert() with the noSplit flag set to true, to avoid an
855  // assertion failure in Index::split. Also do this for domains of
856  // extent ZERO. These are not supposed to happen in IPPL, but when the
857  // do when something like BinaryRepartition creates them, accomodating
858  // them here will prevent subsequent code hangs or other errors and allow
859  // recovery such as rejection of the result of BinaryRepartition (discard
860  // the FieldLayout with zero-size domains and go on with the original
861  // layout unchanged). (tjw)
862  // if (domains_c[v].size() == 1) {
863  if (domains_c[v].size() <= 1) {
864  typedef typename ac_domain_vnodes::value_type v1;
865  remote_ac->insert(v1(domains_c[v], vnode), true);
866  } else {
867  typedef typename ac_domain_vnodes::value_type v1;
868  remote_ac->insert(v1(domains_c[v], vnode), false);
869  }
870  }
871 
872  // Delete the memory we allocated.
873  delete [] domains_c;
874  delete [] vnodeProcs;
875  delete [] sizes;
876  for (d=0; d<Dim; d++) {
877  delete [] base[d];
878  delete [] bound[d];
879  }
880 
881  // Calculate the widths of all the vnodes
882  calcWidths();
883 }
884 //-----------------------------------------------------------------------------
885 
886 
887 //-----------------------------------------------------------------------------
888 // Return number of vnodes along a direction.
889 template<unsigned Dim>
890 unsigned FieldLayout<Dim>::
891 getVnodesPerDirection(unsigned dir) {
892 
893  // If not set, then not appropriate to be calling this function. Example:
894  // now-power-of-two-constructed FieldLayout which did *not* specify numbers
895  // of vnodes along each direction:
896  PAssert(vnodesPerDirection_m);
897 
898  return(vnodesPerDirection_m[dir]);
899 }
900 //-----------------------------------------------------------------------------
901 
902 
903 //---------------------------------------------------------------------
904 
905 template<unsigned Dim>
906 void FieldLayout<Dim>::new_gc_layout(const GuardCellSizes<Dim>& gc)
907 {
908 
909 
910 
911  // Build the guarded domain.
912  NDIndex<Dim> guarded_domain( AddGuardCells(Domain,gc) );
913  // Build a container for vnodes in that domain.
914  ac_domain_vnodes *gr = new ac_domain_vnodes(guarded_domain);
915  // Record pointer to that container using gc as the key.
916  typedef typename ac_gc_domain_vnodes::value_type v1;
917  Remotes_ac.insert(v1(gc,gr) );
918  // Get the container of vnodes stored w/o guard cells.
919  ac_domain_vnodes &v0 = *Remotes_ac[ gc0() ];
920  // Loop over all the remote vnodes.
921  for (iterator_dv v_i = v0.begin(); v_i != v0.end(); ++v_i)
922  {
923  // Build the domain for this vnode with gc guard cells.
924  NDIndex<Dim> domain(AddGuardCells((*v_i).first,gc));
925  // Record pointer to this vnode with these guard cells.
926  typedef typename ac_domain_vnodes::value_type v2;
927  gr->insert(v2(domain,(*v_i).second) );
928  }
929 }
930 
931 //----------------------------------------------------------------------
932 template<unsigned Dim>
933 FieldLayout<Dim>::FieldLayout(const NDIndex<Dim>& domain,
934  const NDIndex<Dim>* idx_begin,
935  const NDIndex<Dim>* idx_end)
936  : Domain(domain)
937 {
938 
939  // we have one more FieldLayout, indicate this
940  //INCIPPLSTAT(incFieldLayouts);
941 
942  // for this kind of construction, we just take it that the user is
943  // requesting all parallel axes
944  for (unsigned int dl=0; dl < Dim; ++dl) RequestedLayout[dl] = PARALLEL;
945 
946  // Build Vnodes for each of the local domains.
947  vnodesPerDirection_m = 0;
948  int mynode = Ippl::Comm->myNode();
949  for (const NDIndex<Dim>* p=idx_begin; p!=idx_end; ++p)
950  {
951  typedef typename ac_id_vnodes::value_type v1;
952  Local_ac.insert(v1(Unique::get(), new Vnode<Dim>(*p,mynode)));
953  }
954 
955  // Everybody broadcasts their new local domains to everybody.
956  // Build a message with the local domains and the id's for them.
957  Message* bcast_mess = new Message;
958  int count = idx_end - idx_begin;
959  bcast_mess->put(count);
960  // ::putMessage(*bcast_mess,idx_begin,idx_end);
961  const NDIndex<Dim>* indxmsg = idx_begin;
962  while (indxmsg != idx_end)
963  (indxmsg++)->putMessage(*bcast_mess);
964  // Send it to everybody except yourself, and record the number sent.
965  int tag = Ippl::Comm->next_tag( F_REPARTITION_BCAST_TAG , F_TAG_CYCLE );
966  int node_count = Ippl::Comm->broadcast_others(bcast_mess,tag);
967  // Create the container for the remote vnodes.
968  ac_domain_vnodes *remote_ac = new ac_domain_vnodes( Domain );
969  typedef typename ac_gc_domain_vnodes::value_type v1;
970  Remotes_ac.insert(v1(gc0(),remote_ac) );
971 
972  // Loop until we receive a message from each other node.
973  while ((--node_count)>=0)
974  {
975  // Receive a broadcast message from any node.
976  int other_node = COMM_ANY_NODE;
977  Message *recv_mess = Ippl::Comm->receive_block(other_node,tag);
978  PAssert(recv_mess);
979  // Extract the number of vnodes coming in.
980  int count = 0;
981  recv_mess->get(count);
982  // Now get the domains for the vnodes from the message.
983  NDIndex<Dim> p;
984  while (count-- > 0) {
985  p.getMessage(*recv_mess);
986  Vnode<Dim> *vnode = new Vnode<Dim>(p, other_node);
987  // For domains of extent 1 in all directions, must call
988  // DomainMap::insert() with the noSplit flag set to true, to avoid an
989  // assertion failure in Index::split. Also do this for domains of
990  // extent ZERO. These are not supposed to happen in IPPL, but when the
991  // do when something like BinaryRepartition creates them, accomodating
992  // them here will prevent subsequent code hangs or other errors and
993  // allow recovery such as rejection of the result of BinaryRepartition
994  // (discard the FieldLayout with zero-size domains and go on with the
995  // original layout unchanged). (tjw)
996  // if (p.size() == 1) {
997  if (p.size() <= 1) {
998  typedef typename ac_domain_vnodes::value_type v1;
999  remote_ac->insert(v1(p, vnode), true);
1000  } else {
1001  typedef typename ac_domain_vnodes::value_type v1;
1002  remote_ac->insert(v1(p, vnode), false);
1003  }
1004  }
1005  delete recv_mess;
1006  }
1007 
1008  // Calculate the widths of all the vnodes
1009  calcWidths();
1010 }
1011 
1012 
1013 //----------------------------------------------------------------------
1014 // This differs from the previous ctor in that it allows preservation of
1015 // global Vnode integer ID numbers associated with the input Vnodes. --tjw
1016 template<unsigned Dim>
1017 FieldLayout<Dim>::FieldLayout(const NDIndex<Dim>& domain,
1018  const Vnode<Dim>* idx_begin,
1019  const Vnode<Dim>* idx_end)
1020  : Domain(domain)
1021 {
1022 
1023  // we have one more FieldLayout, indicate this
1024  //INCIPPLSTAT(incFieldLayouts);
1025 
1026  // for this kind of construction, we just take it that the user is
1027  // requesting all parallel axes
1028  for (unsigned int dl=0; dl < Dim; ++dl) RequestedLayout[dl] = PARALLEL;
1029 
1030  // Build Vnodes for each of the local domains.
1031  vnodesPerDirection_m = 0;
1032  int mynode = Ippl::Comm->myNode();
1033  for (const Vnode<Dim>* p=idx_begin; p!=idx_end; ++p) {
1034  typedef typename ac_id_vnodes::value_type v1;
1035  Local_ac.insert(v1(Unique::get(),
1036  new Vnode<Dim>((*p).getDomain(),mynode,(*p).getVnode())));
1037  }
1038 
1039  // Everybody broadcasts their new local domains to everybody.
1040  // Build a message with the local domains and the id's for them.
1041  Message* bcast_mess = new Message;
1042  int count = idx_end - idx_begin;
1043  bcast_mess->put(count);
1044  // ::putMessage(*bcast_mess,idx_begin,idx_end);
1045  const Vnode<Dim>* indxmsg = idx_begin;
1046  while (indxmsg != idx_end)
1047  (indxmsg++)->putMessage(*bcast_mess);
1048  // Send it to everybody except yourself, and record the number sent.
1049  int tag = Ippl::Comm->next_tag( F_REPARTITION_BCAST_TAG , F_TAG_CYCLE );
1050  int node_count = Ippl::Comm->broadcast_others(bcast_mess,tag);
1051  // Create the container for the remote vnodes.
1052  ac_domain_vnodes *remote_ac = new ac_domain_vnodes( Domain );
1053  typedef typename ac_gc_domain_vnodes::value_type v1;
1054  Remotes_ac.insert(v1(gc0(),remote_ac) );
1055 
1056  // Loop until we receive a message from each other node.
1057  while ((--node_count)>=0)
1058  {
1059  // Receive a broadcast message from any node.
1060  int other_node = COMM_ANY_NODE;
1061  Message *recv_mess = Ippl::Comm->receive_block(other_node,tag);
1062  PAssert(recv_mess);
1063  // Extract the number of vnodes coming in.
1064  int count;
1065  recv_mess->get(count);
1066  // Now get the domains for the vnodes from the message.
1067  Vnode<Dim> p;
1068  while (count-- > 0) {
1069  p.getMessage(*recv_mess);
1070  Vnode<Dim> *vnode =
1071  new Vnode<Dim>(p.getDomain(), other_node, p.getVnode());
1072  // For domains of extent 1 in all directions, must call
1073  // DomainMap::insert() with the noSplit flag set to true, to avoid an
1074  // assertion failure in Index::split. Also do this for domains of
1075  // extent ZERO. These are not supposed to happen in IPPL, but when the
1076  // do when something like BinaryRepartition creates them, accomodating
1077  // them here will prevent subsequent code hangs or other errors and
1078  // allow recovery such as rejection of the result of BinaryRepartition
1079  // (discard the FieldLayout with zero-size domains and go on with the
1080  // original layout unchanged). (tjw)
1081  // if (p.getDomain().size() == 1) {
1082  if (p.getDomain().size() <= 1) {
1083  typedef typename ac_domain_vnodes::value_type v1;
1084  remote_ac->insert(v1(p.getDomain(), vnode), true);
1085  } else {
1086  typedef typename ac_domain_vnodes::value_type v1;
1087  remote_ac->insert(v1(p.getDomain(), vnode), false);
1088  }
1089  }
1090  }
1091 
1092  // Calculate the widths of all the vnodes
1093  calcWidths();
1094 }
1095 
1096 
1097 //----------------------------------------------------------------------
1098 //
1099 // Completely repartition this FieldLayout and all of the Fields
1100 // defined on it.
1101 //
1102 template<unsigned Dim>
1103 void
1104 FieldLayout<Dim>::Repartition(const NDIndex<Dim>* idxBegin,
1105  const NDIndex<Dim>* idxEnd)
1106 {
1107 
1108 
1109 
1110  // Build a temporary FieldLayout to declare the temporary arrays.
1111  FieldLayout<Dim> tempLayout(Domain,idxBegin,idxEnd);
1112 
1113  // (TJW) Copy the vnodesPerDirection information:
1114  if (vnodesPerDirection_m != 0) {
1115  tempLayout.vnodesPerDirection_m = new unsigned[Dim];
1116  for (unsigned int d=0; d<Dim; d++) {
1117  tempLayout.vnodesPerDirection_m[d] = vnodesPerDirection_m[d];
1118  }
1119  }
1120 
1121  // Repartition each field.
1122  iterator_if p, endp=end_if();
1123  for (p=begin_if(); p!=endp; ++p) {
1124  FieldLayoutUser *user = (FieldLayoutUser *)((*p).second);
1125  user->Repartition(&tempLayout);
1126  }
1127 
1128  // Copy back the layout information.
1129  Local_ac = tempLayout.Local_ac;
1130  Remotes_ac = tempLayout.Remotes_ac;
1131 
1132  // Calculate the widths of all the vnodes
1133  calcWidths();
1134 
1135  //INCIPPLSTAT(incRepartitions);
1136 }
1137 
1138 
1139 //----------------------------------------------------------------------
1140 //
1141 // Completely repartition this FieldLayout and all of the Fields
1142 // defined on it.
1143 // This differs from the previous ctor in that it allows preservation of
1144 // global Vnode integer ID numbers associated with the input Vnodes. --tjw
1145 //
1146 template<unsigned Dim>
1147 void
1148 FieldLayout<Dim>::Repartition(const Vnode<Dim>* idxBegin,
1149  const Vnode<Dim>* idxEnd)
1150 {
1151 
1152 
1153 
1154  // Build a temporary FieldLayout to declare the temporary arrays.
1155  FieldLayout<Dim> tempLayout(Domain,idxBegin,idxEnd);
1156 
1157  // (TJW) Copy the vnodesPerDirection information:
1158  if (vnodesPerDirection_m != 0) {
1159  tempLayout.vnodesPerDirection_m = new unsigned[Dim];
1160  for (unsigned int d=0; d<Dim; d++) {
1161  tempLayout.vnodesPerDirection_m[d] = vnodesPerDirection_m[d];
1162  }
1163  }
1164 
1165  // Repartition each field.
1166  iterator_if p, endp=end_if();
1167  for (p=begin_if(); p!=endp; ++p) {
1168  FieldLayoutUser *user = (FieldLayoutUser *)((*p).second);
1169  user->Repartition(&tempLayout);
1170  }
1171 
1172  // Copy back the layout information.
1173  Local_ac = tempLayout.Local_ac;
1174  Remotes_ac = tempLayout.Remotes_ac;
1175 
1176  // Calculate the widths of all the vnodes
1177  calcWidths();
1178 }
1179 
1180 
1181 //----------------------------------------------------------------------
1182 //
1183 // Write out a FieldLayout data to a file. The file is in ascii format,
1184 // with keyword=value. Return success.
1185 template<unsigned Dim>
1186 bool FieldLayout<Dim>::write(const char *filename) {
1187 
1188 
1189 
1190  unsigned int d;
1191 
1192  // only do the read on node 0
1193  if (Ippl::myNode() == 0) {
1194  // create the file, make sure the creation is OK
1195  FILE *f = fopen(filename, "w");
1196  if (f == 0) {
1197  ERRORMSG("Could not create FieldLayout data file '" << filename);
1198  ERRORMSG("'." << endl);
1199  return false;
1200  }
1201 
1202  // write out FieldLayout information
1203  fprintf(f, "dim = %d\n", Dim);
1204  fprintf(f, "vnodes = %d\n", numVnodes());
1205  fprintf(f, "pnodes = %d\n", Ippl::getNodes());
1206 
1207  fprintf(f, "domain =");
1208  for (d=0; d < Dim; ++d)
1209  fprintf(f, " %d %d %d",
1210  Domain[d].first(), Domain[d].last(), Domain[d].stride());
1211  fprintf(f, "\n");
1212 
1213  if (vnodesPerDirection_m != 0) {
1214  fprintf(f, "vnodesperdir =");
1215  for (d=0; d < Dim; ++d)
1216  fprintf(f, " %d", vnodesPerDirection_m[d]);
1217  fprintf(f, "\n");
1218  }
1219 
1220  for (iterator_iv localv = begin_iv(); localv != end_iv(); ++localv) {
1221  Vnode<Dim> *vn = (*localv).second.get();
1222  fprintf(f, "block = %d %d", vn->getVnode(), vn->getNode());
1223  for (d=0; d < Dim; ++d)
1224  fprintf(f, " %d %d %d", vn->getDomain()[d].first(),
1225  vn->getDomain()[d].last(), vn->getDomain()[d].stride());
1226  fprintf(f, "\n");
1227  }
1228 
1229  for (iterator_dv remotev = begin_rdv(); remotev != end_rdv(); ++remotev) {
1230  Vnode<Dim> *vn = (*remotev).second;
1231  fprintf(f, "block = %d %d", vn->getVnode(), vn->getNode());
1232  for (d=0; d < Dim; ++d)
1233  fprintf(f, " %d %d %d", vn->getDomain()[d].first(),
1234  vn->getDomain()[d].last(), vn->getDomain()[d].stride());
1235  fprintf(f, "\n");
1236  }
1237 
1238  fclose(f);
1239  }
1240 
1241  return true;
1242 }
1243 
1244 
1245 //----------------------------------------------------------------------
1246 //
1247 // Read in FieldLayout data from a file. If the
1248 // file contains data for an equal number of nodes as we are running on,
1249 // then that vnode -> pnode mapping will be used. If the file does not
1250 // contain info for the same number of pnodes, the vnodes will be
1251 // distributed in some other manner.
1252 // Note that if this FieldLayout is initially empty (e.g., it has no
1253 // local or remote vnodes), the domain will be changed to that in the
1254 // file. But if we already have some vnodes, we can only repartition
1255 // if the domain matches that in the file, or if we do not have any
1256 // users. If an error occurs, return
1257 // false and leave our own layout unchanged.
1258 template<unsigned Dim>
1259 bool FieldLayout<Dim>::read(const char *filename) {
1260 
1261 
1262 
1263  // generate a tag to use for communication
1264  int tag = Ippl::Comm->next_tag(F_LAYOUT_IO_TAG, F_TAG_CYCLE);
1265 
1266  // storage for data read from file
1267  NDIndex<Dim> fdomain;
1268  Vnode<Dim> *vnlist = 0;
1269  unsigned *vnodesPerDir = 0;
1270  int fdim = 0;
1271  int fvnodes = 0;
1272  int fpnodes = 0;
1273  int vnodesread = 0;
1274  int ok = 1;
1275 
1276  // only do the read on node 0
1277  if (Ippl::myNode() == 0) {
1278  // read the file, make sure the read is OK
1279  DiscMeta f(filename);
1280 
1281  // make sure it is OK
1282  DiscMeta::iterator metaline = f.begin();
1283  for ( ; ok == 1 && metaline != f.end(); ++metaline) {
1284  // get number of tokens and list of tokens in the line
1285  int numtokens = (*metaline).second.first;
1286  std::string *tokens = (*metaline).second.second;
1287 
1288  // check first word
1289  if (tokens[0] == "dim") {
1290  if (fdim != 0) {
1291  ERRORMSG("Repeated 'dim' line in FieldLayout data file '");
1292  ERRORMSG(filename << "'." << endl);
1293  ok = 0;
1294  }
1295  fdim = atoi(tokens[1].c_str());
1296  if (fdim != Dim) {
1297  ERRORMSG("Mismatched dimension in FieldLayout data file '");
1298  ERRORMSG(filename << "'." << endl);
1299  ok = 0;
1300  }
1301  } else if (tokens[0] == "vnodes") {
1302  if (fvnodes != 0 || vnlist != 0) {
1303  ERRORMSG("Repeated 'vnodes' line in FieldLayout data file '");
1304  ERRORMSG(filename << "'." << endl);
1305  ok = 0;
1306  }
1307  fvnodes = atoi(tokens[1].c_str());
1308  if (fvnodes < 1) {
1309  ERRORMSG("Incorrect 'vnodes' line in FieldLayout data file '");
1310  ERRORMSG(filename << "'." << endl);
1311  ok = 0;
1312  } else {
1313  vnlist = new Vnode<Dim>[fvnodes];
1314  }
1315  } else if (tokens[0] == "pnodes") {
1316  if (fpnodes != 0) {
1317  ERRORMSG("Repeated 'pnodes' line in FieldLayout data file '");
1318  ERRORMSG(filename << "'." << endl);
1319  ok = 0;
1320  }
1321  fpnodes = atoi(tokens[1].c_str());
1322  if (fpnodes < 1) {
1323  ERRORMSG("Incorrect 'pnodes' line in FieldLayout data file '");
1324  ERRORMSG(filename << "'." << endl);
1325  ok = 0;
1326  }
1327  } else if (tokens[0] == "vnodesperdir") {
1328  if (vnodesPerDir != 0) {
1329  ERRORMSG("Repeated 'vnodesperdir' line in FieldLayout data file '");
1330  ERRORMSG(filename << "'." << endl);
1331  ok = 0;
1332  }
1333  if (numtokens - 1 != Dim) {
1334  ERRORMSG("Incorrect 'vnodesperdir' line in FieldLayout data file '");
1335  ERRORMSG(filename << "'." << endl);
1336  ok = 0;
1337  }
1338  vnodesPerDir = new unsigned[Dim];
1339  for (unsigned int d=0; d < Dim; ++d) {
1340  vnodesPerDir[d] = atoi(tokens[d+1].c_str());
1341  if (vnodesPerDir[d] < 1) {
1342  ERRORMSG("Illegal vnode per direction value for dim=" << d);
1343  ERRORMSG(" in FieldLayout data file '" << filename << "'."<<endl);
1344  ok = 0;
1345  }
1346  }
1347  } else if (tokens[0] == "domain") {
1348  // make sure we have (first,last,stride) for all dims
1349  if ((numtokens-1) % 3 != 0 || (numtokens-1) / 3 != Dim) {
1350  ERRORMSG("Incorrect 'domain' line in FieldLayout data file '");
1351  ERRORMSG(filename << "'." << endl);
1352  ok = 0;
1353  }
1354  for (unsigned int d=0, dindx=1; d < Dim; ++d) {
1355  fdomain[d] = Index(atoi(tokens[dindx].c_str()),
1356  atoi(tokens[dindx + 1].c_str()),
1357  atoi(tokens[dindx + 2].c_str()));
1358  dindx += 3;
1359  }
1360  } else if (tokens[0] == "block") {
1361  // this is vnode information; it should be of the form
1362  // vnode# pnode# (domain, with 3 numbers/dimension)
1363  // So, the number of tokens should be a multipple of three
1364  if (numtokens % 3 != 0 || numtokens / 3 != (Dim+1)) {
1365  ERRORMSG("Incorrect 'block' line in FieldLayout data file '");
1366  ERRORMSG(filename << "'." << endl);
1367  ok = 0;
1368  } else if (vnlist == 0) {
1369  ERRORMSG("In FieldLayout data file '" << filename << ": You must ");
1370  ERRORMSG("give the number vnodes before any block lines." << endl);
1371  ok = 0;
1372  } else {
1373  int vnum = atoi(tokens[1].c_str());
1374  int pnum = atoi(tokens[2].c_str());
1375  if (pnum < 0 || pnum >= fpnodes) {
1376  ERRORMSG("Illegal pnode number in 'block' line ");
1377  ERRORMSG("of file '" << filename << "'." << endl);
1378  ok = 0;
1379  }
1380  NDIndex<Dim> vdom;
1381  for (unsigned int d=0, dindx=3; d < Dim; ++d) {
1382  vdom[d] = Index(atoi(tokens[dindx].c_str()),
1383  atoi(tokens[dindx + 1].c_str()),
1384  atoi(tokens[dindx + 2].c_str()));
1385  dindx += 3;
1386  }
1387 
1388  // construct a new vnode from this info, and store it at the
1389  // index given by the vnum value
1390  vnlist[vnodesread++] = Vnode<Dim>(vdom, pnum, vnum);
1391  }
1392  } else {
1393  ERRORMSG("Unrecognized '" << tokens[0] << "' line in FieldLayout ");
1394  ERRORMSG(" data file '" << filename << "'." << endl);
1395  ok = 0;
1396  }
1397  }
1398 
1399  // do some final sanity checks
1400  if (ok != 0 && (fvnodes < 1 || fdim != Dim || vnodesread != fvnodes)) {
1401  ERRORMSG("Inconsistent FieldLayout data in file '" << filename);
1402  ERRORMSG("'." << endl);
1403  ok = 0;
1404  }
1405 
1406  // if we have users, we cannot change the global domain any
1407  if ( ok != 0 && getNumUsers() > 0 && !(fdomain == Domain) ) {
1408  ERRORMSG("You cannot change the global domain of a FieldLayout which");
1409  ERRORMSG(" has users." << endl);
1410  ok = 0;
1411  }
1412 
1413  // Send out new layout info to other nodes
1414  if (Ippl::getNodes() > 1) {
1415  Message *msg = new Message;
1416  msg->put(ok);
1417  if (ok != 0) {
1418  msg->put(fvnodes);
1419  msg->put(fpnodes);
1420  msg->put(fdomain);
1421  int vnpdok = (vnodesPerDir != 0 ? 1 : 0);
1422  msg->put(vnpdok);
1423  if (vnpdok != 0)
1424  msg->put(vnodesPerDir, vnodesPerDir + Dim);
1425  for (int v=0; v < fvnodes; ++v) {
1426  msg->put(vnlist[v]);
1427  }
1428  }
1429 
1430  Ippl::Comm->broadcast_others(msg, tag);
1431  }
1432  } else {
1433  // on the client nodes, get the info from node 0 and store it
1434  // first receive the message
1435  int node = 0;
1436  Message *msg = Ippl::Comm->receive_block(node, tag);
1437  PAssert(msg);
1438 
1439  // then get data out of the message
1440  msg->get(ok);
1441  if (ok != 0) {
1442  msg->get(fvnodes);
1443  msg->get(fpnodes);
1444  msg->get(fdomain);
1445  int vnpdok;
1446  msg->get(vnpdok);
1447  if (vnpdok != 0) {
1448  vnodesPerDir = new unsigned[Dim];
1449  msg->get_iter(vnodesPerDir);
1450  }
1451  if ((vnodesread = fvnodes) > 0) {
1452  vnlist = new Vnode<Dim>[fvnodes];
1453  for (int v=0; v < fvnodes; ++v) {
1454  Vnode<Dim> vn;
1455  msg->get(vn);
1456  vnlist[v] = vn;
1457  }
1458  }
1459  }
1460 
1461  // done with the message, we can delete it
1462  delete msg;
1463  }
1464 
1465  // now, on all nodes, figure out which vnodes are ours, and which are
1466  // local. But if an error occurred during reading, we just exit
1467  // without changing anything.
1468  if (ok != 1 || fvnodes < 1) {
1469  if (vnodesPerDir != 0)
1470  delete [] vnodesPerDir;
1471  if (vnlist != 0)
1472  delete [] vnlist;
1473  return false;
1474  }
1475 
1476  // save the new domain
1477  Domain = fdomain;
1478 
1479  // for each vnode, figure out which physical node it belongs on, and
1480  // sort those vnodes to the top of the list.
1481  int localvnodes = 0;
1482  for (int v=0; v < fvnodes; ++v) {
1483  int node = vnlist[v].getNode() % Ippl::getNodes();
1484  if (node == Ippl::myNode()) {
1485  if (v > localvnodes) {
1486  // move this vnode to the top
1487  Vnode<Dim> tempv(vnlist[localvnodes]);
1488  vnlist[localvnodes] = vnlist[v];
1489  vnlist[v] = tempv;
1490  }
1491  localvnodes++;
1492  }
1493  }
1494 
1495  // save the info on how many vnodes we have per direction
1496  if (vnodesPerDirection_m != 0)
1497  delete [] vnodesPerDirection_m;
1498  vnodesPerDirection_m = vnodesPerDir;
1499 
1500  // Repartition the system using our list of local vnodes
1501  Repartition(vnlist, vnlist + localvnodes);
1502 
1503  // success!
1504  delete [] vnlist;
1505  return true;
1506 }
1507 
1508 
1509 //----------------------------------------------------------------------
1510 //
1511 // calculate the minimum vnode sizes in each dimension
1512 template<unsigned Dim>
1513 void FieldLayout<Dim>::calcWidths() {
1514  unsigned int d;
1515 
1516  // initialize widths first
1517  for (d=0; d < Dim; ++d)
1518  MinWidth[d] = getDomain()[d].length();
1519 
1520  // look for minimum width in local vnodes
1521  for (const_iterator_iv v_i = begin_iv() ; v_i != end_iv(); ++v_i) {
1522  const NDIndex<Dim> &dom = (*v_i).second->getDomain();
1523  for (d=0; d < Dim; ++d) {
1524  if ((unsigned int) dom[d].length() < MinWidth[d])
1525  MinWidth[d] = dom[d].length();
1526  }
1527  }
1528 
1529  // look for minimum width in remove vnodes
1530  ac_domain_vnodes *v_ac = Remotes_ac[ gc0() ].get();
1531  for (iterator_dv dv_i = v_ac->begin(); dv_i != v_ac->end(); ++ dv_i) {
1532  const NDIndex<Dim> &dom = (*dv_i).first;
1533  for (d=0; d < Dim; ++d) {
1534  if ((unsigned int) dom[d].length() < MinWidth[d])
1535  MinWidth[d] = dom[d].length();
1536  }
1537  }
1538 
1539  //Inform dbgmsg("calcWidths", INFORM_ALL_NODES);
1540  //dbgmsg << "Calculated minimum widths for layout " << *this << endl;
1541  //for (d=0; d < Dim; ++d) {
1542  // dbgmsg << " ==> MinWidth[" << d << "] = " << MinWidth[d];
1543  // dbgmsg << " ... " << (getDistribution(d)==SERIAL?"serial":"parallel");
1544  // dbgmsg << endl;
1545  //}
1546 }
1547 
1548 
1549 //----------------------------------------------------------------------
1550 //
1551 // Tell the FieldLayout that a FieldLayoutUser is using it
1552 template<unsigned Dim>
1553 void
1554 FieldLayout<Dim>::checkin(FieldLayoutUser& f,
1555  const GuardCellSizes<Dim>& gc)
1556 {
1557 
1558 
1559 
1560  checkinUser(f);
1561  iterator_gdv guarded = Remotes_ac.find(gc);
1562  if ( guarded == Remotes_ac.end() )
1563  new_gc_layout(gc);
1564 }
1565 
1566 //----------------------------------------------------------------------
1567 //
1568 // Tell the FieldLayout that a Field is no longer using it.
1569 template<unsigned Dim>
1570 void
1571 FieldLayout<Dim>::checkout(FieldLayoutUser& f)
1572 {
1573 
1574 
1575  checkoutUser(f);
1576 }
1577 
1578 template<unsigned Dim>
1579 NDIndex<Dim> FieldLayout<Dim>::getLocalNDIndex()
1580 {
1581 
1582  NDIndex<Dim> theId;
1583  for (iterator_iv localv = begin_iv(); localv != end_iv(); ++localv) {
1584  Vnode<Dim> *vn = (*localv).second.get();
1585  if(vn->getNode() == Ippl::myNode())
1586  theId = vn->getDomain();
1587  }
1588  return theId;
1589 }
1590 
1591 
1592 //---------------------------------------------------------------------
1593 // output
1594 //---------------------------------------------------------------------
1595 template<unsigned Dim>
1596 void FieldLayout<Dim>::write(std::ostream& out) const
1597 {
1598 
1599 
1600  int icount;
1601 
1602  // the whole domain, and the number of users
1603  out << "Domain = " << Domain << "\n";
1604  out << "FieldLayoutUsers = " << size_if() << "\n";
1605 
1606  // If applicable, vnodes per direction (tjw):
1607  if (vnodesPerDirection_m != 0) {
1608  out << "vnodesPerDirection_m[] =";
1609  for (unsigned int d=0; d<Dim; d++)
1610  out << " " << vnodesPerDirection_m[d];
1611  out << std::endl;
1612  }
1613 
1614  // iterate over the local vnodes and print them out.
1615  out << "Total number of vnodes = " << numVnodes() << std::endl;
1616  out << "Local Vnodes = " << Local_ac.size() << "\n";
1617  icount = 0;
1618  //tjw can now do operator<<() with whole Vnode objects
1619  //tjw for(const_iterator_iv v_i = begin_iv() ; v_i != end_iv(); ++v_i)
1620  //tjw out << " vnode " << icount++ <<" : "<< (*v_i).second->getDomain() << "\n";
1621  for(const_iterator_iv v_i = begin_iv() ; v_i != end_iv(); ++v_i)
1622  out << " vnode " << icount++ << ": " << *((*v_i).second) << std::endl;
1623 
1624  // iterate over the remote vnodes and print them out.
1625  ac_domain_vnodes *v_ac = Remotes_ac[ gc0() ].get();
1626  out << "Remote Vnodes = " << v_ac->size() << "\n";
1627  icount = 0;
1628  for (iterator_dv dv_i = v_ac->begin(); dv_i != v_ac->end(); ++ dv_i)
1629  out << " vnode " << icount++ << " : " << *((*dv_i).second) << std::endl;
1630 }
1631 
1632 
1633 /***************************************************************************
1634  * $RCSfile: FieldLayout.cpp,v $ $Author: adelmann $
1635  * $Revision: 1.1.1.1 $ $Date: 2003/01/23 07:40:27 $
1636  * IPPL_VERSION_ID: $Id: FieldLayout.cpp,v 1.1.1.1 2003/01/23 07:40:27 adelmann Exp $
1637  ***************************************************************************/
NDIndex::split
bool split(NDIndex< Dim > &l, NDIndex< Dim > &r, unsigned d, double a) const
Definition: NDIndexInlines.h:128
IpplInfo::getNodes
static int getNodes()
Definition: IpplInfo.cpp:773
Unique::get
static type get()
Definition: Unique.h:33
vnodeRecursiveBisection
void vnodeRecursiveBisection(int dim, const int *sizes, int nprocs, int *procs)
Definition: VRB.cpp:69
FieldLayout::Remotes_ac
ac_gc_domain_vnodes Remotes_ac
Definition: FieldLayout.h:441
FieldLayout::calcWidths
void calcWidths()
Definition: FieldLayout.hpp:1513
vmap::insert
std::pair< iterator, bool > insert(const value_type &x)
Definition: vmap.hpp:73
FieldLayout::new_gc_layout
void new_gc_layout(const GuardCellSizes< Dim > &)
Definition: FieldLayout.hpp:906
FieldLayout::iterator_gdv
ac_gc_domain_vnodes::iterator iterator_gdv
Definition: FieldLayout.h:78
INFORM_ALL_NODES
#define INFORM_ALL_NODES
Definition: Inform.h:38
Communicate::myNode
int myNode() const
Definition: Communicate.h:155
max
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:123
ERRORMSG
#define ERRORMSG(msg)
Definition: IpplInfo.h:399
NDIndex::getMessage
Message & getMessage(Message &m)
Definition: NDIndex.h:147
FieldLayout::gc0
static GuardCellSizes< Dim > gc0()
Definition: FieldLayout.h:84
COMM_ANY_NODE
const int COMM_ANY_NODE
Definition: Communicate.h:40
IpplInfo::myNode
static int myNode()
Definition: IpplInfo.cpp:794
DiscMeta::end
iterator end()
Definition: DiscMeta.h:77
NDIndex::size
unsigned size() const
Definition: NDIndexInlines.h:68
FieldLayout::initialize
void initialize(const Index &i1, e_dim_tag p1=PARALLEL, int vnodes=-1)
Definition: FieldLayout.hpp:119
DomainMap::end
iterator end()
Definition: DomainMap.h:502
AddGuardCells
NDIndex< Dim > AddGuardCells(const NDIndex< Dim > &idx, const GuardCellSizes< Dim > &g)
Definition: GuardCellSizes.h:69
FieldLayout::vnodesPerDirection_m
unsigned * vnodesPerDirection_m
Definition: FieldLayout.h:454
Message::get_iter
Message & get_iter(OutputIterator o)
Definition: Message.h:519
TagMaker::next_tag
int next_tag(int t, int s=1000)
Definition: TagMaker.h:43
FieldLayout::Repartition
void Repartition(const NDIndex< Dim > *, const NDIndex< Dim > *)
Definition: FieldLayout.hpp:1104
FieldLayout::RequestedLayout
e_dim_tag RequestedLayout[Dim]
Definition: FieldLayout.h:449
FieldLayout::Local_ac
ac_id_vnodes Local_ac
Definition: FieldLayout.h:440
DomainMap::insert
void insert(const value_type &d, bool noSplit=false)
Definition: DomainMap.hpp:43
FieldLayout::getVnodesPerDirection
unsigned getVnodesPerDirection(unsigned dir)
Definition: FieldLayout.hpp:891
FieldLayout::iterator_if
iterator_user iterator_if
Definition: FieldLayout.h:79
Vnode::getNode
int getNode() const
Definition: Vnode.h:69
PAssert_EQ
#define PAssert_EQ(a, b)
Definition: PAssert.h:119
Index
Definition: Index.h:236
FieldLayout::Domain
NDIndex< Dim > Domain
Definition: FieldLayout.h:444
DiscMeta::begin
iterator begin()
Definition: DiscMeta.h:76
FieldLayout::iterator_iv
ac_id_vnodes::iterator iterator_iv
Definition: FieldLayout.h:73
DomainMap::begin
iterator begin()
Definition: DomainMap.h:501
FieldLayout::ac_domain_vnodes
DomainMap< NDIndex< Dim >, RefCountedP< Vnode< Dim > >, Touches< Dim >, Contains< Dim >, Split< Dim > > ac_domain_vnodes
Definition: FieldLayout.h:68
FieldLayoutUser::Repartition
virtual void Repartition(UserList *)=0
Communicate::broadcast_others
virtual int broadcast_others(Message *, int, bool delmsg=true)
Definition: Communicate.cpp:414
FieldLayout::write
bool write(const char *filename)
Definition: FieldLayout.hpp:1186
FieldLayout::read
bool read(const char *filename)
Definition: FieldLayout.hpp:1259
vmap< GuardCellSizes< Dim >, my_auto_ptr< ac_domain_vnodes > >::value_type
std::pair< GuardCellSizes< Dim >, my_auto_ptr< ac_domain_vnodes > > value_type
Definition: vmap.h:71
Vnode::getMessage
Message & getMessage(Message &m)
Definition: Vnode.h:82
Vnode::getDomain
const NDIndex< Dim > & getDomain() const
Definition: Vnode.h:71
Message::get
Message & get(const T &cval)
Definition: Message.h:484
Message::put
Message & put(const T &val)
Definition: Message.h:414
FieldLayout::checkout
void checkout(FieldLayoutUser &f)
Definition: FieldLayout.hpp:1571
DiscMeta::iterator
container_t::iterator iterator
Definition: DiscMeta.h:41
Vnode
Definition: Vnode.h:22
PAssert
#define PAssert(c)
Definition: PAssert.h:117
FieldLayout::setup
void setup(const NDIndex< Dim > &, e_dim_tag *, int)
Definition: FieldLayout.hpp:523
e_dim_tag
e_dim_tag
Definition: FieldLayout.h:55
PInsist
#define PInsist(c, m)
Definition: PAssert.h:135
Dim
const unsigned Dim
Definition: P3MPoissonSolver.h:23
putMessage
void putMessage(Message &m, const T &t)
Definition: Message.h:557
Communicate::receive_block
Message * receive_block(int &node, int &tag)
Definition: Communicate.cpp:301
Inform
Definition: Inform.h:41
DomainMap::size
size_type size() const
Definition: DomainMap.h:514
IpplInfo::Comm
static Communicate * Comm
Definition: IpplInfo.h:93
FieldLayout::const_iterator_iv
ac_id_vnodes::const_iterator const_iterator_iv
Definition: FieldLayout.h:74
FieldLayout::getLocalNDIndex
NDIndex< Dim > getLocalNDIndex()
Definition: FieldLayout.hpp:1579
FieldLayout::~FieldLayout
virtual ~FieldLayout()
Definition: FieldLayout.hpp:90
FieldLayout::checkin
void checkin(FieldLayoutUser &f, const GuardCellSizes< Dim > &gc=gc0())
Definition: FieldLayout.hpp:1554
Vnode::getVnode
int getVnode() const
Definition: Vnode.h:70
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
F_LAYOUT_IO_TAG
#define F_LAYOUT_IO_TAG
Definition: Tags.h:52
F_REPARTITION_BCAST_TAG
#define F_REPARTITION_BCAST_TAG
Definition: Tags.h:48
F_TAG_CYCLE
#define F_TAG_CYCLE
Definition: Tags.h:53
Generated on Sat Jan 4 2020 14:47:36 for OPAL (Object Oriented Parallel Accelerator Library) by   doxygen 1.8.5