source: XIOS/dev/branch_yushan_merged/src/transformation/grid_transformation.cpp @ 1203

Last change on this file since 1203 was 1203, checked in by yushan, 5 years ago

prep to merge with trunk @1200

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