# source:XIOS/xios_training/XIOS_COMPILED/src/transformation/axis_algorithm_inverse.cpp@2101

Last change on this file since 2101 was 2101, checked in by yushan, 3 years ago

XIOS TRAINING some source updates for HO 11

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