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

Last change on this file since 1784 was 1784, checked in by ymipsl, 2 years ago
  • Preparing coupling functionalities.
  • Make some cleaner things

YM

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