source: XIOS/dev/dev_olga/src/transformation/axis_algorithm_inverse.cpp @ 1620

Last change on this file since 1620 was 1612, checked in by oabramkina, 2 years ago

Dev: adding exception handling.

To activate it, compilation flag -DXIOS_EXCEPTION should be added.

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