source: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/grid_transformation.cpp @ 1985

Last change on this file since 1985 was 1985, checked in by ymipsl, 10 months ago

intermediate commit for new tranformation engine?
YM

  • 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, false,
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, false,
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, false,
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  tmpGridDestination_-> solveElementsRefInheritance() ; // ym not sure
250  tempGridDests_.push_back(tmpGridDestination_);
251}
252CATCH
253
254/*!
255  Assign the current grid destination to the grid source in the new transformation.
256The current grid destination plays the role of grid source in next transformation (if any).
257Only element on which the transformation is performed is modified
258  \param [in] elementPositionInGrid position of element in grid
259  \param [in] transType transformation type
260*/
261void CGridTransformation::setUpGridSource(int elementPositionInGrid)
262TRY
263{
264  if (!tempGridSrcs_.empty() && (getNbAlgo()-1) == tempGridSrcs_.size())
265  {
266    tempGridSrcs_.resize(0);
267  }
268
269  std::vector<CScalar*> scalarListDestP = tmpGridDestination_->getScalars();
270  std::vector<CScalar*> scalarListSrcP = gridSource_->getScalars(), scalarSrc;
271
272  std::vector<CAxis*> axisListDestP = tmpGridDestination_->getAxis();
273  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisSrc;
274
275  std::vector<CDomain*> domListDestP = tmpGridDestination_->getDomains();
276  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainSrc;
277
278  CArray<int,1> axisDomainOrderSrc = gridSource_->axis_domain_order;
279  CArray<int,1> axisDomainOrderDst = tmpGridDestination_->axis_domain_order;
280
281  std::map<int,std::pair<int,int> > elementPositionSrc = getElementPosition(gridSource_);
282  std::map<int,std::pair<int,int> > elementPositionDst = getElementPosition(tmpGridDestination_);
283
284  for (int idx = 0; idx < axisDomainOrderDst.numElements(); ++idx)
285  {   
286    if (elementPositionInGrid == idx)
287    {
288      int dimElementDst = elementPositionDst[idx].first;
289      int elementIndex  = elementPositionDst[idx].second;
290      if (2 == dimElementDst)
291      {
292        CDomain* domain = CDomain::createDomain();
293        domain->domain_ref.setValue(domListDestP[elementIndex]->getId());
294        domain->solveRefInheritance(true);
295        domain->checkAttributes();
296        domainSrc.push_back(domain);
297      }
298      else if (1 == dimElementDst)
299      {
300        CAxis* axis = CAxis::createAxis();
301        axis->axis_ref.setValue(axisListDestP[elementIndex]->getId());
302        axis->solveRefInheritance(true);
303        axis->checkAttributes();
304        axisSrc.push_back(axis); 
305      }
306      else
307      {
308        CScalar* scalar = CScalar::createScalar();
309        scalar->scalar_ref.setValue(scalarListDestP[elementIndex]->getId());
310        scalar->solveRefInheritance(true);
311        scalar->checkAttributes();
312        scalarSrc.push_back(scalar);
313      }
314    }
315    else
316    {     
317      int dimElementDst = elementPositionDst[idx].first;
318      int elementIndex  = elementPositionDst[idx].second;
319      switch (dimElementDst)
320      {
321        case 2:
322          domainSrc.push_back(domListDestP[elementIndex]);
323          break;
324        case 1:
325          axisSrc.push_back(axisListDestP[elementIndex]);
326          break;
327        case 0:
328          scalarSrc.push_back(scalarListDestP[elementIndex]);
329          break;
330        default:
331          break;
332      }
333    }
334  }
335
336  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, scalarSrc, tmpGridDestination_->axis_domain_order); 
337  gridSource_->solveElementsRefInheritance() ; // ym not sure
338
339  tempGridSrcs_.push_back(gridSource_);
340}
341CATCH
342
343/*!
344  Perform all transformations
345  For each transformation, there are some things to do:
346  -) Chose the correct algorithm by transformation type and position of element
347  -) Calculate the mapping of global index between the current grid source and grid destination
348  -) Calculate the mapping of global index between current grid DESTINATION and grid SOURCE
349  -) Make current grid destination become grid source in the next transformation
350*/
351void CGridTransformation::computeAll(const std::vector<CArray<double,1>* >& dataAuxInputs, Time timeStamp)
352TRY
353{
354  if (nbNormalAlgos_ < 1) return;
355  if (!auxInputs_.empty() && !dynamicalTransformation_) { dynamicalTransformation_ = true; return; }
356  if (dynamicalTransformation_)
357  {
358    if (timeStamp_.insert(timeStamp).second)   //Reset map
359    {
360      std::list<SendingIndexGridSourceMap>().swap(localIndexToSendFromGridSource_);
361      std::list<RecvIndexGridDestinationMap>().swap(localIndexToReceiveOnGridDest_);
362      std::list<size_t>().swap(nbLocalIndexOnGridDest_);
363    }
364    else
365      return;
366  }
367
368  CContext* context = CContext::getCurrent();
369
370  ListAlgoType::const_iterator itb = listAlgos_.begin(),
371                               ite = listAlgos_.end(), it;
372
373  CGenericAlgorithmTransformation* algo = 0;
374  int nbAgloTransformation = 0; // Only count for executed transformation. Generate domain is a special one, not executed in the list
375  for (it = itb; it != ite; ++it)
376  {
377    int elementPositionInGrid = it->first;
378    ETranformationType transType = (it->second).first;
379    int transformationOrder = (it->second).second.first;
380    int algoType = ((it->second).second.second); //algoTypes_[std::distance(itb, it)];
381    SourceDestinationIndexMap globaIndexWeightFromSrcToDst;
382   
383
384    // Create a temporary grid destination which contains transformed element of grid destination and
385    // non-transformed elements to grid source
386    setUpGridDestination(elementPositionInGrid, transType);
387
388    // First of all, select an algorithm
389    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size()))
390    {
391      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoType);
392      algo = algoTransformation_.back();
393    }
394    else
395      algo = algoTransformation_[std::distance(itb, it)];
396
397    if ((0 != algo) &&
398        ((CGenericAlgorithmTransformation::ELEMENT_NO_MODIFICATION_WITH_DATA == algo->type()) ||
399        (CGenericAlgorithmTransformation::ELEMENT_MODIFICATION_WITH_DATA == algo->type()))) // Only registered transformation can be executed
400    {
401      CTimer::get("computeIndexSourceMapping").resume() ;
402      algo->computeIndexSourceMapping(dataAuxInputs);
403      CTimer::get("computeIndexSourceMapping").suspend() ;
404     
405      // ComputeTransformation of global index of each element
406      int elementPosition = it->first;
407      bool nonDistributedActivated = CXios::getin<bool>("activate_non_distributed_transformation",false);
408     
409      if (nonDistributedActivated && !algo->isDistributedTransformation(elementPositionInGrid, gridSource_, tmpGridDestination_) )
410      {
411        vector<int> localSrc ;
412        vector<int> localDst ;
413        vector<double> weight ;
414        int nbLocalIndexOnGridDest;
415        CTimer::get("computeTransformationMappingNonDistributed").resume(); 
416        algo->computeTransformationMappingNonDistributed(elementPosition, gridSource_, tmpGridDestination_, 
417                                                         localSrc, localDst, weight, nbLocalIndexOnGridDest) ;
418        CTimer::get("computeTransformationMappingNonDistributed").suspend(); 
419
420        CTimer::get("computeTransformationMappingConvert").resume(); 
421        nbLocalIndexOnGridDest_.push_back(nbLocalIndexOnGridDest) ;
422        int clientRank=context->intraCommRank_ ;
423        {
424          SendingIndexGridSourceMap tmp;
425          localIndexToSendFromGridSource_.push_back(tmp) ;
426          SendingIndexGridSourceMap& src=localIndexToSendFromGridSource_.back() ;
427          CArray<int,1> arrayTmp ;
428          src.insert( pair<int,CArray<int,1> >(clientRank,arrayTmp)) ;
429          CArray<int,1>& array=src[clientRank] ;
430          array.resize(localSrc.size()) ;
431          for(int i=0;i< localSrc.size();++i) array(i)=localSrc[i]  ;   
432        }
433        {
434          RecvIndexGridDestinationMap tmp;
435          localIndexToReceiveOnGridDest_.push_back(tmp) ;
436          RecvIndexGridDestinationMap& dst=localIndexToReceiveOnGridDest_.back() ;
437          vector<pair<int,double> > vectTmp ;
438          dst.insert( pair<int,vector<pair<int,double> >  >(clientRank,vectTmp)) ;
439          vector<pair<int,double> >& vect=dst[clientRank] ;
440          vect.resize(localDst.size()) ;
441          for(int i=0;i< localDst.size();++i) vect[i]=pair<int,double>(localDst[i],weight[i])  ;     
442        }
443        CTimer::get("computeTransformationMappingConvert").suspend(); 
444      }
445      else
446      {
447        CTimer::get("computeGlobalSourceIndex").resume();           
448        algo->computeGlobalSourceIndex(elementPosition,
449                                       gridSource_,
450                                       tmpGridDestination_,
451                                       globaIndexWeightFromSrcToDst);
452                                     
453        CTimer::get("computeGlobalSourceIndex").suspend();           
454        CTimer::get("computeTransformationMapping").resume();     
455        // Compute transformation of global indexes among grids
456        computeTransformationMapping(globaIndexWeightFromSrcToDst);
457        CTimer::get("computeTransformationMapping").suspend(); 
458      } 
459      if (1 < nbNormalAlgos_)
460      {
461        // Now grid destination becomes grid source in a new transformation
462        if (nbAgloTransformation != (nbNormalAlgos_-1)) setUpGridSource(elementPositionInGrid);
463      }
464      ++nbAgloTransformation;
465    }
466  }
467}
468CATCH
469
470/*!
471  Compute exchange index between grid source and grid destination
472  \param [in] globalIndexWeightFromDestToSource global index mapping between grid destination and grid source
473*/
474void CGridTransformation::computeTransformationMapping(const SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
475TRY
476{
477  CContext* context = CContext::getCurrent();
478  int nbClient = context->intraCommSize_;
479  int clientRank = context->intraCommRank_;
480
481  // Recalculate the distribution of grid destination
482  CDistributionClient distributionClientDest(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);
488//  localMaskOnGridDest_.push_back(std::vector<bool>());
489//  std::vector<bool>& tmpMask = localMaskOnGridDest_.back();
490//  tmpMask.resize(nbLocalIndex,false);
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, context->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,context->intraComm_);
524  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,context->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, context->intraComm_, &requests.back());
549    requests.push_back(MPI_Request());
550    MPI_Irecv(recvMaskDst[recvRank], recvSize, MPI_UNSIGNED_CHAR, recvRank, 47, context->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))
576          sendMaskDst[sendRank][countIndex] = 1;
577        else
578          sendMaskDst[sendRank][countIndex] = 0;
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, context->intraComm_, &requests.back());
587    requests.push_back(MPI_Request());
588    MPI_Isend(sendMaskDst[sendRank], sendSize, MPI_UNSIGNED_CHAR, sendRank, 47, context->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, context->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(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;
616    unsigned char* recvMask = recvMaskDst[recvRank];
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
625         *(recvMask+idx) = 0;
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, context->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)
647  localIndexToReceiveOnGridDest_.push_back(RecvIndexGridDestinationMap());
648  RecvIndexGridDestinationMap& recvTmp = localIndexToReceiveOnGridDest_.back();
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;
653    unsigned char* recvMask = sendMaskDst[recvRank];
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;
686  for (itChar = sendMaskDst.begin(); itChar != sendMaskDst.end(); ++itChar)
687    delete [] itChar->second;
688  for (itChar = recvMaskDst.begin(); itChar != recvMaskDst.end(); ++itChar)
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*/
714const std::list<CGridTransformation::RecvIndexGridDestinationMap>& CGridTransformation::getLocalIndexToReceiveOnGridDest() const
715TRY
716{
717  return localIndexToReceiveOnGridDest_;
718}
719CATCH
720
721/*!
722 Number of index will be received on the grid destination
723  \return number of index of data
724*/
725const std::list<size_t>& CGridTransformation::getNbLocalIndexToReceiveOnGridDest() const
726TRY
727{
728  return nbLocalIndexOnGridDest_;
729}
730CATCH
731
732}
Note: See TracBrowser for help on using the repository browser.