source: XIOS/dev/dev_olga/src/transformation/grid_transformation.cpp @ 1599

Last change on this file since 1599 was 1589, checked in by oabramkina, 5 years ago

Backporting r1578 and r1586 to dev, cleaning the code before merging it to XIOS 2.5.

  • Property svn:executable set to *
File size: 31.0 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 <unordered_map>
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    }
349    else
350      return;
351  }
352
353  CContext* context = CContext::getCurrent();
354  CContextClient* client = context->client;
355
356  ListAlgoType::const_iterator itb = listAlgos_.begin(),
357                               ite = listAlgos_.end(), it;
358
359  CGenericAlgorithmTransformation* algo = 0;
360  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
361  for (it = itb; it != ite; ++it)
362  {
363    int elementPositionInGrid = it->first;
364    ETranformationType transType = (it->second).first;
365    int transformationOrder = (it->second).second.first;
366    int algoType = ((it->second).second.second); //algoTypes_[std::distance(itb, it)];
367    SourceDestinationIndexMap globaIndexWeightFromSrcToDst;
368   
369
370    // Create a temporary grid destination which contains transformed element of grid destination and
371    // non-transformed elements to grid source
372    setUpGridDestination(elementPositionInGrid, transType);
373
374    // First of all, select an algorithm
375    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size()))
376    {
377      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoType);
378      algo = algoTransformation_.back();
379    }
380    else
381      algo = algoTransformation_[std::distance(itb, it)];
382
383    if ((0 != algo) &&
384        ((CGenericAlgorithmTransformation::ELEMENT_NO_MODIFICATION_WITH_DATA == algo->type()) ||
385        (CGenericAlgorithmTransformation::ELEMENT_MODIFICATION_WITH_DATA == algo->type()))) // Only registered transformation can be executed
386    {
387      CTimer::get("computeIndexSourceMapping").resume() ;
388      algo->computeIndexSourceMapping(dataAuxInputs);
389      CTimer::get("computeIndexSourceMapping").suspend() ;
390     
391      // ComputeTransformation of global index of each element
392      int elementPosition = it->first;
393      bool nonDistributedActivated = CXios::getin<bool>("activate_non_distributed_transformation",false);
394     
395      if (nonDistributedActivated && !algo->isDistributedTransformation(elementPositionInGrid, gridSource_, tmpGridDestination_) )
396      {
397        vector<int> localSrc ;
398        vector<int> localDst ;
399        vector<double> weight ;
400        int nbLocalIndexOnGridDest;
401        CTimer::get("computeTransformationMappingNonDistributed").resume(); 
402//        algo->computeTransformationMappingNonDistributed(elementPosition, gridSource_, tmpGridDestination_,
403//                                                         localSrc, localDst, weight, localMaskOnGridDest_.back()) ;
404        algo->computeTransformationMappingNonDistributed(elementPosition, gridSource_, tmpGridDestination_, 
405                                                         localSrc, localDst, weight, nbLocalIndexOnGridDest) ;
406        CTimer::get("computeTransformationMappingNonDistributed").suspend(); 
407
408        CTimer::get("computeTransformationMappingConvert").resume(); 
409        nbLocalIndexOnGridDest_.push_back(nbLocalIndexOnGridDest) ;
410        int clientRank=client->clientRank ;
411        {
412          SendingIndexGridSourceMap tmp;
413          localIndexToSendFromGridSource_.push_back(tmp) ;
414          SendingIndexGridSourceMap& src=localIndexToSendFromGridSource_.back() ;
415          CArray<int,1> arrayTmp ;
416          src.insert( pair<int,CArray<int,1> >(clientRank,arrayTmp)) ;
417          CArray<int,1>& array=src[clientRank] ;
418          array.resize(localSrc.size()) ;
419          for(int i=0;i< localSrc.size();++i) array(i)=localSrc[i]  ;   
420        }
421        {
422          RecvIndexGridDestinationMap tmp;
423          localIndexToReceiveOnGridDest_.push_back(tmp) ;
424          RecvIndexGridDestinationMap& dst=localIndexToReceiveOnGridDest_.back() ;
425          vector<pair<int,double> > vectTmp ;
426          dst.insert( pair<int,vector<pair<int,double> >  >(clientRank,vectTmp)) ;
427          vector<pair<int,double> >& vect=dst[clientRank] ;
428          vect.resize(localDst.size()) ;
429          for(int i=0;i< localDst.size();++i) vect[i]=pair<int,double>(localDst[i],weight[i])  ;     
430        }
431        CTimer::get("computeTransformationMappingConvert").suspend(); 
432      }
433      else
434      {
435        CTimer::get("computeGlobalSourceIndex").resume();           
436        algo->computeGlobalSourceIndex(elementPosition,
437                                       gridSource_,
438                                       tmpGridDestination_,
439                                       globaIndexWeightFromSrcToDst);
440                                     
441        CTimer::get("computeGlobalSourceIndex").suspend();           
442        CTimer::get("computeTransformationMapping").resume();     
443        // Compute transformation of global indexes among grids
444        computeTransformationMapping(globaIndexWeightFromSrcToDst);
445        CTimer::get("computeTransformationMapping").suspend(); 
446      } 
447      if (1 < nbNormalAlgos_)
448      {
449        // Now grid destination becomes grid source in a new transformation
450        if (nbAgloTransformation != (nbNormalAlgos_-1)) setUpGridSource(elementPositionInGrid);
451      }
452      ++nbAgloTransformation;
453    }
454  }
455}
456
457/*!
458  Compute exchange index between grid source and grid destination
459  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
460*/
461void CGridTransformation::computeTransformationMapping(const SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
462{
463  CContext* context = CContext::getCurrent();
464  CContextClient* client = context->client;
465  int nbClient = client->clientSize;
466  int clientRank = client->clientRank;
467
468  // Recalculate the distribution of grid destination
469  CDistributionClient distributionClientDest(client->clientRank, tmpGridDestination_);
470  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridDestSendToServer = distributionClientDest.getGlobalLocalDataSendToServer();
471
472  // Update number of local index on each transformation
473  size_t nbLocalIndex = globalLocalIndexGridDestSendToServer.size();
474  nbLocalIndexOnGridDest_.push_back(nbLocalIndex);
475//  localMaskOnGridDest_.push_back(std::vector<bool>());
476//  std::vector<bool>& tmpMask = localMaskOnGridDest_.back();
477//  tmpMask.resize(nbLocalIndex,false);
478
479  // Find out number of index sent from grid source and number of index received on grid destination
480  SourceDestinationIndexMap::const_iterator itbIndex = globaIndexWeightFromSrcToDst.begin(),
481                                            iteIndex = globaIndexWeightFromSrcToDst.end(), itIndex;
482  typedef std::unordered_map<size_t, std::vector<std::pair<size_t,double> > > SendIndexMap;
483  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
484  int connectedClient = globaIndexWeightFromSrcToDst.size();
485  int* recvCount=new int[nbClient];
486  int* displ=new int[nbClient];
487  int* sendRankBuff=new int[connectedClient];
488  int* sendSizeBuff=new int[connectedClient];
489  int n = 0;
490  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
491  {
492    sendRankBuff[n] = itIndex->first;
493    const SendIndexMap& sendIndexMap = itIndex->second;
494    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
495    int sendSize = 0;
496    for (itSend = itbSend; itSend != iteSend; ++itSend)
497    {
498      sendSize += itSend->second.size();
499    }
500    sendSizeBuff[n] = sendSize;
501    sendRankSizeMap[itIndex->first] = sendSize;
502  }
503  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
504
505  displ[0]=0 ;
506  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
507  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
508  int* recvRankBuff=new int[recvSize];
509  int* recvSizeBuff=new int[recvSize];
510  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,client->intraComm);
511  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,client->intraComm);
512  for (int i = 0; i < nbClient; ++i)
513  {
514    int currentPos = displ[i];
515    for (int j = 0; j < recvCount[i]; ++j)
516      if (recvRankBuff[currentPos+j] == clientRank)
517      {
518        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
519      }
520  }
521
522  // Sending global index of grid source to corresponding process as well as the corresponding mask
523  std::vector<MPI_Request> requests;
524  std::vector<MPI_Status> status;
525  std::unordered_map<int, unsigned char* > recvMaskDst;
526  std::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
527  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
528  {
529    int recvRank = itRecv->first;
530    int recvSize = itRecv->second;
531    recvMaskDst[recvRank] = new unsigned char [recvSize];
532    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
533
534    requests.push_back(MPI_Request());
535    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, client->intraComm, &requests.back());
536    requests.push_back(MPI_Request());
537    MPI_Irecv(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 47, client->intraComm, &requests.back());
538  }
539
540  std::unordered_map<int, CArray<size_t,1> > globalIndexDst;
541  std::unordered_map<int, CArray<double,1> > weightDst;
542  std::unordered_map<int, unsigned char* > sendMaskDst;
543  std::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
544  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
545  {
546    int sendRank = itIndex->first;
547    int sendSize = sendRankSizeMap[sendRank];
548    const SendIndexMap& sendIndexMap = itIndex->second;
549    SendIndexMap::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
550    globalIndexDst[sendRank].resize(sendSize);
551    weightDst[sendRank].resize(sendSize);
552    sendMaskDst[sendRank] = new unsigned char [sendSize];
553    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
554    int countIndex = 0;
555    for (itSend = itbSend; itSend != iteSend; ++itSend)
556    {
557      const std::vector<std::pair<size_t,double> >& dstWeight = itSend->second;
558      for (int idx = 0; idx < dstWeight.size(); ++idx)
559      {
560        globalIndexDst[sendRank](countIndex) = dstWeight[idx].first;
561        weightDst[sendRank](countIndex) = dstWeight[idx].second;
562        if (0 < globalLocalIndexGridDestSendToServer.count(dstWeight[idx].first))
563          sendMaskDst[sendRank][countIndex] = 1;
564        else
565          sendMaskDst[sendRank][countIndex] = 0;
566        sendGlobalIndexSrc[sendRank][countIndex] = itSend->first;
567        ++countIndex;
568      }
569    }
570
571    // Send global index source and mask
572    requests.push_back(MPI_Request());
573    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, client->intraComm, &requests.back());
574    requests.push_back(MPI_Request());
575    MPI_Isend(sendMaskDst[sendRank], sendSize, MPI_UNSIGNED_CHAR, sendRank, 47, client->intraComm, &requests.back());
576  }
577
578  status.resize(requests.size());
579  MPI_Waitall(requests.size(), &requests[0], &status[0]);
580
581  // 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
582  std::vector<MPI_Request>().swap(requests);
583  std::vector<MPI_Status>().swap(status);
584  // Okie, on destination side, we will wait for information of masked index of source
585  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
586  {
587    int recvRank = itSend->first;
588    int recvSize = itSend->second;
589
590    requests.push_back(MPI_Request());
591    MPI_Irecv(sendMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
592  }
593
594  // Ok, now we fill in local index of grid source (we even count for masked index)
595  CDistributionClient distributionClientSrc(client->clientRank, gridSource_);
596  CDistributionClient::GlobalLocalDataMap& globalLocalIndexGridSrcSendToServer = distributionClientSrc.getGlobalLocalDataSendToServer();
597  localIndexToSendFromGridSource_.push_back(SendingIndexGridSourceMap());
598  SendingIndexGridSourceMap& tmpSend = localIndexToSendFromGridSource_.back();
599  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
600  {
601    int recvRank = itRecv->first;
602    int recvSize = itRecv->second;
603    unsigned char* recvMask = recvMaskDst[recvRank];
604    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
605    int realSendSize = 0;
606    for (int idx = 0; idx < recvSize; ++idx)
607    {
608      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
609        if (0 < globalLocalIndexGridSrcSendToServer.count(*(recvIndexSrc+idx))) // check whether index source is masked
610         ++realSendSize;
611        else // inform the destination that this index is masked
612         *(recvMask+idx) = 0;
613    }
614
615    tmpSend[recvRank].resize(realSendSize);
616    realSendSize = 0;
617    for (int idx = 0; idx < recvSize; ++idx)
618    {
619      if (0 != (*(recvMask+idx))) // OKie, now we have a demand from non-masked index destination
620      {
621        tmpSend[recvRank](realSendSize) = globalLocalIndexGridSrcSendToServer[*(recvIndexSrc+idx)];
622         ++realSendSize;
623      }
624    }
625
626    // Okie, now inform the destination which source index are masked
627    requests.push_back(MPI_Request());
628    MPI_Isend(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 48, client->intraComm, &requests.back());
629  }
630  status.resize(requests.size());
631  MPI_Waitall(requests.size(), &requests[0], &status[0]);
632
633  // Cool, now we can fill in local index of grid destination (counted for masked index)
634  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
635  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
636  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
637  {
638    int recvRank = itSend->first;
639    int recvSize = itSend->second;
640    unsigned char* recvMask = sendMaskDst[recvRank];
641
642    CArray<size_t,1>& recvIndexDst = globalIndexDst[recvRank];
643    CArray<double,1>& recvWeightDst = weightDst[recvRank];
644    int realRecvSize = 0;
645    for (int idx = 0; idx < recvSize; ++idx)
646    {
647      if (0 != *(recvMask+idx)) // OKie, now we have a non-masked index destination
648         ++realRecvSize;
649    }
650
651    int localIndexDst;
652    recvTmp[recvRank].resize(realRecvSize);
653    realRecvSize = 0;
654    for (int idx = 0; idx < recvSize; ++idx)
655    {
656      if (0 != *(recvMask+idx)) // OKie, now we have a demand from non-masked index destination
657      {
658        recvTmp[recvRank][realRecvSize].first = globalLocalIndexGridDestSendToServer[recvIndexDst(idx)];
659        recvTmp[recvRank][realRecvSize].second = recvWeightDst(idx);
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  std::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  std::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}
Note: See TracBrowser for help on using the repository browser.