source: XIOS/dev/XIOS_DEV_CMIP6/src/transformation/grid_transformation.cpp @ 1399

Last change on this file since 1399 was 1399, checked in by ymipsl, 4 years ago

2 update for spatial tranformations :

  • One bug fix introduce in r1389 concerning domain interpolation or other distributed spatial distransfomation
  • Specific optimisation for transformation on element wich is not distributed across processes

This optimisation is experimental and will be activated only when using a variable in xios context :

<variable id="activate_non_distributed_transformation" type="bool">true</variable>

Later when more succesfull test, the variable use will removed...

YM

  • Property svn:executable set to *
File size: 31.1 KB
Line 
1/*!
2   \file grid_transformation.cpp
3   \author Ha NGUYEN
4   \since 14 May 2015
5   \date 02 Jul 2015
6
7   \brief Interface for all transformations.
8 */
9#include "grid_transformation.hpp"
10#include "grid_transformation_factory_impl.hpp"
11#include "algo_types.hpp"
12#include "context.hpp"
13#include "context_client.hpp"
14#include "distribution_client.hpp"
15#include "mpi_tag.hpp"
16#include "grid.hpp"
17#include <boost/unordered_map.hpp>
18#include "timer.hpp"
19
20namespace xios {
21CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source)
22 : CGridTransformationSelector(destination, source),
23  tmpGridDestination_(destination), originalGridSource_(source),
24  tempGridSrcs_(), tempGridDests_(),
25  dynamicalTransformation_(false), timeStamp_()
26{
27}
28
29CGridTransformation::~CGridTransformation()
30{
31}
32
33/*!
34  Select algorithm of a scalar corresponding to its transformation type and its position in each element
35  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 0 and position of axis is 1
36  \param [in] transType transformation type, for now we have
37  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
38*/
39void CGridTransformation::selectScalarAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
40{
41  std::vector<CScalar*> scaListDestP = gridDestination_->getScalars();
42  int scalarDstIndex =  elementPositionInGridDst2ScalarPosition_[elementPositionInGrid];
43  CScalar::TransMapTypes trans = scaListDestP[scalarDstIndex]->getAllTransformations();
44  CScalar::TransMapTypes::const_iterator it = trans.begin();
45
46  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
47  CGenericAlgorithmTransformation* algo = 0;
48  algo = CGridTransformationFactory<CScalar>::createTransformation(transType,
49                                                                  gridDestination_,
50                                                                  gridSource_,
51                                                                  it->second,
52                                                                  elementPositionInGrid,
53                                                                  elementPositionInGridSrc2ScalarPosition_,
54                                                                  elementPositionInGridSrc2AxisPosition_,
55                                                                  elementPositionInGridSrc2DomainPosition_,
56                                                                  elementPositionInGridDst2ScalarPosition_,
57                                                                  elementPositionInGridDst2AxisPosition_,
58                                                                  elementPositionInGridDst2DomainPosition_);
59  algoTransformation_.push_back(algo);
60}
61
62/*!
63  Select algorithm of an axis corresponding to its transformation type and its position in each element
64  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 0 and position of axis is 1
65  \param [in] transType transformation type, for now we have zoom_axis, inverse_axis, interpolate_axis
66  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
67*/
68void CGridTransformation::selectAxisAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
69{
70  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
71  int axisDstIndex =  elementPositionInGridDst2AxisPosition_[elementPositionInGrid];
72  CAxis::TransMapTypes trans = axisListDestP[axisDstIndex]->getAllTransformations();
73  CAxis::TransMapTypes::const_iterator it = trans.begin();
74  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
75
76  CGenericAlgorithmTransformation* algo = 0;
77  algo = CGridTransformationFactory<CAxis>::createTransformation(transType,
78                                                                 gridDestination_,
79                                                                 gridSource_,
80                                                                 it->second,
81                                                                 elementPositionInGrid,
82                                                                 elementPositionInGridSrc2ScalarPosition_,
83                                                                 elementPositionInGridSrc2AxisPosition_,
84                                                                 elementPositionInGridSrc2DomainPosition_,
85                                                                 elementPositionInGridDst2ScalarPosition_,
86                                                                 elementPositionInGridDst2AxisPosition_,
87                                                                 elementPositionInGridDst2DomainPosition_);
88  algoTransformation_.push_back(algo);
89}
90
91/*!
92  Select algorithm of a domain corresponding to its transformation type and its position in each element
93  \param [in] elementPositionInGrid position of element in grid. e.g: a grid has 1 domain and 1 axis, then position of domain is 0 and position of axis is 1
94  \param [in] transType transformation type, for now we have zoom_domain, interpolate_domain
95  \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
96*/
97void CGridTransformation::selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder)
98{
99  std::vector<CDomain*> domainListDestP = gridDestination_->getDomains();
100  int domainIndex =  elementPositionInGridDst2DomainPosition_[elementPositionInGrid];
101  CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
102  CDomain::TransMapTypes::const_iterator it = trans.begin();
103  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation 
104
105  CGenericAlgorithmTransformation* algo = 0;
106  algo = CGridTransformationFactory<CDomain>::createTransformation(transType,
107                                                                   gridDestination_,
108                                                                   gridSource_,
109                                                                   it->second,
110                                                                   elementPositionInGrid,
111                                                                   elementPositionInGridSrc2ScalarPosition_,
112                                                                   elementPositionInGridSrc2AxisPosition_,
113                                                                   elementPositionInGridSrc2DomainPosition_,
114                                                                   elementPositionInGridDst2ScalarPosition_,
115                                                                   elementPositionInGridDst2AxisPosition_,
116                                                                   elementPositionInGridDst2DomainPosition_);
117  algoTransformation_.push_back(algo);
118}
119
120/*!
121  Find position of element in a grid as well as its type (domain, axis, scalar) and position in its own element list
122  \return element position: map<int,<int,int> > corresponds to <element position in grid, <element type, element position in element list> >
123*/
124std::map<int,std::pair<int,int> > CGridTransformation::getElementPosition(CGrid* grid)
125{
126  std::vector<CScalar*> scalarListP = grid->getScalars(); 
127  std::vector<CAxis*> axisListP = grid->getAxis();
128  std::vector<CDomain*> domListP = grid->getDomains(); 
129  CArray<int,1> axisDomainOrder = grid->axis_domain_order;
130  int scalarIndex = 0, axisIndex = 0, domainIndex = 0;
131  int nbElement = axisDomainOrder.numElements(), elementDim;
132  std::map<int,std::pair<int,int> > elementPosition;
133  for (int idx = 0; idx < nbElement; ++idx)
134  {
135    elementDim = axisDomainOrder(idx);
136    switch (elementDim)
137    {
138      case 2:
139        elementPosition[idx] = std::make_pair(elementDim, domainIndex);
140        ++domainIndex;
141        break;
142      case 1:
143        elementPosition[idx] = std::make_pair(elementDim, axisIndex);
144        ++axisIndex;
145        break;
146      case 0:
147        elementPosition[idx] = std::make_pair(elementDim, scalarIndex);
148        ++scalarIndex;
149        break;
150      default:
151        break;       
152    }
153  }
154
155  return elementPosition; 
156}
157
158/*!
159  If there are more than one transformation, a new temporary grid will be created and it will play the role of grid destination.
160This new created one keeps a pointer to the real transformed element of grid destination and generate new copies of other elements from grid source.
161  \param [in] elementPositionInGrid position of element in grid
162  \param [in] transType transformation type
163*/
164void CGridTransformation::setUpGridDestination(int elementPositionInGrid, ETranformationType transType)
165{
166  if (isSpecialTransformation(transType)) return;
167
168  if (!tempGridDests_.empty() && (getNbAlgo() == tempGridDests_.size()))
169  {
170    tempGridDests_.resize(0);
171  }
172
173  if (1 == getNbAlgo()) 
174  {
175    tmpGridDestination_ = gridDestination_;
176    return;
177  }
178
179  std::vector<CScalar*> scalarListDestP = gridDestination_->getScalars();
180  std::vector<CScalar*> scalarListSrcP = gridSource_->getScalars(), scalarDst;
181
182  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis();
183  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisDst;
184
185  std::vector<CDomain*> domListDestP = gridDestination_->getDomains();
186  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainDst;
187
188  CArray<int,1> axisDomainOrderSrc = gridSource_->axis_domain_order;
189  CArray<int,1> axisDomainOrderDst = gridDestination_->axis_domain_order;
190
191  std::map<int,std::pair<int,int> > elementPositionSrc = getElementPosition(gridSource_);
192  std::map<int,std::pair<int,int> > elementPositionDst = getElementPosition(gridDestination_);
193
194  CArray<int,1> elementOrder(axisDomainOrderDst.numElements());
195  for (int idx = 0; idx < axisDomainOrderDst.numElements(); ++idx)
196  {
197    if (elementPositionInGrid == idx)
198    {
199      int dimElementDst = elementPositionDst[idx].first;
200      int elementIndex  = elementPositionDst[idx].second;
201      switch (dimElementDst) 
202      {
203        case 2:
204          domainDst.push_back(domListDestP[elementIndex]);
205          break;
206        case 1:
207          axisDst.push_back(axisListDestP[elementIndex]);
208          break;
209        case 0:
210          scalarDst.push_back(scalarListDestP[elementIndex]);
211          break;
212        default:
213          break;
214      }
215      elementOrder(idx) = dimElementDst;
216    }
217    else
218    {     
219      int dimElementSrc = elementPositionSrc[idx].first;
220      int elementIndex  = elementPositionSrc[idx].second;
221      switch (dimElementSrc)
222      {
223        case 2:
224          domainDst.push_back(domListSrcP[elementIndex]);
225          break;
226        case 1:
227          axisDst.push_back(axisListSrcP[elementIndex]);
228          break;
229        case 0:
230          scalarDst.push_back(scalarListSrcP[elementIndex]);
231          break;
232        default:
233          break;
234      }
235      elementOrder(idx) = dimElementSrc;
236    }
237  }
238
239  tmpGridDestination_ = CGrid::createGrid(domainDst, axisDst, scalarDst, elementOrder); 
240  tempGridDests_.push_back(tmpGridDestination_);
241}
242
243/*!
244  Assign the current grid destination to the grid source in the new transformation.
245The current grid destination plays the role of grid source in next transformation (if any).
246Only element on which the transformation is performed is modified
247  \param [in] elementPositionInGrid position of element in grid
248  \param [in] transType transformation type
249*/
250void CGridTransformation::setUpGridSource(int elementPositionInGrid)
251{
252  if (!tempGridSrcs_.empty() && (getNbAlgo()-1) == tempGridSrcs_.size())
253  {
254    tempGridSrcs_.resize(0);
255  }
256
257  std::vector<CScalar*> scalarListDestP = tmpGridDestination_->getScalars();
258  std::vector<CScalar*> scalarListSrcP = gridSource_->getScalars(), scalarSrc;
259
260  std::vector<CAxis*> axisListDestP = tmpGridDestination_->getAxis();
261  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisSrc;
262
263  std::vector<CDomain*> domListDestP = tmpGridDestination_->getDomains();
264  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainSrc;
265
266  CArray<int,1> axisDomainOrderSrc = gridSource_->axis_domain_order;
267  CArray<int,1> axisDomainOrderDst = tmpGridDestination_->axis_domain_order;
268
269  std::map<int,std::pair<int,int> > elementPositionSrc = getElementPosition(gridSource_);
270  std::map<int,std::pair<int,int> > elementPositionDst = getElementPosition(tmpGridDestination_);
271
272  for (int idx = 0; idx < axisDomainOrderDst.numElements(); ++idx)
273  {   
274    if (elementPositionInGrid == idx)
275    {
276      int dimElementDst = elementPositionDst[idx].first;
277      int elementIndex  = elementPositionDst[idx].second;
278      if (2 == dimElementDst)
279      {
280        CDomain* domain = CDomain::createDomain();
281        domain->domain_ref.setValue(domListDestP[elementIndex]->getId());
282        domain->solveRefInheritance(true);
283        domain->checkAttributesOnClient();
284        domainSrc.push_back(domain);
285      }
286      else if (1 == dimElementDst)
287      {
288        CAxis* axis = CAxis::createAxis();
289        axis->axis_ref.setValue(axisListDestP[elementIndex]->getId());
290        axis->solveRefInheritance(true);
291        axis->checkAttributesOnClient();
292        axisSrc.push_back(axis); 
293      }
294      else
295      {
296        CScalar* scalar = CScalar::createScalar();
297        scalar->scalar_ref.setValue(scalarListDestP[elementIndex]->getId());
298        scalar->solveRefInheritance(true);
299        scalar->checkAttributesOnClient();
300        scalarSrc.push_back(scalar);
301      }
302    }
303    else
304    {     
305      int dimElementDst = elementPositionDst[idx].first;
306      int elementIndex  = elementPositionDst[idx].second;
307      switch (dimElementDst)
308      {
309        case 2:
310          domainSrc.push_back(domListDestP[elementIndex]);
311          break;
312        case 1:
313          axisSrc.push_back(axisListDestP[elementIndex]);
314          break;
315        case 0:
316          scalarSrc.push_back(scalarListDestP[elementIndex]);
317          break;
318        default:
319          break;
320      }
321    }
322  }
323
324  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, scalarSrc, tmpGridDestination_->axis_domain_order); 
325
326  tempGridSrcs_.push_back(gridSource_);
327}
328
329/*!
330  Perform all transformations
331  For each transformation, there are some things to do:
332  -) Chose the correct algorithm by transformation type and position of element
333  -) Calculate the mapping of global index between the current grid source and grid destination
334  -) Calculate the mapping of global index between current grid DESTINATION and grid SOURCE
335  -) Make current grid destination become grid source in the next transformation
336*/
337void CGridTransformation::computeAll(const std::vector<CArray<double,1>* >& dataAuxInputs, Time timeStamp)
338{
339  if (nbNormalAlgos_ < 1) return;
340  if (!auxInputs_.empty() && !dynamicalTransformation_) { dynamicalTransformation_ = true; return; }
341  if (dynamicalTransformation_)
342  {
343    if (timeStamp_.insert(timeStamp).second)   //Reset map
344    {
345      std::list<SendingIndexGridSourceMap>().swap(localIndexToSendFromGridSource_);
346      std::list<RecvIndexGridDestinationMap>().swap(localIndexToReceiveOnGridDest_);
347      std::list<size_t>().swap(nbLocalIndexOnGridDest_);
348      std::list<std::vector<bool> >().swap(localMaskOnGridDest_);
349    }
350    else
351      return;
352  }
353
354  CContext* context = CContext::getCurrent();
355  CContextClient* client = context->client;
356
357  ListAlgoType::const_iterator itb = listAlgos_.begin(),
358                               ite = listAlgos_.end(), it;
359
360  CGenericAlgorithmTransformation* algo = 0;
361  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
362  for (it = itb; it != ite; ++it)
363  {
364    int elementPositionInGrid = it->first;
365    ETranformationType transType = (it->second).first;
366    int transformationOrder = (it->second).second.first;
367    int algoType = ((it->second).second.second); //algoTypes_[std::distance(itb, it)];
368    SourceDestinationIndexMap globaIndexWeightFromSrcToDst;
369   
370
371    // Create a temporary grid destination which contains transformed element of grid destination and
372    // non-transformed elements to grid source
373    setUpGridDestination(elementPositionInGrid, transType);
374
375    // First of all, select an algorithm
376    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size()))
377    {
378      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoType);
379      algo = algoTransformation_.back();
380    }
381    else
382      algo = algoTransformation_[std::distance(itb, it)];
383
384    if ((0 != algo) &&
385        ((CGenericAlgorithmTransformation::ELEMENT_NO_MODIFICATION_WITH_DATA == algo->type()) ||
386        (CGenericAlgorithmTransformation::ELEMENT_MODIFICATION_WITH_DATA == algo->type()))) // Only registered transformation can be executed
387    {
388      algo->computeIndexSourceMapping(dataAuxInputs);
389
390      // ComputeTransformation of global index of each element
391      int elementPosition = it->first;
392      bool nonDistributedActivated = CXios::getin<bool>("activate_non_distributed_transformation",false);
393     
394      if (nonDistributedActivated && !algo->isDistributedTransformation(elementPositionInGrid, gridSource_, tmpGridDestination_) )
395      {
396        vector<int> localSrc ;
397        vector<int> localDst ;
398        vector<double> weight ;
399        localMaskOnGridDest_.push_back(vector<bool>()) ;
400        CTimer::get("computeTransformationMappingNonDistributed").resume(); 
401        algo->computeTransformationMappingNonDistributed(elementPosition, gridSource_, tmpGridDestination_, 
402                                                         localSrc, localDst, weight, localMaskOnGridDest_.back()) ;
403        CTimer::get("computeTransformationMappingNonDistributed").suspend(); 
404
405        CTimer::get("computeTransformationMappingConvert").resume(); 
406        nbLocalIndexOnGridDest_.push_back(localMaskOnGridDest_.back().size()) ;
407        int clientRank=client->clientRank ;
408        {
409          SendingIndexGridSourceMap tmp;
410          localIndexToSendFromGridSource_.push_back(tmp) ;
411          SendingIndexGridSourceMap& src=localIndexToSendFromGridSource_.back() ;
412          CArray<int,1> arrayTmp ;
413          src.insert( pair<int,CArray<int,1> >(clientRank,arrayTmp)) ;
414          CArray<int,1>& array=src[clientRank] ;
415          array.resize(localSrc.size()) ;
416          for(int i=0;i< localSrc.size();++i) array(i)=localSrc[i]  ;   
417        }
418        {
419          RecvIndexGridDestinationMap tmp;
420          localIndexToReceiveOnGridDest_.push_back(tmp) ;
421          RecvIndexGridDestinationMap& dst=localIndexToReceiveOnGridDest_.back() ;
422          vector<pair<int,double> > vectTmp ;
423          dst.insert( pair<int,vector<pair<int,double> >  >(clientRank,vectTmp)) ;
424          vector<pair<int,double> >& vect=dst[clientRank] ;
425          vect.resize(localDst.size()) ;
426          for(int i=0;i< localDst.size();++i) vect[i]=pair<int,double>(localDst[i],weight[i])  ;     
427        }
428        CTimer::get("computeTransformationMappingConvert").suspend(); 
429      }
430      else
431      {
432        CTimer::get("computeGlobalSourceIndex").resume();           
433        algo->computeGlobalSourceIndex(elementPosition,
434                                       gridSource_,
435                                       tmpGridDestination_,
436                                       globaIndexWeightFromSrcToDst);
437                                     
438        CTimer::get("computeGlobalSourceIndex").suspend();           
439        CTimer::get("computeTransformationMapping").resume();     
440        // Compute transformation of global indexes among grids
441        computeTransformationMapping(globaIndexWeightFromSrcToDst);
442        CTimer::get("computeTransformationMapping").suspend(); 
443      } 
444      if (1 < nbNormalAlgos_)
445      {
446        // Now grid destination becomes grid source in a new transformation
447        if (nbAgloTransformation != (nbNormalAlgos_-1)) setUpGridSource(elementPositionInGrid);
448      }
449      ++nbAgloTransformation;
450    }
451  }
452}
453
454/*!
455  Compute exchange index between grid source and grid destination
456  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
457*/
458void CGridTransformation::computeTransformationMapping(const SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
459{
460  CContext* context = CContext::getCurrent();
461  CContextClient* client = context->client;
462  int nbClient = client->clientSize;
463  int clientRank = client->clientRank;
464
465  // Recalculate the distribution of grid destination
466  CDistributionClient distributionClientDest(client->clientRank, tmpGridDestination_);
467  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
468
469  // Update number of local index on each transformation
470  size_t nbLocalIndex = globalLocalIndexGridDestSendToServer.size();
471  nbLocalIndexOnGridDest_.push_back(nbLocalIndex);
472  localMaskOnGridDest_.push_back(std::vector<bool>());
473  std::vector<bool>& tmpMask = localMaskOnGridDest_.back();
474  tmpMask.resize(nbLocalIndex,false);
475
476  // Find out number of index sent from grid source and number of index received on grid destination
477  SourceDestinationIndexMap::const_iterator itbIndex = globaIndexWeightFromSrcToDst.begin(),
478                                            iteIndex = globaIndexWeightFromSrcToDst.end(), itIndex;
479  typedef boost::unordered_map<size_t, std::vector<std::pair<size_t,double> > > SendIndexMap;
480  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
481  int connectedClient = globaIndexWeightFromSrcToDst.size();
482  int* recvCount=new int[nbClient];
483  int* displ=new int[nbClient];
484  int* sendRankBuff=new int[connectedClient];
485  int* sendSizeBuff=new int[connectedClient];
486  int n = 0;
487  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
488  {
489    sendRankBuff[n] = itIndex->first;
490    const SendIndexMap& sendIndexMap = itIndex->second;
491    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
492    int sendSize = 0;
493    for (itSend = itbSend; itSend != iteSend; ++itSend)
494    {
495      sendSize += itSend->second.size();
496    }
497    sendSizeBuff[n] = sendSize;
498    sendRankSizeMap[itIndex->first] = sendSize;
499  }
500  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
501
502  displ[0]=0 ;
503  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
504  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
505  int* recvRankBuff=new int[recvSize];
506  int* recvSizeBuff=new int[recvSize];
507  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,client->intraComm);
508  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,client->intraComm);
509  for (int i = 0; i < nbClient; ++i)
510  {
511    int currentPos = displ[i];
512    for (int j = 0; j < recvCount[i]; ++j)
513      if (recvRankBuff[currentPos+j] == clientRank)
514      {
515        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
516      }
517  }
518
519  // Sending global index of grid source to corresponding process as well as the corresponding mask
520  std::vector<MPI_Request> requests;
521  std::vector<MPI_Status> status;
522  boost::unordered_map<int, unsigned char* > recvMaskDst;
523  boost::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
524  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
525  {
526    int recvRank = itRecv->first;
527    int recvSize = itRecv->second;
528    recvMaskDst[recvRank] = new unsigned char [recvSize];
529    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
530
531    requests.push_back(MPI_Request());
532    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, client->intraComm, &requests.back());
533    requests.push_back(MPI_Request());
534    MPI_Irecv(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 47, client->intraComm, &requests.back());
535  }
536
537  boost::unordered_map<int, CArray<size_t,1> > globalIndexDst;
538  boost::unordered_map<int, CArray<double,1> > weightDst;
539  boost::unordered_map<int, unsigned char* > sendMaskDst;
540  boost::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
541  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
542  {
543    int sendRank = itIndex->first;
544    int sendSize = sendRankSizeMap[sendRank];
545    const SendIndexMap& sendIndexMap = itIndex->second;
546    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
547    globalIndexDst[sendRank].resize(sendSize);
548    weightDst[sendRank].resize(sendSize);
549    sendMaskDst[sendRank] = new unsigned char [sendSize];
550    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
551    int countIndex = 0;
552    for (itSend = itbSend; itSend != iteSend; ++itSend)
553    {
554      const std::vector<std::pair<size_t,double> >& dstWeight = itSend->second;
555      for (int idx = 0; idx < dstWeight.size(); ++idx)
556      {
557        globalIndexDst[sendRank](countIndex) = dstWeight[idx].first;
558        weightDst[sendRank](countIndex) = dstWeight[idx].second;
559        if (0 < globalLocalIndexGridDestSendToServer.count(dstWeight[idx].first))
560          sendMaskDst[sendRank][countIndex] = 1;
561        else
562          sendMaskDst[sendRank][countIndex] = 0;
563        sendGlobalIndexSrc[sendRank][countIndex] = itSend->first;
564        ++countIndex;
565      }
566    }
567
568    // Send global index source and mask
569    requests.push_back(MPI_Request());
570    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, client->intraComm, &requests.back());
571    requests.push_back(MPI_Request());
572    MPI_Isend(sendMaskDst[sendRank], sendSize, MPI_UNSIGNED_CHAR, sendRank, 47, client->intraComm, &requests.back());
573  }
574
575  status.resize(requests.size());
576  MPI_Waitall(requests.size(), &requests[0], &status[0]);
577
578  // Okie, now use the mask to identify which index source we need to send, then also signal the destination which masked index we will return
579  std::vector<MPI_Request>().swap(requests);
580  std::vector<MPI_Status>().swap(status);
581  // Okie, on destination side, we will wait for information of masked index of source
582  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
583  {
584    int recvRank = itSend->first;
585    int recvSize = itSend->second;
586
587    requests.push_back(MPI_Request());
588    MPI_Irecv(sendMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
589  }
590
591  // Ok, now we fill in local index of grid source (we even count for masked index)
592  CDistributionClient distributionClientSrc(client->clientRank, gridSource_);
593  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridSrcSendToServer = distributionClientSrc.getGlobalLocalDataSendToServer();
594  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
595  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
596  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
597  {
598    int recvRank = itRecv->first;
599    int recvSize = itRecv->second;
600    unsigned char* recvMask = recvMaskDst[recvRank];
601    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
602    int realSendSize = 0;
603    for (int idx = 0; idx < recvSize; ++idx)
604    {
605      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
606        if (0 < globalLocalIndexGridSrcSendToServer.count(*(recvIndexSrc+idx))) // check whether index source is masked
607         ++realSendSize;
608        else // inform the destination that this index is masked
609         *(recvMask+idx) = 0;
610    }
611
612    tmpSend[recvRank].resize(realSendSize);
613    realSendSize = 0;
614    for (int idx = 0; idx < recvSize; ++idx)
615    {
616      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
617      {
618        tmpSend[recvRank](realSendSize) = globalLocalIndexGridSrcSendToServer[*(recvIndexSrc+idx)];
619         ++realSendSize;
620      }
621    }
622
623    // Okie, now inform the destination which source index are masked
624    requests.push_back(MPI_Request());
625    MPI_Isend(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
626  }
627  status.resize(requests.size());
628  MPI_Waitall(requests.size(), &requests[0], &status[0]);
629
630  // Cool, now we can fill in local index of grid destination (counted for masked index)
631  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
632  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
633  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
634  {
635    int recvRank = itSend->first;
636    int recvSize = itSend->second;
637    unsigned char* recvMask = sendMaskDst[recvRank];
638
639    CArray<size_t,1>& recvIndexDst = globalIndexDst[recvRank];
640    CArray<double,1>& recvWeightDst = weightDst[recvRank];
641    int realRecvSize = 0;
642    for (int idx = 0; idx < recvSize; ++idx)
643    {
644      if (0 != *(recvMask+idx)) // OKie, now we have a non-masked index destination
645         ++realRecvSize;
646    }
647
648    int localIndexDst;
649    recvTmp[recvRank].resize(realRecvSize);
650    realRecvSize = 0;
651    for (int idx = 0; idx < recvSize; ++idx)
652    {
653      if (0 != *(recvMask+idx)) // OKie, now we have a demand from non-masked index destination
654      {
655        recvTmp[recvRank][realRecvSize].first = globalLocalIndexGridDestSendToServer[recvIndexDst(idx)];
656        recvTmp[recvRank][realRecvSize].second = recvWeightDst(idx);
657        tmpMask[globalLocalIndexGridDestSendToServer[recvIndexDst(idx)]] = true;
658         ++realRecvSize;
659      }
660    }
661  }
662
663  delete [] recvCount;
664  delete [] displ;
665  delete [] sendRankBuff;
666  delete [] recvRankBuff;
667  delete [] sendSizeBuff;
668  delete [] recvSizeBuff;
669
670  boost::unordered_map<int, unsigned char* >::const_iterator itChar;
671  for (itChar = sendMaskDst.begin(); itChar != sendMaskDst.end(); ++itChar)
672    delete [] itChar->second;
673  for (itChar = recvMaskDst.begin(); itChar != recvMaskDst.end(); ++itChar)
674    delete [] itChar->second;
675  boost::unordered_map<int, unsigned long* >::const_iterator itLong;
676  for (itLong = sendGlobalIndexSrc.begin(); itLong != sendGlobalIndexSrc.end(); ++itLong)
677    delete [] itLong->second;
678  for (itLong = recvGlobalIndexSrc.begin(); itLong != recvGlobalIndexSrc.end(); ++itLong)
679    delete [] itLong->second;
680
681}
682
683/*!
684  Local index of data which need sending from the grid source
685  \return local index of data
686*/
687const std::list<CGridTransformation::SendingIndexGridSourceMap>& CGridTransformation::getLocalIndexToSendFromGridSource() const
688{
689  return localIndexToSendFromGridSource_;
690}
691
692/*!
693  Local index of data which will be received on the grid destination
694  \return local index of data
695*/
696const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
697{
698  return localIndexToReceiveOnGridDest_;
699}
700
701/*!
702 Number of index will be received on the grid destination
703  \return number of index of data
704*/
705const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
706{
707  return nbLocalIndexOnGridDest_;
708}
709
710/*!
711  Local mask of data which will be received on the grid destination
712  \return local mask of data
713*/
714const std::list<std::vector<bool> >& CGridTransformation::getLocalMaskIndexOnGridDest() const
715{
716  return localMaskOnGridDest_;
717}
718
719}
Note: See TracBrowser for help on using the repository browser.