source: XIOS/trunk/src/transformation/axis_algorithm_inverse.cpp @ 833

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

Improvements for dht

+) Implement adaptive hierarchy for dht, level of hierarchy depends on number of processes
+) Remove some redundant codes

Test
+) On Curie
+) All tests are correct

File size: 5.7 KB
Line 
1/*!
2   \file axis_algorithm_inverse.hpp
3   \author Ha NGUYEN
4   \since 14 May 2015
5   \date 29 June 2015
6
7   \brief Algorithm for inversing an axis..
8 */
9#include "axis_algorithm_inverse.hpp"
10#include "transformation_mapping.hpp"
11#include "context.hpp"
12#include "context_client.hpp"
13
14namespace xios {
15
16CAxisAlgorithmInverse::CAxisAlgorithmInverse(CAxis* axisDestination, CAxis* axisSource)
17 : CAxisAlgorithmTransformation(axisDestination, axisSource)
18{
19  if (axisDestination->n_glo.getValue() != axisSource->n_glo.getValue())
20  {
21    ERROR("CAxisAlgorithmInverse::CAxisAlgorithmInverse(CAxis* axisDestination, CAxis* axisSource)",
22           << "Two axis have different global size"
23           << "Size of axis source " <<axisSource->getId() << " is " << axisSource->n_glo.getValue()  << std::endl
24           << "Size of axis destionation " <<axisDestination->getId() << " is " << axisDestination->n_glo.getValue());
25  }
26}
27
28void CAxisAlgorithmInverse::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
29{
30  this->transformationMapping_.resize(1);
31  this->transformationWeight_.resize(1);
32
33  TransformationIndexMap& transMap = this->transformationMapping_[0];
34  TransformationWeightMap& transWeight = this->transformationWeight_[0];
35
36  int globalIndexSize = axisDestGlobalIndex_.size();
37  for (int idx = 0; idx < globalIndexSize; ++idx)
38  {
39    transMap[axisDestGlobalIndex_[idx]].push_back(axisDestGlobalSize_-axisDestGlobalIndex_[idx]-1);
40    transWeight[axisDestGlobalIndex_[idx]].push_back(1.0);
41  }
42
43  int niSrc   = axisSrc_->n.getValue();
44  int sizeSrc = axisSrc_->n_glo.getValue();
45  if (niSrc != sizeSrc) updateAxisValue();
46  else
47  {
48    for (int idx = 0; idx < sizeSrc; ++idx)
49    {
50      axisDest_->value(idx) = axisSrc_->value(sizeSrc-idx-1);
51    }
52  }
53}
54
55/*!
56  Update value on axis after inversing
57  After an axis is inversed, not only the data on it must be inversed but also the value
58*/
59void CAxisAlgorithmInverse::updateAxisValue()
60{
61  CContext* context = CContext::getCurrent();
62  CContextClient* client=context->client;
63
64  int niSrc     = axisSrc_->n.getValue();
65  int ibeginSrc = axisSrc_->begin.getValue();
66
67  CTransformationMapping transformationMap(axisDest_, axisSrc_);
68
69  TransformationIndexMap& transMap = this->transformationMapping_[0];
70  TransformationWeightMap& transWeight = this->transformationWeight_[0];
71
72  CTransformationMapping::DestinationIndexMap globaIndexMapFromDestToSource;
73  TransformationIndexMap::const_iterator it = transMap.begin(), ite = transMap.end();
74  int localIndex = 0;
75  for (; it != ite; ++it)
76  {
77    globaIndexMapFromDestToSource[it->first].push_back(make_pair(localIndex,make_pair((it->second)[0], (transWeight[it->first])[0])));
78    ++localIndex;
79  }
80
81  transformationMap.computeTransformationMapping(globaIndexMapFromDestToSource);
82
83  const CTransformationMapping::ReceivedIndexMap& globalIndexToReceive = transformationMap.getGlobalIndexReceivedOnGridDestMapping();
84  const CTransformationMapping::SentIndexMap& globalIndexToSend = transformationMap.getGlobalIndexSendToGridDestMapping();
85
86 // Sending global index of original grid source
87  CTransformationMapping::SentIndexMap::const_iterator itbSend = globalIndexToSend.begin(), itSend,
88                                                       iteSend = globalIndexToSend.end();
89 int sendBuffSize = 0;
90 for (itSend = itbSend; itSend != iteSend; ++itSend) sendBuffSize += (itSend->second).size();
91
92 typedef double Scalar;
93 Scalar* sendBuff, *currentSendBuff;
94 if (0 != sendBuffSize) sendBuff = new Scalar[sendBuffSize];
95 for (StdSize idx = 0; idx < sendBuffSize; ++idx) sendBuff[idx] = NumTraits<Scalar>::sfmax();
96
97 int currentBuffPosition = 0;
98 for (itSend = itbSend; itSend != iteSend; ++itSend)
99 {
100   int destRank = itSend->first;
101   const std::vector<std::pair<int,size_t> >& globalIndexOfCurrentGridSourceToSend = itSend->second;
102   int countSize = globalIndexOfCurrentGridSourceToSend.size();
103   for (int idx = 0; idx < (countSize); ++idx)
104   {
105     int index = globalIndexOfCurrentGridSourceToSend[idx].first;
106     sendBuff[idx+currentBuffPosition] = (axisSrc_->value)(index);
107   }
108   currentSendBuff = sendBuff + currentBuffPosition;
109   MPI_Send(currentSendBuff, countSize, MPI_DOUBLE, destRank, 14, client->intraComm);
110   currentBuffPosition += countSize;
111 }
112
113 // Receiving global index of grid source sending from current grid source
114 CTransformationMapping::ReceivedIndexMap::const_iterator itbRecv = globalIndexToReceive.begin(), itRecv,
115                                                          iteRecv = globalIndexToReceive.end();
116 int recvBuffSize = 0;
117 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv) recvBuffSize += (itRecv->second).size();
118
119 Scalar* recvBuff, *currentRecvBuff;
120 if (0 != recvBuffSize) recvBuff = new Scalar [recvBuffSize];
121 for (StdSize idx = 0; idx < recvBuffSize; ++idx) recvBuff[idx] = NumTraits<Scalar>::sfmax();
122
123 int currentRecvBuffPosition = 0;
124 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
125 {
126   MPI_Status status;
127   int srcRank = itRecv->first;
128   int countSize = (itRecv->second).size();
129   currentRecvBuff = recvBuff + currentRecvBuffPosition;
130   MPI_Recv(currentRecvBuff, countSize, MPI_DOUBLE, srcRank, 14, client->intraComm, &status);
131   currentRecvBuffPosition += countSize;
132 }
133
134 int ibeginDest = axisDest_->begin.getValue();
135 currentRecvBuff = recvBuff;
136 for (itRecv = itbRecv; itRecv != iteRecv; ++itRecv)
137 {
138   int countSize = (itRecv->second).size();
139   for (int idx = 0; idx < countSize; ++idx, ++currentRecvBuff)
140   {
141     int ssize = (itRecv->second)[idx].size();
142     for (int i = 0; i < ssize; ++i)
143     {
144       int index = ((itRecv->second)[idx][i]).first - ibeginDest;
145       (axisDest_->value)(index) = *currentRecvBuff;
146     }
147   }
148 }
149
150 if (0 != sendBuffSize) delete [] sendBuff;
151 if (0 != recvBuffSize) delete [] recvBuff;
152}
153
154}
Note: See TracBrowser for help on using the repository browser.