OPAL (Object Oriented Parallel Accelerator Library) 2022.1
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)
43template<unsigned Dim>
44FieldLayout<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.
66template<unsigned Dim>
67FieldLayout<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.
89template<unsigned Dim>
90FieldLayout<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
107template<unsigned Dim>
108void
109FieldLayout<Dim>::initialize(const NDIndex<Dim>& domain,
110 e_dim_tag *p, int vnodes) {
111
112
113 setup(domain, p, vnodes);
114}
115
116
117template<unsigned Dim>
118void
119FieldLayout<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
129template<unsigned Dim>
130void
131FieldLayout<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}
146template<unsigned Dim>
147void
148FieldLayout<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}
164template<unsigned Dim>
165void
166FieldLayout<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}
187template<unsigned Dim>
188void
189FieldLayout<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}
211template<unsigned Dim>
212void
213FieldLayout<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
245template<unsigned Dim>
246void
247FieldLayout<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
268template<unsigned Dim>
269void
270FieldLayout<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
281template<unsigned Dim>
282void
283FieldLayout<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}
301template<unsigned Dim>
302void
303FieldLayout<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}
325template<unsigned Dim>
326void
327FieldLayout<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}
354template<unsigned Dim>
355void
356FieldLayout<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}
387template<unsigned Dim>
388void
389FieldLayout<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
434template<unsigned Dim>
435void
436FieldLayout<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
521template<unsigned Dim>
522void
523FieldLayout<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
713template<unsigned Dim>
714void
715FieldLayout<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.
889template<unsigned Dim>
890unsigned FieldLayout<Dim>::
891getVnodesPerDirection(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
905template<unsigned Dim>
906void 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//----------------------------------------------------------------------
932template<unsigned Dim>
933FieldLayout<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
1016template<unsigned Dim>
1017FieldLayout<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//
1102template<unsigned Dim>
1103void
1104FieldLayout<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//
1146template<unsigned Dim>
1147void
1148FieldLayout<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.
1185template<unsigned Dim>
1186bool 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.
1258template<unsigned Dim>
1259bool 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
1512template<unsigned Dim>
1513void 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
1552template<unsigned Dim>
1553void
1554FieldLayout<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.
1569template<unsigned Dim>
1570void
1571FieldLayout<Dim>::checkout(FieldLayoutUser& f)
1572{
1573
1574
1575 checkoutUser(f);
1576}
1577
1578template<unsigned Dim>
1579NDIndex<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//---------------------------------------------------------------------
1595template<unsigned Dim>
1596void 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 ***************************************************************************/
Dim
const unsigned Dim
Definition: P3MPoissonSolver.h:26
e_dim_tag
e_dim_tag
Definition: FieldLayout.h:55
PARALLEL
@ PARALLEL
Definition: FieldLayout.h:55
SERIAL
@ SERIAL
Definition: FieldLayout.h:55
vnodeRecursiveBisection
void vnodeRecursiveBisection(int dim, const int *sizes, int nprocs, int *procs)
Definition: VRB.cpp:69
putMessage
void putMessage(Message &m, const T &t)
Definition: Message.h:549
COMM_ANY_NODE
const int COMM_ANY_NODE
Definition: Communicate.h:40
F_REPARTITION_BCAST_TAG
#define F_REPARTITION_BCAST_TAG
Definition: Tags.h:48
F_TAG_CYCLE
#define F_TAG_CYCLE
Definition: Tags.h:53
F_LAYOUT_IO_TAG
#define F_LAYOUT_IO_TAG
Definition: Tags.h:52
max
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:84
AddGuardCells
NDIndex< Dim > AddGuardCells(const NDIndex< Dim > &idx, const GuardCellSizes< Dim > &g)
Definition: GuardCellSizes.h:69
a
std::complex< double > a
Definition: IpplExpressions.h:56
endl
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
INFORM_ALL_NODES
#define INFORM_ALL_NODES
Definition: Inform.h:39
PInsist
#define PInsist(c, m)
Definition: PAssert.h:120
PAssert_EQ
#define PAssert_EQ(a, b)
Definition: PAssert.h:104
PAssert
#define PAssert(c)
Definition: PAssert.h:102
ERRORMSG
#define ERRORMSG(msg)
Definition: IpplInfo.h:350
NDIndex::size
unsigned size() const
Definition: NDIndexInlines.h:68
NDIndex::getMessage
Message & getMessage(Message &m)
Definition: NDIndex.h:138
NDIndex::split
bool split(NDIndex< Dim > &l, NDIndex< Dim > &r, unsigned d, double a) const
Definition: NDIndexInlines.h:128
DomainMap::end
iterator end()
Definition: DomainMap.h:497
DomainMap::size
size_type size() const
Definition: DomainMap.h:509
DomainMap::insert
void insert(const value_type &d, bool noSplit=false)
Definition: DomainMap.hpp:42
DomainMap::begin
iterator begin()
Definition: DomainMap.h:496
FieldLayout::write
bool write(const char *filename)
Definition: FieldLayout.hpp:1186
FieldLayout::Domain
NDIndex< Dim > Domain
Definition: FieldLayout.h:444
FieldLayout::new_gc_layout
void new_gc_layout(const GuardCellSizes< Dim > &)
Definition: FieldLayout.hpp:906
FieldLayout::Repartition
void Repartition(const NDIndex< Dim > *, const NDIndex< Dim > *)
Definition: FieldLayout.hpp:1104
FieldLayout::read
bool read(const char *filename)
Definition: FieldLayout.hpp:1259
FieldLayout::getLocalNDIndex
NDIndex< Dim > getLocalNDIndex()
Definition: FieldLayout.hpp:1579
FieldLayout::Remotes_ac
ac_gc_domain_vnodes Remotes_ac
Definition: FieldLayout.h:441
FieldLayout::checkout
void checkout(FieldLayoutUser &f)
Definition: FieldLayout.hpp:1571
FieldLayout::getVnodesPerDirection
unsigned getVnodesPerDirection(unsigned dir)
Definition: FieldLayout.hpp:891
FieldLayout::iterator_gdv
ac_gc_domain_vnodes::iterator iterator_gdv
Definition: FieldLayout.h:78
FieldLayout::calcWidths
void calcWidths()
Definition: FieldLayout.hpp:1513
FieldLayout::const_iterator_iv
ac_id_vnodes::const_iterator const_iterator_iv
Definition: FieldLayout.h:74
FieldLayout::gc0
static GuardCellSizes< Dim > gc0()
Definition: FieldLayout.h:84
FieldLayout::RequestedLayout
e_dim_tag RequestedLayout[Dim]
Definition: FieldLayout.h:449
FieldLayout::~FieldLayout
virtual ~FieldLayout()
Definition: FieldLayout.hpp:90
FieldLayout::vnodesPerDirection_m
unsigned * vnodesPerDirection_m
Definition: FieldLayout.h:454
FieldLayout::ac_domain_vnodes
DomainMap< NDIndex< Dim >, RefCountedP< Vnode< Dim > >, Touches< Dim >, Contains< Dim >, Split< Dim > > ac_domain_vnodes
Definition: FieldLayout.h:68
FieldLayout::initialize
void initialize(const Index &i1, e_dim_tag p1=PARALLEL, int vnodes=-1)
Definition: FieldLayout.hpp:119
FieldLayout::Local_ac
ac_id_vnodes Local_ac
Definition: FieldLayout.h:440
FieldLayout::setup
void setup(const NDIndex< Dim > &, e_dim_tag *, int)
Definition: FieldLayout.hpp:523
FieldLayout::iterator_if
iterator_user iterator_if
Definition: FieldLayout.h:79
FieldLayout::checkin
void checkin(FieldLayoutUser &f, const GuardCellSizes< Dim > &gc=gc0())
Definition: FieldLayout.hpp:1554
FieldLayout::iterator_iv
ac_id_vnodes::iterator iterator_iv
Definition: FieldLayout.h:73
FieldLayoutUser::Repartition
virtual void Repartition(UserList *)=0
Vnode
Definition: Vnode.h:38
Vnode::getVnode
int getVnode() const
Definition: Vnode.h:66
Vnode::getDomain
const NDIndex< Dim > & getDomain() const
Definition: Vnode.h:67
Vnode::getMessage
Message & getMessage(Message &m)
Definition: Vnode.h:78
Vnode::getNode
int getNode() const
Definition: Vnode.h:65
Index
Definition: Index.h:237
Communicate::broadcast_others
virtual int broadcast_others(Message *, int, bool delmsg=true)
Definition: Communicate.cpp:413
Communicate::receive_block
Message * receive_block(int &node, int &tag)
Definition: Communicate.cpp:300
Communicate::myNode
int myNode() const
Definition: Communicate.h:155
Message::get_iter
Message & get_iter(OutputIterator o)
Definition: Message.h:511
Message::put
Message & put(const T &val)
Definition: Message.h:406
Message::get
Message & get(const T &cval)
Definition: Message.h:476
TagMaker::next_tag
int next_tag(int t, int s=1000)
Definition: TagMaker.h:39
DiscMeta::begin
iterator begin()
Definition: DiscMeta.h:76
DiscMeta::end
iterator end()
Definition: DiscMeta.h:77
DiscMeta::iterator
container_t::iterator iterator
Definition: DiscMeta.h:41
Inform
Definition: Inform.h:42
IpplInfo::getNodes
static int getNodes()
Definition: IpplInfo.cpp:670
IpplInfo::myNode
static int myNode()
Definition: IpplInfo.cpp:691
IpplInfo::Comm
static Communicate * Comm
Definition: IpplInfo.h:84
Unique::get
static type get()
Definition: Unique.h:33
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:64
vmap::insert
std::pair< iterator, bool > insert(const value_type &x)
Definition: vmap.hpp:73
Generated on Thu Oct 20 2022 17:40:08 for OPAL (Object Oriented Parallel Accelerator Library) by doxygen 1.9.3