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

Last change on this file since 1620 was 1612, checked in by oabramkina, 22 months ago

Dev: adding exception handling.

To activate it, compilation flag -DXIOS_EXCEPTION should be added.

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