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

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

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