source: XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp @ 1328

Last change on this file since 1328 was 1328, checked in by yushan, 6 years ago

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