source:XIOS/trunk/src/transformation/grid_transformation.cpp@1639

Last change on this file since 1639 was 1639, checked in by yushan, 3 years ago

revert erroneous commit on trunk

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