source: XIOS/trunk/src/transformation/generic_algorithm_transformation.cpp @ 979

Last change on this file since 979 was 979, checked in by mhnguyen, 5 years ago

Implementing AVERAGE operation

+) Add average operation for reduction transformation

Test
+) On Curie
+) Correct

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