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

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

More timer instrumentation...

YM

  • Property svn:executable set to *
File size: 31.2 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      CTimer::get("computeIndexSourceMapping").resume() ;
389      algo->computeIndexSourceMapping(dataAuxInputs);
390      CTimer::get("computeIndexSourceMapping").suspend() ;
391     
392      // ComputeTransformation of global index of each element
393      int elementPosition = it->first;
394      bool nonDistributedActivated = CXios::getin<bool>("activate_non_distributed_transformation",false);
395     
396      if (nonDistributedActivated && !algo->isDistributedTransformation(elementPositionInGrid, gridSource_, tmpGridDestination_) )
397      {
398        vector<int> localSrc ;
399        vector<int> localDst ;
400        vector<double> weight ;
401        localMaskOnGridDest_.push_back(vector<bool>()) ;
402        CTimer::get("computeTransformationMappingNonDistributed").resume(); 
403        algo->computeTransformationMappingNonDistributed(elementPosition, gridSource_, tmpGridDestination_, 
404                                                         localSrc, localDst, weight, localMaskOnGridDest_.back()) ;
405        CTimer::get("computeTransformationMappingNonDistributed").suspend(); 
406
407        CTimer::get("computeTransformationMappingConvert").resume(); 
408        nbLocalIndexOnGridDest_.push_back(localMaskOnGridDest_.back().size()) ;
409        int clientRank=client->clientRank ;
410        {
411          SendingIndexGridSourceMap tmp;
412          localIndexToSendFromGridSource_.push_back(tmp) ;
413          SendingIndexGridSourceMap& src=localIndexToSendFromGridSource_.back() ;
414          CArray<int,1> arrayTmp ;
415          src.insert( pair<int,CArray<int,1> >(clientRank,arrayTmp)) ;
416          CArray<int,1>& array=src[clientRank] ;
417          array.resize(localSrc.size()) ;
418          for(int i=0;i< localSrc.size();++i) array(i)=localSrc[i]  ;   
419        }
420        {
421          RecvIndexGridDestinationMap tmp;
422          localIndexToReceiveOnGridDest_.push_back(tmp) ;
423          RecvIndexGridDestinationMap& dst=localIndexToReceiveOnGridDest_.back() ;
424          vector<pair<int,double> > vectTmp ;
425          dst.insert( pair<int,vector<pair<int,double> >  >(clientRank,vectTmp)) ;
426          vector<pair<int,double> >& vect=dst[clientRank] ;
427          vect.resize(localDst.size()) ;
428          for(int i=0;i< localDst.size();++i) vect[i]=pair<int,double>(localDst[i],weight[i])  ;     
429        }
430        CTimer::get("computeTransformationMappingConvert").suspend(); 
431      }
432      else
433      {
434        CTimer::get("computeGlobalSourceIndex").resume();           
435        algo->computeGlobalSourceIndex(elementPosition,
436                                       gridSource_,
437                                       tmpGridDestination_,
438                                       globaIndexWeightFromSrcToDst);
439                                     
440        CTimer::get("computeGlobalSourceIndex").suspend();           
441        CTimer::get("computeTransformationMapping").resume();     
442        // Compute transformation of global indexes among grids
443        computeTransformationMapping(globaIndexWeightFromSrcToDst);
444        CTimer::get("computeTransformationMapping").suspend(); 
445      } 
446      if (1 < nbNormalAlgos_)
447      {
448        // Now grid destination becomes grid source in a new transformation
449        if (nbAgloTransformation != (nbNormalAlgos_-1)) setUpGridSource(elementPositionInGrid);
450      }
451      ++nbAgloTransformation;
452    }
453  }
454}
455
456/*!
457  Compute exchange index between grid source and grid destination
458  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
459*/
460void CGridTransformation::computeTransformationMapping(const SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
461{
462  CContext* context = CContext::getCurrent();
463  CContextClient* client = context->client;
464  int nbClient = client->clientSize;
465  int clientRank = client->clientRank;
466
467  // Recalculate the distribution of grid destination
468  CDistributionClient distributionClientDest(client->clientRank, tmpGridDestination_);
469  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
470
471  // Update number of local index on each transformation
472  size_t nbLocalIndex = globalLocalIndexGridDestSendToServer.size();
473  nbLocalIndexOnGridDest_.push_back(nbLocalIndex);
474  localMaskOnGridDest_.push_back(std::vector<bool>());
475  std::vector<bool>& tmpMask = localMaskOnGridDest_.back();
476  tmpMask.resize(nbLocalIndex,false);
477
478  // Find out number of index sent from grid source and number of index received on grid destination
479  SourceDestinationIndexMap::const_iterator itbIndex = globaIndexWeightFromSrcToDst.begin(),
480                                            iteIndex = globaIndexWeightFromSrcToDst.end(), itIndex;
481  typedef boost::unordered_map<size_t, std::vector<std::pair<size_t,double> > > SendIndexMap;
482  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
483  int connectedClient = globaIndexWeightFromSrcToDst.size();
484  int* recvCount=new int[nbClient];
485  int* displ=new int[nbClient];
486  int* sendRankBuff=new int[connectedClient];
487  int* sendSizeBuff=new int[connectedClient];
488  int n = 0;
489  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
490  {
491    sendRankBuff[n] = itIndex->first;
492    const SendIndexMap& sendIndexMap = itIndex->second;
493    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
494    int sendSize = 0;
495    for (itSend = itbSend; itSend != iteSend; ++itSend)
496    {
497      sendSize += itSend->second.size();
498    }
499    sendSizeBuff[n] = sendSize;
500    sendRankSizeMap[itIndex->first] = sendSize;
501  }
502  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
503
504  displ[0]=0 ;
505  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
506  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
507  int* recvRankBuff=new int[recvSize];
508  int* recvSizeBuff=new int[recvSize];
509  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,client->intraComm);
510  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,client->intraComm);
511  for (int i = 0; i < nbClient; ++i)
512  {
513    int currentPos = displ[i];
514    for (int j = 0; j < recvCount[i]; ++j)
515      if (recvRankBuff[currentPos+j] == clientRank)
516      {
517        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
518      }
519  }
520
521  // Sending global index of grid source to corresponding process as well as the corresponding mask
522  std::vector<MPI_Request> requests;
523  std::vector<MPI_Status> status;
524  boost::unordered_map<int, unsigned char* > recvMaskDst;
525  boost::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
526  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
527  {
528    int recvRank = itRecv->first;
529    int recvSize = itRecv->second;
530    recvMaskDst[recvRank] = new unsigned char [recvSize];
531    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
532
533    requests.push_back(MPI_Request());
534    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, client->intraComm, &requests.back());
535    requests.push_back(MPI_Request());
536    MPI_Irecv(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 47, client->intraComm, &requests.back());
537  }
538
539  boost::unordered_map<int, CArray<size_t,1> > globalIndexDst;
540  boost::unordered_map<int, CArray<double,1> > weightDst;
541  boost::unordered_map<int, unsigned char* > sendMaskDst;
542  boost::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
543  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
544  {
545    int sendRank = itIndex->first;
546    int sendSize = sendRankSizeMap[sendRank];
547    const SendIndexMap& sendIndexMap = itIndex->second;
548    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
549    globalIndexDst[sendRank].resize(sendSize);
550    weightDst[sendRank].resize(sendSize);
551    sendMaskDst[sendRank] = new unsigned char [sendSize];
552    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
553    int countIndex = 0;
554    for (itSend = itbSend; itSend != iteSend; ++itSend)
555    {
556      const std::vector<std::pair<size_t,double> >& dstWeight = itSend->second;
557      for (int idx = 0; idx < dstWeight.size(); ++idx)
558      {
559        globalIndexDst[sendRank](countIndex) = dstWeight[idx].first;
560        weightDst[sendRank](countIndex) = dstWeight[idx].second;
561        if (0 < globalLocalIndexGridDestSendToServer.count(dstWeight[idx].first))
562          sendMaskDst[sendRank][countIndex] = 1;
563        else
564          sendMaskDst[sendRank][countIndex] = 0;
565        sendGlobalIndexSrc[sendRank][countIndex] = itSend->first;
566        ++countIndex;
567      }
568    }
569
570    // Send global index source and mask
571    requests.push_back(MPI_Request());
572    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, client->intraComm, &requests.back());
573    requests.push_back(MPI_Request());
574    MPI_Isend(sendMaskDst[sendRank], sendSize, MPI_UNSIGNED_CHAR, sendRank, 47, client->intraComm, &requests.back());
575  }
576
577  status.resize(requests.size());
578  MPI_Waitall(requests.size(), &requests[0], &status[0]);
579
580  // 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
581  std::vector<MPI_Request>().swap(requests);
582  std::vector<MPI_Status>().swap(status);
583  // Okie, on destination side, we will wait for information of masked index of source
584  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
585  {
586    int recvRank = itSend->first;
587    int recvSize = itSend->second;
588
589    requests.push_back(MPI_Request());
590    MPI_Irecv(sendMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
591  }
592
593  // Ok, now we fill in local index of grid source (we even count for masked index)
594  CDistributionClient distributionClientSrc(client->clientRank, gridSource_);
595  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridSrcSendToServer = distributionClientSrc.getGlobalLocalDataSendToServer();
596  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
597  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
598  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
599  {
600    int recvRank = itRecv->first;
601    int recvSize = itRecv->second;
602    unsigned char* recvMask = recvMaskDst[recvRank];
603    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
604    int realSendSize = 0;
605    for (int idx = 0; idx < recvSize; ++idx)
606    {
607      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
608        if (0 < globalLocalIndexGridSrcSendToServer.count(*(recvIndexSrc+idx))) // check whether index source is masked
609         ++realSendSize;
610        else // inform the destination that this index is masked
611         *(recvMask+idx) = 0;
612    }
613
614    tmpSend[recvRank].resize(realSendSize);
615    realSendSize = 0;
616    for (int idx = 0; idx < recvSize; ++idx)
617    {
618      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
619      {
620        tmpSend[recvRank](realSendSize) = globalLocalIndexGridSrcSendToServer[*(recvIndexSrc+idx)];
621         ++realSendSize;
622      }
623    }
624
625    // Okie, now inform the destination which source index are masked
626    requests.push_back(MPI_Request());
627    MPI_Isend(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
628  }
629  status.resize(requests.size());
630  MPI_Waitall(requests.size(), &requests[0], &status[0]);
631
632  // Cool, now we can fill in local index of grid destination (counted for masked index)
633  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
634  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
635  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
636  {
637    int recvRank = itSend->first;
638    int recvSize = itSend->second;
639    unsigned char* recvMask = sendMaskDst[recvRank];
640
641    CArray<size_t,1>& recvIndexDst = globalIndexDst[recvRank];
642    CArray<double,1>& recvWeightDst = weightDst[recvRank];
643    int realRecvSize = 0;
644    for (int idx = 0; idx < recvSize; ++idx)
645    {
646      if (0 != *(recvMask+idx)) // OKie, now we have a non-masked index destination
647         ++realRecvSize;
648    }
649
650    int localIndexDst;
651    recvTmp[recvRank].resize(realRecvSize);
652    realRecvSize = 0;
653    for (int idx = 0; idx < recvSize; ++idx)
654    {
655      if (0 != *(recvMask+idx)) // OKie, now we have a demand from non-masked index destination
656      {
657        recvTmp[recvRank][realRecvSize].first = globalLocalIndexGridDestSendToServer[recvIndexDst(idx)];
658        recvTmp[recvRank][realRecvSize].second = recvWeightDst(idx);
659        tmpMask[globalLocalIndexGridDestSendToServer[recvIndexDst(idx)]] = true;
660         ++realRecvSize;
661      }
662    }
663  }
664
665  delete [] recvCount;
666  delete [] displ;
667  delete [] sendRankBuff;
668  delete [] recvRankBuff;
669  delete [] sendSizeBuff;
670  delete [] recvSizeBuff;
671
672  boost::unordered_map<int, unsigned char* >::const_iterator itChar;
673  for (itChar = sendMaskDst.begin(); itChar != sendMaskDst.end(); ++itChar)
674    delete [] itChar->second;
675  for (itChar = recvMaskDst.begin(); itChar != recvMaskDst.end(); ++itChar)
676    delete [] itChar->second;
677  boost::unordered_map<int, unsigned long* >::const_iterator itLong;
678  for (itLong = sendGlobalIndexSrc.begin(); itLong != sendGlobalIndexSrc.end(); ++itLong)
679    delete [] itLong->second;
680  for (itLong = recvGlobalIndexSrc.begin(); itLong != recvGlobalIndexSrc.end(); ++itLong)
681    delete [] itLong->second;
682
683}
684
685/*!
686  Local index of data which need sending from the grid source
687  \return local index of data
688*/
689const std::list<CGridTransformation::SendingIndexGridSourceMap>& CGridTransformation::getLocalIndexToSendFromGridSource() const
690{
691  return localIndexToSendFromGridSource_;
692}
693
694/*!
695  Local index of data which will be received on the grid destination
696  \return local index of data
697*/
698const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
699{
700  return localIndexToReceiveOnGridDest_;
701}
702
703/*!
704 Number of index will be received on the grid destination
705  \return number of index of data
706*/
707const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
708{
709  return nbLocalIndexOnGridDest_;
710}
711
712/*!
713  Local mask of data which will be received on the grid destination
714  \return local mask of data
715*/
716const std::list<std::vector<bool> >& CGridTransformation::getLocalMaskIndexOnGridDest() const
717{
718  return localMaskOnGridDest_;
719}
720
721}
Note: See TracBrowser for help on using the repository browser.