source: XIOS/dev/dev_trunk_omp/extern/remap/src/parallel_tree.cpp @ 1646

Last change on this file since 1646 was 1646, checked in by yushan, 5 years ago

branch merged with trunk @1645. arch file (ep&mpi) added for ADA

File size: 16.4 KB
Line 
1#include <cassert>
2#include "node.hpp"
3#include "timerRemap.hpp"
4#include "circle.hpp"
5#include "meshutil.hpp"
6#include "polyg.hpp"
7#include "gridRemap.hpp"
8#include "intersect.hpp"
9#include "errhandle.hpp"
10#include "mpi_routing.hpp"
11#include "misc.hpp"
12
13#include "parallel_tree.hpp"
14#ifdef _usingEP
15using namespace ep_lib;
16#endif
17
18namespace sphereRemap {
19
20extern CRemapGrid srcGrid;
21#pragma omp threadprivate(srcGrid)
22
23extern CRemapGrid tgtGrid;
24#pragma omp threadprivate(tgtGrid)
25
26static const int assignLevel = 2;
27
28// only the circle is packed, rest of node will be initialized on unpacking
29static void packNode(Node& node, char *buffer, int& index)
30{
31        if (buffer == NULL) // compute size only
32                index += 4 * sizeof(double);
33        else
34        {
35                *(Coord *)(&buffer[index]) = node.centre;
36                index += sizeof(Coord);
37
38                *(double *)(&buffer[index]) = node.radius;
39                index += sizeof(double);
40        }
41}
42
43static void unpackNode(Node& node, char *buffer, int& index)
44{
45        Coord centre;
46        double r;
47
48        centre = *(Coord *)(&buffer[index]);
49        index += sizeof(Coord);
50
51        r = *(double *)(&buffer[index]);
52        index += sizeof(double);
53
54        node.centre = centre;
55        node.radius = r;
56}
57
58
59static void packElement(Elt *ptElement, char *buffer, int& index)
60{
61        if (buffer == NULL) index += packedPolygonSize(*ptElement);
62        else packPolygon(*ptElement, buffer, index);
63}
64
65static void unpackElement(Elt *ptElement, char *buffer, int& index)
66{
67        unpackPolygon(*ptElement, buffer, index);
68}
69
70static void packVector(const vector<int>& vec, char *buf, int& pos)
71{
72        if (buf == NULL)
73                pos += sizeof(int) + vec.size()*sizeof(int);
74        else
75        {
76                *((int *) &(buf[pos])) = vec.size();
77                pos += sizeof(int);
78                for (int i = 0; i < vec.size(); i++)
79                {
80                        *((int *) &(buf[pos])) = vec[i];
81                        pos += sizeof(int);
82                }
83        }
84}
85
86static void unpackVector(vector<int>& vec, const char *buf, int& pos)
87{
88        vec.resize(*((int *) &(buf[pos])));
89        pos += sizeof(int);
90        for (int i = 0; i < vec.size(); i++)
91        {
92                vec[i] = *((int *) &(buf[pos]));
93                pos += sizeof(int);
94        }
95}
96
97static void assignRoute(CSampleTree& tree, const CCascadeLevel& cl)  // newroot <- root
98{
99        vector<int> routeRank(cl.group_size);
100        for (int i = 0; i < cl.group_size; i++)
101                routeRank[i] = i; //(cl.group_size + cl.polour() - i) % cl.group_size;
102        std::vector<int>::iterator rank = routeRank.begin();
103        tree.root->assignRoute(rank, assignLevel);
104}
105
106static void transferRoutedNodes(CMPIRouting& MPIRoute, /*const*/ vector<Node>& nodeSend, const vector<int>& route, vector<Node>& nodeRecv)
107{
108        /* `route` knows what goes where */
109        MPIRoute.init(route);
110        int nRecv = MPIRoute.getTotalSourceElement();
111        nodeRecv.resize(nRecv);
112        MPIRoute.transferToTarget(&nodeSend[0], &nodeRecv[0], packNode, unpackNode);
113}
114
115static void transferRoutedIntersections(CMPIRouting& MPIRoute, const vector<Node>& nodeSend, const vector<vector<int> >& route, vector<Node>& nodeRecv)
116{
117        // `route` knows what goes where
118        MPIRoute.init(route);
119        int nRecv = MPIRoute.getTotalSourceElement();
120        nodeRecv.resize(nRecv);
121        MPIRoute.transferToTarget((Node * /*mpi wants non-const*/) &nodeSend[0], &nodeRecv[0], packNode, unpackNode);
122//cout << MPIRoute.mpiRank << " ROUTE " << nRecv << ": " << nodeSend.size() << " " << nodeRecv.size() << "    ";
123}
124
125CParallelTree::CParallelTree(MPI_Comm comm) : communicator(comm), cascade(MIN_NODE_SZ*MIN_NODE_SZ*2, comm)
126{
127  treeCascade.reserve(cascade.num_levels);
128  for (int lev = 0; lev < cascade.num_levels; lev++)
129    treeCascade.push_back(CSampleTree(cascade.level[lev].group_size, assignLevel));
130}
131
132void CParallelTree::buildSampleTreeCascade(vector<Node>& sampleNodes /*route field will be modified*/, int level)
133{
134        buildSampleTree(treeCascade[level], sampleNodes, cascade.level[level]);
135        assignRoute(treeCascade[level] /*out*/, cascade.level[level] /*in*/);
136
137        if (level+1 < cascade.num_levels)
138        {
139                vector<int> route(sampleNodes.size());
140                treeCascade[level].routeNodes(route, sampleNodes, assignLevel);
141
142                vector<Node> routedNodes;
143                CMPIRouting mpiRoute(cascade.level[level].pg_comm);
144                transferRoutedNodes(mpiRoute, sampleNodes, route, routedNodes);
145                buildSampleTreeCascade(routedNodes, level+1);
146        }
147}
148
149void buildSampleTree(CSampleTree& tree, const vector<Node>& node, const CCascadeLevel& comm)
150/*
151        In the beginning all the sample elements are distributed
152        -> communicate to make available at each rank
153           so that each rank can build the same sample tree
154*/
155{
156        int n = node.size(); // number of samples intially on this rank (local)
157
158        int blocSize = comm.group_size * ipow(MAX_NODE_SZ, assignLevel);
159
160        int nrecv; // global number of samples  THIS WILL BE THE NUMBER OF LEAFS IN THE SAMPLE TREE
161        MPI_Allreduce(&n, &nrecv, 1, MPI_INT, MPI_SUM, comm.comm); // => size of sample tree does not depend on keepNodes!
162        double ratio = blocSize / (1.0 * nrecv);
163        int nsend = ratio * n + 1; // nsend = n_local_samples / n_global_samples * blocksize + 1 = blocksize/comm.size
164        if (nsend > n) nsend = n;
165
166        int *counts = new int[comm.size];
167        MPI_Allgather(&nsend, 1, MPI_INT, counts, 1, MPI_INT, comm.comm);
168
169        nrecv = 0;
170        int *displs = new int[comm.size];
171        for (int i = 0; i < comm.size; i++)
172        {
173                displs[i] = 4 * nrecv;
174                nrecv += counts[i];
175                counts[i] = 4 * counts[i];
176        }
177
178        /* pack circles around sample elements */
179        double *sendBuffer = new double[nsend*4];
180        int index = 0;
181        vector<int> randomArray(n);
182        randomizeArray(randomArray);
183        for (int i = 0; i < nsend; i++)
184        {
185                const Node& no = node[randomArray[i]];
186                *((Coord *) (sendBuffer + index)) = no.centre;
187                index += sizeof(Coord)/sizeof(*sendBuffer);
188                sendBuffer[index++] = no.radius;
189        }
190
191        /* each process needs the sample elements from all processes */
192        double *recvBuffer = new double[nrecv*4];
193        MPI_Allgatherv(sendBuffer, 4 * nsend, MPI_DOUBLE, recvBuffer, counts, displs, MPI_DOUBLE, comm.comm);
194        delete[] sendBuffer;
195        delete[] counts;
196        delete[] displs;
197
198        /* unpack */
199/*
200        randomArray.resize(nrecv);
201        randomizeArray(randomArray);
202        tree.leafs.resize(nrecv);
203        index = 0;
204
205
206  for (int i = 0; i < nrecv; i++)
207        {
208                Coord x = *(Coord *)(&recvBuffer[index]);
209                index += sizeof(Coord)/sizeof(*recvBuffer);
210                double radius = recvBuffer[index++];
211                tree.leafs[randomArray[i]].centre = x;
212                tree.leafs[randomArray[i]].radius = radius;
213
214        }
215*/
216
217  randomArray.resize(blocSize);
218        randomizeArray(randomArray);
219        tree.leafs.resize(blocSize);
220        index = 0;
221 
222  size_t s=(sizeof(Coord)/sizeof(*recvBuffer)+1)*nrecv ;
223 
224  for (int i = 0; i < blocSize; i++)
225        {
226                Coord x = *(Coord *)(&recvBuffer[index%s]);
227                index += sizeof(Coord)/sizeof(*recvBuffer);
228                double radius = recvBuffer[index%s];
229    index++ ;
230                tree.leafs[randomArray[i]].centre = x;
231                tree.leafs[randomArray[i]].radius = radius;
232
233        }
234
235
236        delete [] recvBuffer;
237
238        CTimer::get("buildSampleTree(local)").resume();
239        tree.build(tree.leafs);
240//      cout << "SampleTree build : assign Level " << assignLevel << " nb Nodes : " << tree.levelSize[assignLevel] << endl;
241        CTimer::get("buildSampleTree(local)").suspend();
242        CTimer::get("buildSampleTree(local)").print();
243
244        /* End gracefully if sample tree could not be constructed with desired number of nodes on assignLevel */
245        int allok, ok = (tree.levelSize[assignLevel] == comm.group_size);
246        if (!ok)
247  {
248    cerr << comm.rank << ": PROBLEM: (node assign)" << tree.levelSize[assignLevel] << " != " << comm.group_size << " (keepNodes)" 
249         << "   node size : "<<node.size()<<"   bloc size : "<<blocSize<<"  total number of leaf : "<<tree.leafs.size()<<endl ;
250/*
251        MPI_Allreduce(&ok, &allok, 1, MPI_INT, MPI_PROD, communicator);
252        if (!allok) {
253                MPI_Finalize();
254                exit(1);
255        }
256*/
257    MPI_Abort(MPI_COMM_WORLD,-1) ;
258  }
259/*
260  cout<<"*******************************************"<<endl ;
261  cout<<"******* Sample Tree output        *********"<<endl ;
262  cout<<"*******************************************"<<endl ;
263  tree.output(cout,1) ;
264
265  cout<<"*******************************************"<<endl ;
266*/
267
268  assert(tree.root->incluCheck() == 0);
269}
270
271
272void CParallelTree::buildLocalTree(const vector<Node>& node, const vector<int>& route)
273{
274        CMPIRouting MPIRoute(communicator);
275        MPI_Barrier(communicator);
276        CTimer::get("buildLocalTree(initRoute)").resume();
277        MPIRoute.init(route);
278        CTimer::get("buildLocalTree(initRoute)").suspend();
279        CTimer::get("buildLocalTree(initRoute)").print();
280
281        nbLocalElements = MPIRoute.getTotalSourceElement();
282        localElements = new Elt[nbLocalElements];
283
284        vector<Elt*> ptElement(node.size());
285        for (int i = 0; i < node.size(); i++)
286                ptElement[i] = (Elt *) (node[i].data);
287
288        vector<Elt*> ptLocalElement(nbLocalElements);
289        for (int i = 0; i < nbLocalElements; i++)
290                ptLocalElement[i] = &localElements[i];
291
292        CTimer::get("buildLocalTree(transfer)").resume();
293        MPIRoute.transferToTarget(&ptElement[0], &ptLocalElement[0], packElement, unpackElement);
294        CTimer::get("buildLocalTree(transfer)").suspend();
295        CTimer::get("buildLocalTree(transfer)").print();
296
297        CTimer::get("buildLocalTree(local)").resume();
298
299        int mpiRank;
300        MPI_Comm_rank(communicator, &mpiRank);
301        localTree.leafs.reserve(nbLocalElements);
302        for (int i = 0; i < nbLocalElements; i++)
303        {
304                Elt& elt = localElements[i];
305                elt.id.ind = i;
306                elt.id.rank = mpiRank;
307                localTree.leafs.push_back(Node(elt.x, cptRadius(elt), &localElements[i]));
308        }
309        localTree.build(localTree.leafs);
310
311        cptAllEltsGeom(localElements, nbLocalElements, srcGrid.pole);
312        CTimer::get("buildLocalTree(local)").suspend();
313        CTimer::get("buildLocalTree(local)").print();
314}
315
316void CParallelTree::build(vector<Node>& node, vector<Node>& node2)
317{
318
319        int assignLevel = 2;
320        int nbSampleNodes = 2*ipow(MAX_NODE_SZ + 1, assignLevel);
321
322
323  long int nb1, nb2, nb, nbTot ;
324  nb1=node.size() ; nb2=node2.size() ;
325  nb=nb1+nb2 ;
326  MPI_Allreduce(&nb, &nbTot, 1, MPI_LONG, MPI_SUM, communicator) ;
327  int commSize ;
328  MPI_Comm_size(communicator,&commSize) ;
329 
330        // make multiple of two
331        nbSampleNodes /= 2;
332        nbSampleNodes *= 2;
333//  assert( nbTot > nbSampleNodes*commSize) ;
334   
335  int nbSampleNodes1 = nbSampleNodes * (nb1*commSize)/(1.*nbTot) ;
336  int nbSampleNodes2 = nbSampleNodes * (nb2*commSize)/(1.*nbTot) ;
337 
338
339//      assert(node.size() > nbSampleNodes);
340//      assert(node2.size() > nbSampleNodes);
341//      assert(node.size() + node2.size() > nbSampleNodes);
342        vector<Node> sampleNodes; sampleNodes.reserve(nbSampleNodes1+nbSampleNodes2);
343
344        vector<int> randomArray1(node.size());
345        randomizeArray(randomArray1);
346        vector<int> randomArray2(node2.size());
347        randomizeArray(randomArray2);
348
349/*     
350        int s1,s2 ;
351       
352        if (node.size()< nbSampleNodes/2)
353        {
354          s1 = node.size() ;
355          s2 = nbSampleNodes-s1 ;
356        }
357        else if (node2.size()< nbSampleNodes/2)
358        {
359          s2 = node.size() ;
360          s1 = nbSampleNodes-s2 ;
361        }
362        else
363        {
364          s1=nbSampleNodes/2 ;
365          s2=nbSampleNodes/2 ;
366        }
367*/
368        for (int i = 0; i <nbSampleNodes1; i++) sampleNodes.push_back(Node(node[randomArray1[i%nb1]].centre,  node[randomArray1[i%nb1]].radius, NULL));
369        for (int i = 0; i <nbSampleNodes2; i++) sampleNodes.push_back(Node(node2[randomArray2[i%nb2]].centre, node2[randomArray2[i%nb2]].radius, NULL));
370
371/*         
372        for (int i = 0; i < nbSampleNodes/2; i++)
373        {
374          sampleNodes.push_back(Node(node[randomArray1[i]].centre,  node[randomArray1[i]].radius, NULL));
375          sampleNodes.push_back(Node(node2[randomArray2[i]].centre, node2[randomArray2[i]].radius, NULL));
376        }
377*/
378        CTimer::get("buildParallelSampleTree").resume();
379        //sampleTree.buildParallelSampleTree(sampleNodes, cascade);
380        buildSampleTreeCascade(sampleNodes);
381        CTimer::get("buildParallelSampleTree").suspend();
382        CTimer::get("buildParallelSampleTree").print();
383
384        //route source mesh
385        CTimer::get("parallelRouteNode").resume();
386        vector<int> route(node.size());
387        routeNodes(route /*out*/, node);
388        CTimer::get("parallelRouteNode").suspend();
389        CTimer::get("parallelRouteNode").print();
390
391        CTimer::get("buildLocalTree").resume();
392        buildLocalTree(node, route);
393        CTimer::get("buildLocalTree").suspend();
394        CTimer::get("buildLocalTree").print();
395
396        CTimer::get("buildRouteTree").resume();
397        /* update circles of tree cascade so it can be used for routing */
398        updateCirclesForRouting(localTree.root->centre, localTree.root->radius);
399        CTimer::get("buildRouteTree").suspend();
400        CTimer::get("buildRouteTree").print();
401}
402
403void CParallelTree::routeNodes(vector<int>& route, vector<Node>& nodes /*route field used*/, int level)
404{
405        treeCascade[level].routeNodes(route /*assign*/, nodes, assignLevel);
406
407        if (level+1 < cascade.num_levels)
408        {
409                vector<Node> routedNodes;
410                CMPIRouting MPIRoute(cascade.level[level].pg_comm);
411                transferRoutedNodes(MPIRoute, nodes, route /*use*/, routedNodes);
412                vector<int> globalRank(routedNodes.size());
413                routeNodes(globalRank, routedNodes, level + 1);
414                MPIRoute.transferFromSource(&route[0] /*override*/, &globalRank[0]);
415        }
416        else
417        {
418                CMPIRouting MPIRoute(cascade.level[level].comm); // or use pg_comm, on last cascade level identical
419                MPIRoute.init(route);
420                int nbRecvNode = MPIRoute.getTotalSourceElement();
421                vector<int> globalRank(nbRecvNode);
422                for (int i = 0; i < globalRank.size(); i++)
423                        globalRank[i] = cascade.level[0].rank;
424                MPIRoute.transferFromSource(&route[0] /*override*/, &globalRank[0]);
425        }
426}
427
428/* assume `to` to be empty vector at entry */
429void linearize(const vector<vector<int> >& from, vector<int>& to)
430{
431        int cnt = 0;
432        for (int i = 0; i < from.size(); ++i)
433                cnt += from[i].size();
434        to.resize(cnt);
435        vector<int>::iterator dest = to.begin();
436        for (int i = 0; i < from.size(); ++i)
437                dest = copy(from[i].begin(), from[i].end(), dest);
438}
439
440/* at entry `to` must already have it's final shape and only values are overriden */
441void delinearize(const vector<int>& from, vector<vector<int> >& to)
442{
443        vector<int>::const_iterator end, src = from.begin();
444        for (int i = 0; i < to.size(); ++i, src=end)
445                copy(src, end = src + to[i].size(), to[i].begin());
446}
447
448void CParallelTree::routeIntersections(vector<vector<int> >& routes, vector<Node>& nodes, int level)
449{
450        treeCascade[level].routeIntersections(routes /*assign*/, nodes);
451
452        if (level+1 < cascade.num_levels)
453        {
454                vector<Node> routedNodes;
455                CMPIRouting MPIRoute(cascade.level[level].pg_comm);
456
457                vector<int> flattenedRoutes1;
458                linearize(routes, flattenedRoutes1);
459                vector<Node> double_nodes(flattenedRoutes1.size());
460                int j = 0;
461                for (int i = 0; i < routes.size(); ++i)
462                        for (int k = 0; k < routes[i].size(); ++k, ++j)
463                                double_nodes[j] = nodes[i];
464                transferRoutedNodes(MPIRoute, double_nodes, flattenedRoutes1 /*use*/, routedNodes);
465                vector<vector<int> > globalRanks(routedNodes.size());
466                routeIntersections(globalRanks /*out*/, routedNodes /*in*/, level + 1);
467                vector<vector<int> > flattenedRoutes(flattenedRoutes1.size());
468                // transferFromSource expects sizes (nbSendNode=flattenedRoutes, nbRecvNode=routedNodes.size())
469                MPIRoute.transferFromSource(&flattenedRoutes[0], &globalRanks[0], packVector, unpackVector);
470                for (int i = 0, j = 0; i < routes.size(); ++i)
471                {
472                        int old_size = routes[i].size();
473                        routes[i].resize(0);
474                        for (int k = 0; k < old_size; ++k, ++j)
475                                for (int l = 0; l < flattenedRoutes[j].size(); ++l)
476                                        routes[i].push_back(flattenedRoutes[j][l]);
477                }
478                assert(j == flattenedRoutes1.size());
479
480        }
481        else
482        {
483                CMPIRouting MPIRoute(cascade.level[level].comm);
484                MPIRoute.init(routes);
485                int nbRecvNode = MPIRoute.getTotalSourceElement();
486                vector<int> globalRanks(nbRecvNode, cascade.level[0].rank);
487                vector<int> flattened_routes;
488                linearize(routes, flattened_routes);
489                MPIRoute.transferFromSource(&flattened_routes[0], &globalRanks[0]);
490                delinearize(flattened_routes, routes);
491        }
492        if (level!=level)
493        {
494                for (int i = 0; i < routes.size(); ++i)
495                        for (int k = 0; k < routes[i].size(); ++k)
496                                if (routes[i][k] == cascade.level[0].rank) routes[i].erase(routes[i].begin() + (k--));
497        }
498}
499
500void CParallelTree::updateCirclesForRouting(Coord rootCentre, double rootRadius, int level)
501{
502        if (level + 1 < cascade.num_levels) // children in the cascade have to update first
503        {
504                updateCirclesForRouting(rootCentre, rootRadius, level + 1);
505                rootCentre = treeCascade[level+1].root->centre;
506                rootRadius = treeCascade[level+1].root->radius;
507        }
508
509        // gather circles on this level of the cascade
510        int pg_size;
511        MPI_Comm_size(cascade.level[level].pg_comm, &pg_size);
512        vector<Coord> allRootCentres(pg_size);
513        vector<double> allRootRadia(pg_size);
514        MPI_Allgather(&rootCentre, 3, MPI_DOUBLE, &allRootCentres[0], 3, MPI_DOUBLE, cascade.level[level].pg_comm);
515        MPI_Allgather(&rootRadius, 1, MPI_DOUBLE, &allRootRadia[0],   1, MPI_DOUBLE, cascade.level[level].pg_comm);
516
517        // now allRootsRadia and allRootCentres must be inserted into second levels of us and propagated to root
518        treeCascade[level].root->assignCircleAndPropagateUp(&allRootCentres[0], &allRootRadia[0], assignLevel);
519}
520
521CParallelTree::~CParallelTree()
522{
523        delete [] localElements;
524}
525
526}
Note: See TracBrowser for help on using the repository browser.