source: XIOS/dev/XIOS_DEV_CMIP6/src/transformation/generic_algorithm_transformation.cpp @ 1221

Last change on this file since 1221 was 1221, checked in by mhnguyen, 3 years ago

Fixing a minor bug in transformation: No source index at all for transformation

Test
+) On Curie
+) OK

File size: 27.6 KB
Line 
1/*!
2   \file generic_algorithm_transformation.hpp
3   \author Ha NGUYEN
4   \since 14 May 2015
5   \date 21 Mars 2016
6
7   \brief Interface for all transformation algorithms.
8 */
9#include "generic_algorithm_transformation.hpp"
10#include "context.hpp"
11#include "context_client.hpp"
12#include "client_client_dht_template.hpp"
13#include "utils.hpp"
14#include "timer.hpp"
15
16namespace xios {
17
18CGenericAlgorithmTransformation::CGenericAlgorithmTransformation()
19 : transformationMapping_(), transformationWeight_(), transformationPosition_(),
20   idAuxInputs_(), type_(ELEMENT_NO_MODIFICATION_WITH_DATA), indexElementSrc_(),
21   computedProcSrcNonTransformedElement_(false)
22{
23}
24
25void CGenericAlgorithmTransformation::updateData(CArray<double,1>& dataOut)
26{
27
28}
29
30void CGenericAlgorithmTransformation::apply(const std::vector<std::pair<int,double> >& localIndex,
31                                            const double* dataInput,
32                                            CArray<double,1>& dataOut,
33                                            std::vector<bool>& flagInitial,
34                                            bool ignoreMissingValue)
35{
36  int nbLocalIndex = localIndex.size();   
37  double defaultValue = std::numeric_limits<double>::quiet_NaN();
38  if (ignoreMissingValue)
39  {
40    for (int idx = 0; idx < nbLocalIndex; ++idx)
41    {
42      if (NumTraits<double>::isnan(*(dataInput + idx)))
43      {
44        flagInitial[localIndex[idx].first] = false;
45      }
46      else
47      {
48        dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
49        flagInitial[localIndex[idx].first] = true; // Reset flag to indicate not all data source are nan
50      }
51    }
52
53    // If all data source are nan then data destination must be nan
54    for (int idx = 0; idx < nbLocalIndex; ++idx)
55    {
56      if (!flagInitial[localIndex[idx].first])
57        dataOut(localIndex[idx].first) = defaultValue;
58    }
59  }
60  else
61  {
62    for (int idx = 0; idx < nbLocalIndex; ++idx)
63    {
64      dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
65    }
66  }
67}
68
69void CGenericAlgorithmTransformation::computePositionElements(CGrid* dst, CGrid* src)
70{
71  int idxScalar = 0, idxAxis = 0, idxDomain = 0;
72  CArray<int,1> axisDomainOrderDst = dst->axis_domain_order;
73  for (int i = 0; i < axisDomainOrderDst.numElements(); ++i)
74  {
75    int dimElement = axisDomainOrderDst(i);
76    if (2 == dimElement)
77    {
78      elementPositionInGridDst2DomainPosition_[i] = idxDomain;
79      ++idxDomain;
80    }
81    else if (1 == dimElement)
82    {
83      elementPositionInGridDst2AxisPosition_[i] = idxAxis;
84      ++idxAxis;
85    }
86    else
87    {
88      elementPositionInGridDst2ScalarPosition_[i] = idxScalar;
89      ++idxScalar;
90    }
91  }
92
93  idxScalar = idxAxis = idxDomain = 0;
94  CArray<int,1> axisDomainOrderSrc = src->axis_domain_order;
95  for (int i = 0; i < axisDomainOrderSrc.numElements(); ++i)
96  {
97    int dimElement = axisDomainOrderSrc(i);
98    if (2 == dimElement)
99    {
100      elementPositionInGridSrc2DomainPosition_[i] = idxDomain;
101      ++idxDomain;
102    }
103    else if (1 == dimElement)
104    {
105      elementPositionInGridSrc2AxisPosition_[i] = idxAxis;
106      ++idxAxis;
107    }
108    else
109    {
110      elementPositionInGridSrc2ScalarPosition_[i] = idxScalar;
111      ++idxScalar;
112    }
113  }
114}
115
116/*!
117  This function computes the global indexes of grid source, which the grid destination is in demand.
118  \param[in] elementPositionInGrid position of an element in a grid .E.g: if grid is composed of domain and axis (in order),
119                then position of axis in grid is 1 and domain is positioned at 0.
120  \param[in] gridSrc Grid source
121  \param[in] gridDst Grid destination
122  \param[in\out] globaIndexWeightFromSrcToDst mapping of each global index source and weight to index destination
123*/
124void CGenericAlgorithmTransformation::computeGlobalSourceIndex(int elementPositionInGrid,
125                                                               CGrid* gridSrc,
126                                                               CGrid* gridDst,
127                                                               SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
128 {
129  CContext* context = CContext::getCurrent();
130  CContextClient* client = context->client;
131  int nbClient = client->clientSize;
132
133  typedef boost::unordered_map<int, std::vector<std::pair<int,double> > > SrcToDstMap;
134  int idx;
135
136  // compute position of elements on grids
137  computePositionElements(gridDst, gridSrc);
138  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
139  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
140  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
141  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
142  std::vector<CScalar*> scalarListSrcP  = gridSrc->getScalars();
143  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
144  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
145  CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order; 
146 
147  bool isTransPosEmpty = transformationPosition_.empty();
148  CArray<size_t,1> transPos;
149  if (!isTransPosEmpty) transPos.resize(transformationMapping_.size());
150  std::set<size_t> allIndexSrc; // All index of source, which can be scattered among processes, need for doing transformation
151 
152  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
153  {
154    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
155                                           iteTransMap = transformationMapping_[idxTrans].end();
156    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
157
158    // Build mapping between global source element index and global destination element index.
159    itTransWeight = itbTransWeight;
160    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
161    {
162      const std::vector<int>& srcIndex = itTransMap->second;
163      for (idx = 0; idx < srcIndex.size(); ++idx)
164        allIndexSrc.insert(srcIndex[idx]);
165    }
166
167    if (!isTransPosEmpty)
168    {
169      TransformationPositionMap::const_iterator itPosMap = transformationPosition_[idxTrans].begin();
170      transPos(idxTrans) = itPosMap->second[0];
171    }
172  }
173
174  size_t indexSrcSize = 0;
175  CArray<size_t,1> indexSrc(allIndexSrc.size());
176  for (std::set<size_t>::iterator it = allIndexSrc.begin(); it != allIndexSrc.end(); ++it, ++indexSrcSize)
177    indexSrc(indexSrcSize) = *it;
178
179  // Flag to indicate whether we will recompute distribution of source global index  on processes
180  bool computeGlobalIndexOnProc = false; 
181  if (indexElementSrc_.size() != allIndexSrc.size())
182    computeGlobalIndexOnProc = true; 
183  else
184  {
185    for (std::set<size_t>::iterator it = allIndexSrc.begin(); it != allIndexSrc.end(); ++it)
186      if (0 == indexElementSrc_.count(*it))
187      {
188        computeGlobalIndexOnProc = true;
189        break;       
190      }
191  }
192
193  if (computeGlobalIndexOnProc)
194    indexElementSrc_.swap(allIndexSrc);
195     
196  int sendValue = (computeGlobalIndexOnProc) ? 1 : 0;
197  int recvValue = 0;
198  MPI_Allreduce(&sendValue, &recvValue, 1, MPI_INT, MPI_SUM, client->intraComm);
199  computeGlobalIndexOnProc = (0 < recvValue);
200
201  if (computeGlobalIndexOnProc || !computedProcSrcNonTransformedElement_)
202  {   
203    // Find out global index source of transformed element on corresponding process.   
204    if (globalElementIndexOnProc_.empty())
205      globalElementIndexOnProc_.resize(axisDomainDstOrder.numElements());
206    CClientClientDHTInt::Index2VectorInfoTypeMap globalIndexOfTransformedElementOnProc; 
207    for (idx = 0; idx < axisDomainDstOrder.numElements(); ++idx)
208    {
209      if (idx == elementPositionInGrid)
210        computeExchangeGlobalIndex(indexSrc, axisDomainSrcOrder(idx), globalIndexOfTransformedElementOnProc); //globalElementIndexOnProc[idx]);
211      if (!computedProcSrcNonTransformedElement_)
212      {
213        if (2 == axisDomainDstOrder(idx)) // It's domain
214        {
215          if (idx != elementPositionInGrid)
216            computeExchangeDomainIndex(domainListDestP[elementPositionInGridDst2DomainPosition_[idx]],
217                                       domainListSrcP[elementPositionInGridSrc2DomainPosition_[idx]],
218                                       transPos,
219                                       globalElementIndexOnProc_[idx]);     
220
221        }
222        else if (1 == axisDomainDstOrder(idx))//it's an axis
223        {
224          if (idx != elementPositionInGrid)
225            computeExchangeAxisIndex(axisListDestP[elementPositionInGridDst2AxisPosition_[idx]],
226                                     axisListSrcP[elementPositionInGridSrc2AxisPosition_[idx]],
227                                     transPos,
228                                     globalElementIndexOnProc_[idx]);
229        }
230        else //it's a scalar
231        {
232          if (idx != elementPositionInGrid)
233            computeExchangeScalarIndex(scalarListDestP[elementPositionInGridDst2ScalarPosition_[idx]],
234                                       scalarListSrcP[elementPositionInGridSrc2ScalarPosition_[idx]],
235                                       transPos,
236                                       globalElementIndexOnProc_[idx]);
237
238        }
239      }
240    }
241
242    if (!isTransPosEmpty && !computedProcSrcNonTransformedElement_)
243    {
244      for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx)
245      {
246        if (idx != elementPositionInGrid)
247        {
248          boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,
249                                                                   ite = globalElementIndexOnProc_[idx].end();
250          for (it = itb; it != ite; ++it) it->second.resize(1);
251        }
252      }
253    }
254     
255    if (!computedProcSrcNonTransformedElement_)
256    {
257      for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx)
258      {
259        if (idx != elementPositionInGrid)
260        {
261          boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,
262                                                                   ite = globalElementIndexOnProc_[idx].end();
263          for (it = itb; it != ite; ++it) procOfNonTransformedElements_.insert(it->first);
264          if (procOfNonTransformedElements_.size() == nbClient)
265            break;
266        }
267      }
268    }
269   
270    // Processes contain the source index of transformed element
271    std::set<int> procOfTransformedElement; 
272    CClientClientDHTInt::Index2VectorInfoTypeMap::iterator itIdxb = globalIndexOfTransformedElementOnProc.begin(),
273                                                           itIdxe = globalIndexOfTransformedElementOnProc.end(), itIdx;
274    for (itIdx = itIdxb; itIdx != itIdxe; ++itIdx)
275    {
276      std::vector<int>& tmp = itIdx->second;
277      for (int i = 0; i < tmp.size(); ++i)
278        procOfTransformedElement.insert(tmp[i]);
279      if (tmp.size() == nbClient)
280        break;
281    }                                                           
282   
283    std::set<int>& commonProc = (procOfTransformedElement.size() < procOfNonTransformedElements_.size()) ? procOfTransformedElement
284                              : (!procOfNonTransformedElements_.empty() ? procOfNonTransformedElements_ : procOfTransformedElement);
285   
286    std::vector<int> procContainSrcElementIdx(commonProc.size());
287    int count = 0;
288    for (std::set<int>::iterator it = commonProc.begin(); it != commonProc.end(); ++it)
289      procContainSrcElementIdx[count++] = *it;
290
291    procContainSrcElementIdx_.swap(procContainSrcElementIdx);   
292
293        // For the first time, surely we calculate proc containing non transformed elements
294    if (!computedProcSrcNonTransformedElement_)
295      computedProcSrcNonTransformedElement_ = true;
296  }
297 
298  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans)
299  {
300    TransformationIndexMap::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap,
301                                           iteTransMap = transformationMapping_[idxTrans].end();
302    TransformationWeightMap::const_iterator itbTransWeight = transformationWeight_[idxTrans].begin(), itTransWeight;
303    SrcToDstMap src2DstMap;
304    src2DstMap.rehash(std::ceil(transformationMapping_[idxTrans].size()/src2DstMap.max_load_factor()));
305
306    // Build mapping between global source element index and global destination element index.
307    boost::unordered_map<int,std::vector<size_t> >().swap(globalElementIndexOnProc_[elementPositionInGrid]);
308    std::set<int> tmpCounter;
309    itTransWeight = itbTransWeight;
310    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight)
311    {
312      const std::vector<int>& srcIndex = itTransMap->second;
313      const std::vector<double>& weight = itTransWeight->second;
314      for (idx = 0; idx < srcIndex.size(); ++idx)
315      {
316        src2DstMap[srcIndex[idx]].push_back(make_pair(itTransMap->first, weight[idx]));
317        if (0 == tmpCounter.count(srcIndex[idx]))
318        {         
319          tmpCounter.insert(srcIndex[idx]);
320          for (int j = 0; j < procContainSrcElementIdx_.size(); ++j)
321            globalElementIndexOnProc_[elementPositionInGrid][procContainSrcElementIdx_[j]].push_back(srcIndex[idx]);
322        }
323      }
324    }
325 
326    if (!isTransPosEmpty)
327    {
328      for (idx = 0; idx < globalElementIndexOnProc_.size(); ++idx)
329      {
330        if (idx != elementPositionInGrid)
331        {
332          boost::unordered_map<int,std::vector<size_t> >::iterator itb = globalElementIndexOnProc_[idx].begin(), it,
333                                                                   ite = globalElementIndexOnProc_[idx].end();
334          for (it = itb; it != ite; ++it) it->second[0] = transPos(idxTrans);
335        }
336      }
337    }
338
339    // Ok, now compute global index of grid source and ones of grid destination
340    computeGlobalGridIndexMapping(elementPositionInGrid,
341                                  procContainSrcElementIdx_, //srcRank,
342                                  src2DstMap,
343                                  gridSrc,
344                                  gridDst,
345                                  globalElementIndexOnProc_,
346                                  globaIndexWeightFromSrcToDst);
347  } 
348 }
349
350/*!
351  Compute mapping of global index of grid source and grid destination
352  \param [in] elementPositionInGrid position of element in grid. E.x: grid composed of domain and axis, domain has position 0 and axis 1.
353  \param [in] srcRank rank of client from which we demand global index of element source
354  \param [in] src2DstMap mapping of global index of element source and global index of element destination
355  \param[in] gridSrc Grid source
356  \param[in] gridDst Grid destination
357  \param[in] globalElementIndexOnProc Global index of element source on different client rank
358  \param[out] globaIndexWeightFromSrcToDst Mapping of global index of grid source and grid destination
359*/
360void CGenericAlgorithmTransformation::computeGlobalGridIndexMapping(int elementPositionInGrid,
361                                                                   const std::vector<int>& srcRank,
362                                                                   boost::unordered_map<int, std::vector<std::pair<int,double> > >& src2DstMap,
363                                                                   CGrid* gridSrc,
364                                                                   CGrid* gridDst,
365                                                                   std::vector<boost::unordered_map<int,std::vector<size_t> > >& globalElementIndexOnProc,
366                                                                   SourceDestinationIndexMap& globaIndexWeightFromSrcToDst)
367{
368  std::vector<CDomain*> domainListSrcP = gridSrc->getDomains();
369  std::vector<CAxis*> axisListSrcP = gridSrc->getAxis();
370  std::vector<CScalar*> scalarListSrcP = gridSrc->getScalars();
371  CArray<int,1> axisDomainSrcOrder = gridSrc->axis_domain_order;
372
373  size_t nbElement = axisDomainSrcOrder.numElements();
374  std::vector<size_t> nGlobSrc(nbElement);
375  size_t globalSrcSize = 1;
376  int domainIndex = 0, axisIndex = 0, scalarIndex = 0;
377  for (int idx = 0; idx < nbElement; ++idx)
378  {
379    nGlobSrc[idx] = globalSrcSize;
380    int elementDimension = axisDomainSrcOrder(idx);
381
382    // If this is a domain
383    if (2 == elementDimension)
384    {
385      globalSrcSize *= domainListSrcP[domainIndex]->nj_glo.getValue() * domainListSrcP[domainIndex]->ni_glo.getValue();
386      ++domainIndex;
387    }
388    else if (1 == elementDimension) // So it's an axis
389    {
390      globalSrcSize *= axisListSrcP[axisIndex]->n_glo.getValue();
391      ++axisIndex;
392    }
393    else
394    {
395      globalSrcSize *= 1;
396      ++scalarIndex;
397    }
398  }
399
400  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
401  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
402  std::vector<CScalar*> scalarListDestP = gridDst->getScalars();
403  CArray<int,1> axisDomainDstOrder = gridDst->axis_domain_order;
404
405  std::vector<size_t> nGlobDst(nbElement);
406  size_t globalDstSize = 1;
407  domainIndex = axisIndex = scalarIndex = 0;
408  for (int idx = 0; idx < nbElement; ++idx)
409  {
410    nGlobDst[idx] = globalDstSize;
411    int elementDimension = axisDomainDstOrder(idx);
412
413    // If this is a domain
414    if (2 == elementDimension)
415    {
416      globalDstSize *= domainListDestP[domainIndex]->nj_glo.getValue() * domainListDestP[domainIndex]->ni_glo.getValue();
417      ++domainIndex;
418    }
419    else if (1 == elementDimension) // So it's an axis
420    {
421      globalDstSize *= axisListDestP[axisIndex]->n_glo.getValue();
422      ++axisIndex;
423    }
424    else
425    {
426      globalDstSize *= 1;
427      ++scalarIndex;
428    }
429  }
430
431  for (int i = 0; i < srcRank.size(); ++i)
432  {
433    size_t ssize = 1;
434    int rankSrc = srcRank[i];
435    for (int idx = 0; idx < nbElement; ++idx)
436    {
437      ssize *= (globalElementIndexOnProc[idx][rankSrc]).size();
438    }
439
440    std::vector<int> idxLoop(nbElement,0);
441    std::vector<int> currentIndexSrc(nbElement, 0);
442    std::vector<int> currentIndexDst(nbElement, 0);
443    int innnerLoopSize = (globalElementIndexOnProc[0])[rankSrc].size();
444    size_t idx = 0;
445    while (idx < ssize)
446    {
447      for (int ind = 0; ind < nbElement; ++ind)
448      {
449        if (idxLoop[ind] == (globalElementIndexOnProc[ind])[rankSrc].size())
450        {
451          idxLoop[ind] = 0;
452          ++idxLoop[ind+1];
453        }
454
455        currentIndexDst[ind] = currentIndexSrc[ind] = (globalElementIndexOnProc[ind])[rankSrc][idxLoop[ind]];
456      }
457
458      for (int ind = 0; ind < innnerLoopSize; ++ind)
459      {
460        currentIndexDst[0] = currentIndexSrc[0] = (globalElementIndexOnProc[0])[rankSrc][ind];
461        int globalElementDstIndexSize = 0;
462        if (1 == src2DstMap.count(currentIndexSrc[elementPositionInGrid]))
463        {
464          globalElementDstIndexSize = src2DstMap[currentIndexSrc[elementPositionInGrid]].size();
465        }
466
467        std::vector<size_t> globalDstVecIndex(globalElementDstIndexSize,0);
468        size_t globalSrcIndex = 0;
469        for (int idxElement = 0; idxElement < nbElement; ++idxElement)
470        {
471          if (idxElement == elementPositionInGrid)
472          {
473            for (int k = 0; k < globalElementDstIndexSize; ++k)
474            {
475              globalDstVecIndex[k] += src2DstMap[currentIndexSrc[elementPositionInGrid]][k].first * nGlobDst[idxElement];
476            }
477          }
478          else
479          {
480            for (int k = 0; k < globalElementDstIndexSize; ++k)
481            {
482              globalDstVecIndex[k] += currentIndexDst[idxElement] * nGlobDst[idxElement];
483            }
484          }
485          globalSrcIndex += currentIndexSrc[idxElement] * nGlobSrc[idxElement];
486        }
487
488        for (int k = 0; k < globalElementDstIndexSize; ++k)
489        {
490          globaIndexWeightFromSrcToDst[rankSrc][globalSrcIndex].push_back(make_pair(globalDstVecIndex[k],src2DstMap[currentIndexSrc[elementPositionInGrid]][k].second ));
491        }
492        ++idxLoop[0];
493      }
494      idx += innnerLoopSize;
495    }
496  }
497}
498
499/*!
500  Find out proc and global index of axis source which axis destination is on demande
501  \param[in] scalar Scalar destination
502  \param[in] scalar Scalar source
503  \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.
504  \param[out] globalScalarIndexOnProc Global index of axis source on different procs
505*/
506void CGenericAlgorithmTransformation::computeExchangeScalarIndex(CScalar* scalarDst,
507                                                                 CScalar* scalarSrc,
508                                                                 CArray<size_t,1>& destGlobalIndexPositionInGrid,
509                                                                 boost::unordered_map<int,std::vector<size_t> >& globalScalarIndexOnProc)
510{
511  CContext* context = CContext::getCurrent();
512  CContextClient* client=context->client;
513  int clientRank = client->clientRank;
514  int clientSize = client->clientSize;
515
516  globalScalarIndexOnProc.rehash(std::ceil(clientSize/globalScalarIndexOnProc.max_load_factor()));
517  for (int idx = 0; idx < clientSize; ++idx)
518  {
519    globalScalarIndexOnProc[idx].push_back(0);
520  }
521}
522
523/*!
524  Find out proc and global index of axis source which axis destination is on demande
525  \param[in] axisDst Axis destination
526  \param[in] axisSrc Axis source
527  \param[in] destGlobalIndexPositionInGrid Relative position of axis corresponds to other element of grid.
528  \param[out] globalAxisIndexOnProc Global index of axis source on different procs
529*/
530void CGenericAlgorithmTransformation::computeExchangeAxisIndex(CAxis* axisDst,
531                                                               CAxis* axisSrc,
532                                                               CArray<size_t,1>& destGlobalIndexPositionInGrid,
533                                                               boost::unordered_map<int,std::vector<size_t> >& globalAxisIndexOnProc)
534{
535  CContext* context = CContext::getCurrent();
536  CContextClient* client=context->client;
537  int clientRank = client->clientRank;
538  int clientSize = client->clientSize;
539
540  size_t globalIndex;
541  int nIndexSize = axisSrc->index.numElements();
542  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
543  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
544  for (int idx = 0; idx < nIndexSize; ++idx)
545  {
546    globalIndex = axisSrc->index(idx);
547    globalIndex2ProcRank[globalIndex].push_back(clientRank);
548  }
549
550  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
551  CArray<size_t,1> globalAxisIndex(axisDst->index.numElements());
552  for (int idx = 0; idx < globalAxisIndex.numElements(); ++idx)
553  {
554    globalAxisIndex(idx) = axisDst->index(idx);
555  }
556  dhtIndexProcRank.computeIndexInfoMapping(globalAxisIndex);
557
558  std::vector<int> countIndex(clientSize,0);
559  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
560  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
561                                                               ite = computedGlobalIndexOnProc.end();
562  for (it = itb; it != ite; ++it)
563  {
564    const std::vector<int>& procList = it->second;
565    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
566  }
567
568  globalAxisIndexOnProc.rehash(std::ceil(clientSize/globalAxisIndexOnProc.max_load_factor()));
569  for (int idx = 0; idx < clientSize; ++idx)
570  {
571    if (0 != countIndex[idx])
572    {
573      globalAxisIndexOnProc[idx].resize(countIndex[idx]);
574      countIndex[idx] = 0;
575    }
576  }
577
578  for (it = itb; it != ite; ++it)
579  {
580    const std::vector<int>& procList = it->second;
581    for (int idx = 0; idx < procList.size(); ++idx)
582    {
583      globalAxisIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
584      ++countIndex[procList[idx]];
585    }
586  }
587}
588
589/*!
590  Find out proc and global index of domain source which domain destination is on demande
591  \param[in] domainDst Domain destination
592  \param[in] domainSrc Domain source
593  \param[in] destGlobalIndexPositionInGrid Relative position of domain corresponds to other element of grid.
594  \param[out] globalDomainIndexOnProc Global index of domain source on different procs
595*/
596void CGenericAlgorithmTransformation::computeExchangeDomainIndex(CDomain* domainDst,
597                                                                 CDomain* domainSrc,
598                                                                 CArray<size_t,1>& destGlobalIndexPositionInGrid,
599                                                                 boost::unordered_map<int,std::vector<size_t> >& globalDomainIndexOnProc)
600{
601  CContext* context = CContext::getCurrent();
602  CContextClient* client=context->client;
603  int clientRank = client->clientRank;
604  int clientSize = client->clientSize;
605
606  int niGlobSrc = domainSrc->ni_glo.getValue();
607  size_t globalIndex;
608  int i_ind, j_ind;
609  int nIndexSize = (destGlobalIndexPositionInGrid.isEmpty()) ? domainSrc->i_index.numElements()
610                                                             : destGlobalIndexPositionInGrid.numElements();
611  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
612  globalIndex2ProcRank.rehash(std::ceil(nIndexSize/globalIndex2ProcRank.max_load_factor()));
613  if (destGlobalIndexPositionInGrid.isEmpty())
614  {
615    for (int idx = 0; idx < nIndexSize; ++idx)
616    {
617      i_ind=domainSrc->i_index(idx) ;
618      j_ind=domainSrc->j_index(idx) ;
619
620      globalIndex = i_ind + j_ind * niGlobSrc;
621      globalIndex2ProcRank[globalIndex].resize(1);
622      globalIndex2ProcRank[globalIndex][0] = clientRank;
623    }
624  }
625  else
626  {
627    for (int idx = 0; idx < nIndexSize; ++idx)
628    {
629      globalIndex2ProcRank[destGlobalIndexPositionInGrid(idx)].push_back(clientRank);
630    }
631  }
632
633  CArray<size_t,1> globalDomainIndex;
634  if (destGlobalIndexPositionInGrid.isEmpty())
635  {
636    int niGlobDst = domainDst->ni_glo.getValue();
637    globalDomainIndex.resize(domainDst->i_index.numElements());
638    nIndexSize = domainDst->i_index.numElements();
639
640    for (int idx = 0; idx < nIndexSize; ++idx)
641    {
642      i_ind=domainDst->i_index(idx) ;
643      j_ind=domainDst->j_index(idx) ;
644
645      globalDomainIndex(idx) = i_ind + j_ind * niGlobDst;
646    }
647  }
648  else
649  {
650    globalDomainIndex.reference(destGlobalIndexPositionInGrid);
651  }
652
653  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
654  dhtIndexProcRank.computeIndexInfoMapping(globalDomainIndex);
655
656  std::vector<int> countIndex(clientSize,0);
657  const CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
658  CClientClientDHTInt::Index2VectorInfoTypeMap::const_iterator itb = computedGlobalIndexOnProc.begin(), it,
659                                                               ite = computedGlobalIndexOnProc.end();
660  for (it = itb; it != ite; ++it)
661  {
662    const std::vector<int>& procList = it->second;
663    for (int idx = 0; idx < procList.size(); ++idx) ++countIndex[procList[idx]];
664  }
665
666  globalDomainIndexOnProc.rehash(std::ceil(clientSize/globalDomainIndexOnProc.max_load_factor()));
667  for (int idx = 0; idx < clientSize; ++idx)
668  {
669    if (0 != countIndex[idx])
670    {
671      globalDomainIndexOnProc[idx].resize(countIndex[idx]);
672      countIndex[idx] = 0;
673    }
674  }
675
676  for (it = itb; it != ite; ++it)
677  {
678    const std::vector<int>& procList = it->second;
679    for (int idx = 0; idx < procList.size(); ++idx)
680    {
681      globalDomainIndexOnProc[procList[idx]][countIndex[procList[idx]]] = it->first;
682      ++countIndex[procList[idx]];
683    }
684  }
685}
686
687/*!
688  Compute index mapping between element source and element destination with an auxiliary inputs which determine
689position of each mapped index in global index of grid destination.
690  \param [in] dataAuxInputs auxiliary inputs
691*/
692void CGenericAlgorithmTransformation::computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs)
693{
694  computeIndexSourceMapping_(dataAuxInputs);
695}
696
697std::vector<StdString> CGenericAlgorithmTransformation::getIdAuxInputs()
698{
699  return idAuxInputs_;
700}
701
702CGenericAlgorithmTransformation::AlgoTransType CGenericAlgorithmTransformation::type()
703{
704  return type_;
705}
706
707}
Note: See TracBrowser for help on using the repository browser.