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

Last change on this file since 1542 was 1542, checked in by oabramkina, 3 years ago

Replacing Boost's unordered_map and shared_pointer by its STL counterparts.

Two notes for Curie:

  • one can see the content of unordered_map with ddt only if XIOS has been compiled with gnu
  • XIOS will not compile any more with pgi (all available versions use old STL which are not up to the c++11 norms)
File size: 11.4 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 "context.hpp"
11#include "context_client.hpp"
12#include "axis.hpp"
13#include "grid.hpp"
14#include "grid_transformation_factory_impl.hpp"
15#include "inverse_axis.hpp"
16#include "client_client_dht_template.hpp"
17
18namespace xios {
19
20CGenericAlgorithmTransformation* CAxisAlgorithmInverse::create(CGrid* gridDst, CGrid* gridSrc,
21                                                               CTransformation<CAxis>* transformation,
22                                                               int elementPositionInGrid,
23                                                               std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
24                                                               std::map<int, int>& elementPositionInGridSrc2AxisPosition,
25                                                               std::map<int, int>& elementPositionInGridSrc2DomainPosition,
26                                                               std::map<int, int>& elementPositionInGridDst2ScalarPosition,
27                                                               std::map<int, int>& elementPositionInGridDst2AxisPosition,
28                                                               std::map<int, int>& elementPositionInGridDst2DomainPosition)
29{
30  std::vector<CAxis*> axisListDestP = gridDst->getAxis();
31  std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
32
33  CInverseAxis* inverseAxis = dynamic_cast<CInverseAxis*> (transformation);
34  int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
35  int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
36
37  return (new CAxisAlgorithmInverse(axisListDestP[axisDstIndex], axisListSrcP[axisSrcIndex], inverseAxis));
38}
39
40bool CAxisAlgorithmInverse::registerTrans()
41{
42  CGridTransformationFactory<CAxis>::registerTransformation(TRANS_INVERSE_AXIS, create);
43}
44
45
46CAxisAlgorithmInverse::CAxisAlgorithmInverse(CAxis* axisDestination, CAxis* axisSource, CInverseAxis* inverseAxis)
47 : CAxisAlgorithmTransformation(axisDestination, axisSource)
48{
49  if (axisDestination->n_glo.getValue() != axisSource->n_glo.getValue())
50  {
51    ERROR("CAxisAlgorithmInverse::CAxisAlgorithmInverse(CAxis* axisDestination, CAxis* axisSource)",
52           << "Two axis have different global size"
53           << "Size of axis source " <<axisSource->getId() << " is " << axisSource->n_glo.getValue()  << std::endl
54           << "Size of axis destination " <<axisDestination->getId() << " is " << axisDestination->n_glo.getValue());
55  }
56}
57
58void CAxisAlgorithmInverse::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
59{
60  this->transformationMapping_.resize(1);
61  this->transformationWeight_.resize(1);
62
63  TransformationIndexMap& transMap = this->transformationMapping_[0];
64  TransformationWeightMap& transWeight = this->transformationWeight_[0];
65
66  int globalIndexSize = axisDestGlobalIndex_.size();
67  for (int idx = 0; idx < globalIndexSize; ++idx)
68  {
69    transMap[axisDestGlobalIndex_[idx]].push_back(axisDestGlobalSize_-axisDestGlobalIndex_[idx]-1);
70    transWeight[axisDestGlobalIndex_[idx]].push_back(1.0);
71  }
72
73  int niSrc   = axisSrc_->n.getValue();
74  int sizeSrc = axisSrc_->n_glo.getValue();
75  if (niSrc != sizeSrc) updateAxisValue();
76  else
77  {
78    for (int idx = 0; idx < sizeSrc; ++idx)
79    {
80      axisDest_->value(idx) = axisSrc_->value(sizeSrc-idx-1);
81    }
82  }
83}
84
85/*!
86  Update value on axis after inversing
87  After an axis is inversed, not only the data on it must be inversed but also the value
88*/
89void CAxisAlgorithmInverse::updateAxisValue()
90{
91  CContext* context = CContext::getCurrent();
92  CContextClient* client=context->client;
93  int clientRank = client->clientRank;
94  int nbClient = client->clientSize;
95
96  int niSrc     = axisSrc_->n.getValue();
97  int ibeginSrc = axisSrc_->begin.getValue();
98  int nSrc = axisSrc_->index.numElements();
99
100  CClientClientDHTInt::Index2VectorInfoTypeMap globalIndex2ProcRank;
101  for (int idx = 0; idx < nSrc; ++idx)
102  {
103    if ((axisSrc_->mask)(idx))
104    {
105      globalIndex2ProcRank[(axisSrc_->index)(idx)].resize(1);
106      globalIndex2ProcRank[(axisSrc_->index)(idx)][0] = clientRank;
107    }
108  }
109
110  typedef std::unordered_map<size_t, std::vector<double> > GlobalIndexMapFromSrcToDest;
111  GlobalIndexMapFromSrcToDest globalIndexMapFromSrcToDest;
112  TransformationIndexMap& transMap = this->transformationMapping_[0];
113  TransformationIndexMap::const_iterator itb = transMap.begin(), ite = transMap.end(), it;
114  CArray<size_t,1> globalSrcIndex(transMap.size());
115  int localIndex = 0;
116  for (it = itb; it != ite; ++it)
117  {
118    size_t srcIndex = it->second[0];
119    globalIndexMapFromSrcToDest[srcIndex].resize(1);
120    globalIndexMapFromSrcToDest[srcIndex][0] = it->first;
121    globalSrcIndex(localIndex) = srcIndex;
122    ++localIndex;
123  }
124
125  CClientClientDHTInt dhtIndexProcRank(globalIndex2ProcRank, client->intraComm);
126  dhtIndexProcRank.computeIndexInfoMapping(globalSrcIndex);
127  CClientClientDHTInt::Index2VectorInfoTypeMap& computedGlobalIndexOnProc = dhtIndexProcRank.getInfoIndexMap();
128  std::unordered_map<int, std::vector<size_t> > globalSrcIndexSendToProc;
129  for (int idx = 0; idx < localIndex; ++idx)
130  {
131    size_t tmpIndex = globalSrcIndex(idx);
132    if (1 == computedGlobalIndexOnProc.count(tmpIndex))
133    {
134      std::vector<int>& tmpVec = computedGlobalIndexOnProc[tmpIndex];
135      globalSrcIndexSendToProc[tmpVec[0]].push_back(tmpIndex);
136    }
137  }
138
139  std::unordered_map<int, std::vector<size_t> >::const_iterator itbIndex = globalSrcIndexSendToProc.begin(), itIndex,
140                                                                  iteIndex = globalSrcIndexSendToProc.end();
141  std::map<int,int> sendRankSizeMap,recvRankSizeMap;
142  int connectedClient = globalSrcIndexSendToProc.size();
143  int* recvCount=new int[nbClient];
144  int* displ=new int[nbClient];
145  int* sendRankBuff=new int[connectedClient];
146  int* sendSizeBuff=new int[connectedClient];
147  int n = 0;
148  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex, ++n)
149  {
150    sendRankBuff[n] = itIndex->first;
151    int sendSize = itIndex->second.size();
152    sendSizeBuff[n] = sendSize;
153    sendRankSizeMap[itIndex->first] = sendSize;
154  }
155  MPI_Allgather(&connectedClient,1,MPI_INT,recvCount,1,MPI_INT,client->intraComm);
156
157  displ[0]=0 ;
158  for(int n=1;n<nbClient;n++) displ[n]=displ[n-1]+recvCount[n-1];
159  int recvSize=displ[nbClient-1]+recvCount[nbClient-1];
160  int* recvRankBuff=new int[recvSize];
161  int* recvSizeBuff=new int[recvSize];
162  MPI_Allgatherv(sendRankBuff,connectedClient,MPI_INT,recvRankBuff,recvCount,displ,MPI_INT,client->intraComm);
163  MPI_Allgatherv(sendSizeBuff,connectedClient,MPI_INT,recvSizeBuff,recvCount,displ,MPI_INT,client->intraComm);
164  for (int i = 0; i < nbClient; ++i)
165  {
166    int currentPos = displ[i];
167    for (int j = 0; j < recvCount[i]; ++j)
168      if (recvRankBuff[currentPos+j] == clientRank)
169      {
170        recvRankSizeMap[i] = recvSizeBuff[currentPos+j];
171      }
172  }
173
174  // Sending global index of grid source to corresponding process as well as the corresponding mask
175  std::vector<MPI_Request> requests;
176  std::vector<MPI_Status> status;
177  std::unordered_map<int, unsigned long* > recvGlobalIndexSrc;
178  std::unordered_map<int, double* > sendValueToDest;
179  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
180  {
181    int recvRank = itRecv->first;
182    int recvSize = itRecv->second;
183    recvGlobalIndexSrc[recvRank] = new unsigned long [recvSize];
184    sendValueToDest[recvRank] = new double [recvSize];
185
186    requests.push_back(MPI_Request());
187    MPI_Irecv(recvGlobalIndexSrc[recvRank], recvSize, MPI_UNSIGNED_LONG, recvRank, 46, client->intraComm, &requests.back());
188  }
189
190  std::unordered_map<int, unsigned long* > sendGlobalIndexSrc;
191  std::unordered_map<int, double* > recvValueFromSrc;
192  for (itIndex = itbIndex; itIndex != iteIndex; ++itIndex)
193  {
194    int sendRank = itIndex->first;
195    int sendSize = sendRankSizeMap[sendRank];
196    const std::vector<size_t>& sendIndexMap = itIndex->second;
197    std::vector<size_t>::const_iterator itbSend = sendIndexMap.begin(), iteSend = sendIndexMap.end(), itSend;
198    sendGlobalIndexSrc[sendRank] = new unsigned long [sendSize];
199    recvValueFromSrc[sendRank] = new double [sendSize];
200    int countIndex = 0;
201    for (itSend = itbSend; itSend != iteSend; ++itSend)
202    {
203      sendGlobalIndexSrc[sendRank][countIndex] = *itSend;
204      ++countIndex;
205    }
206
207    // Send global index source and mask
208    requests.push_back(MPI_Request());
209    MPI_Isend(sendGlobalIndexSrc[sendRank], sendSize, MPI_UNSIGNED_LONG, sendRank, 46, client->intraComm, &requests.back());
210  }
211
212  status.resize(requests.size());
213  MPI_Waitall(requests.size(), &requests[0], &status[0]);
214
215
216  std::vector<MPI_Request>().swap(requests);
217  std::vector<MPI_Status>().swap(status);
218
219  // Okie, on destination side, we will wait for information of masked index of source
220  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
221  {
222    int recvRank = itSend->first;
223    int recvSize = itSend->second;
224
225    requests.push_back(MPI_Request());
226    MPI_Irecv(recvValueFromSrc[recvRank], recvSize, MPI_DOUBLE, recvRank, 48, client->intraComm, &requests.back());
227  }
228
229  for (std::map<int,int>::const_iterator itRecv = recvRankSizeMap.begin(); itRecv != recvRankSizeMap.end(); ++itRecv)
230  {
231    int recvRank = itRecv->first;
232    int recvSize = itRecv->second;
233    double* sendValue = sendValueToDest[recvRank];
234    unsigned long* recvIndexSrc = recvGlobalIndexSrc[recvRank];
235    int realSendSize = 0;
236    for (int idx = 0; idx < recvSize; ++idx)
237    {
238      size_t globalIndex = *(recvIndexSrc+idx);
239      int localIndex = globalIndex - ibeginSrc;
240      *(sendValue + idx) = axisSrc_->value(localIndex);
241    }
242    // Okie, now inform the destination which source index are masked
243    requests.push_back(MPI_Request());
244    MPI_Isend(sendValueToDest[recvRank], recvSize, MPI_DOUBLE, recvRank, 48, client->intraComm, &requests.back());
245  }
246  status.resize(requests.size());
247  MPI_Waitall(requests.size(), &requests[0], &status[0]);
248
249
250  size_t nGloAxisDest = axisDest_->n_glo.getValue() - 1;
251  for (std::map<int,int>::const_iterator itSend = sendRankSizeMap.begin(); itSend != sendRankSizeMap.end(); ++itSend)
252  {
253    int recvRank = itSend->first;
254    int recvSize = itSend->second;
255
256    double* recvValue = recvValueFromSrc[recvRank];
257    unsigned long* recvIndex = sendGlobalIndexSrc[recvRank];
258    for (int idx = 0; idx < recvSize; ++idx)
259    {
260      size_t globalIndex = *(recvIndex+idx);
261      int localIndex = ((nGloAxisDest-globalIndex) - axisDest_->begin);
262      axisDest_->value(localIndex) = *(recvValue + idx);
263    }
264  }
265
266  delete [] recvCount;
267  delete [] displ;
268  delete [] sendRankBuff;
269  delete [] recvRankBuff;
270  delete [] sendSizeBuff;
271  delete [] recvSizeBuff;
272
273  std::unordered_map<int, double* >::const_iterator itChar;
274  for (itChar = sendValueToDest.begin(); itChar != sendValueToDest.end(); ++itChar)
275    delete [] itChar->second;
276  for (itChar = recvValueFromSrc.begin(); itChar != recvValueFromSrc.end(); ++itChar)
277    delete [] itChar->second;
278  std::unordered_map<int, unsigned long* >::const_iterator itLong;
279  for (itLong = sendGlobalIndexSrc.begin(); itLong != sendGlobalIndexSrc.end(); ++itLong)
280    delete [] itLong->second;
281  for (itLong = recvGlobalIndexSrc.begin(); itLong != recvGlobalIndexSrc.end(); ++itLong)
282    delete [] itLong->second;
283}
284
285}
Note: See TracBrowser for help on using the repository browser.