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

Last change on this file since 1295 was 1295, checked in by yushan, 4 years ago

EP update all

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