src/ParaMEDMEM/OverlapMapping.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "OverlapMapping.hxx"
  22 #include "MPIProcessorGroup.hxx"
  23
  24 #include "MEDCouplingFieldDouble.hxx"
  25 #include "MCAuto.hxx"
  26
  27 #include "InterpKernelAutoPtr.hxx"
  28
  29 #include <numeric>
  30 #include <algorithm>
  31
  32 using namespace MEDCoupling;
  33
  34 OverlapMapping::OverlapMapping(const ProcessorGroup& group, const OverlapElementLocator & loc):
  35     _group(group),_locator(loc)
  36 {
  37 }
  38
  39 /*!
  40  * Keeps the link between a given a proc holding source mesh data, and the corresponding cell IDs.
  41  */
  42 void OverlapMapping::keepTracksOfSourceIds(int procId, DataArrayIdType *ids)
  43 {
  44   ids->incrRef();
  45   _sent_src_ids[procId] = ids;
  46 }
  47
  48 /*!
  49  * Same as keepTracksOfSourceIds() but for target mesh data.
  50 */
  51 void OverlapMapping::keepTracksOfTargetIds(int procId, DataArrayIdType *ids)
  52 {
  53   ids->incrRef();
  54   _sent_trg_ids[procId] = ids;
  55 }
  56
  57 /*!
  58  * This method stores in the local members the contribution coming from a matrix in format
  59  * Target(rows)/Source(cols) for a source procId 'srcProcId' and for a target procId 'trgProcId'.
  60  * All IDs received here (source and target) are in the format of local IDs.
  61  *
  62  * @param srcIds is null if the source mesh is on the local proc
  63  * @param trgIds is null if the source mesh is on the local proc
  64  *
  65  * One of the 2 is necessarily null (the two can be null together)
  66  */
  67 void OverlapMapping::addContributionST(const std::vector< SparseDoubleVec >& matrixST, const DataArrayIdType *srcIds, int srcProcId, const DataArrayIdType *trgIds, int trgProcId)
  68 {
  69   _matrixes_st.push_back(matrixST);
  70   _source_proc_id_st.push_back(srcProcId);
  71   _target_proc_id_st.push_back(trgProcId);
  72   if(srcIds)  // source mesh part is remote <=> srcProcId != myRank
  73       _nb_of_rcv_src_ids[srcProcId] = srcIds->getNumberOfTuples();
  74   else        // source mesh part is local
  75     {
  76       std::set<mcIdType> s;
  77       // For all source IDs (=col indices) in the sparse matrix:
  78       for(std::vector< SparseDoubleVec >::const_iterator it1=matrixST.begin();it1!=matrixST.end();it1++)
  79         for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
  80           s.insert((*it2).first);
  81       vector<mcIdType> v(s.begin(), s.end());  // turn set into vector
  82       _src_ids_zip_comp[trgProcId] = v;
  83     }
  84 }
  85
  86 /*!
  87  * This method is in charge to send matrices in AlltoAll mode.
  88  *
  89  * 'procsToSendField' gives the list of procs field data has to be sent to.
  90  * See OverlapElementLocator::computeBoundingBoxesAndTodoList()
  91  *
  92  * After the call of this method, 'this' contains the matrixST for all source cells of the current proc
  93  */
  94 void OverlapMapping::prepare(const std::vector< int >& procsToSendField, mcIdType nbOfTrgElems)
  95 {
  96 #ifdef DEC_DEBUG
  97   printMatrixesST();
  98 #endif
  99
 100   CommInterface commInterface=_group.getCommInterface();
 101   const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 102   const MPI_Comm *comm=group->getComm();
 103   std::size_t grpSize=_group.size();
 104   INTERP_KERNEL::AutoPtr<mcIdType> nbsend=new mcIdType[grpSize];
 105   INTERP_KERNEL::AutoPtr<int> nbsend2=new int[grpSize];
 106   INTERP_KERNEL::AutoPtr<int> nbsend3=new int[grpSize];
 107   std::fill<mcIdType *>(nbsend,nbsend+grpSize,0);
 108   int myProcId=_group.myRank();
 109   for(std::size_t i=0;i<_matrixes_st.size();i++)
 110     if(_source_proc_id_st[i]==myProcId)
 111       nbsend[_target_proc_id_st[i]]=(int)_matrixes_st[i].size();
 112   INTERP_KERNEL::AutoPtr<mcIdType> nbrecv=new mcIdType[grpSize];
 113   commInterface.allToAll(nbsend,1,MPI_ID_TYPE,nbrecv,1,MPI_ID_TYPE,*comm);
 114   //exchanging matrix
 115   //first exchanging offsets+ids_source
 116   INTERP_KERNEL::AutoPtr<int> nbrecv1=new int[grpSize];
 117   INTERP_KERNEL::AutoPtr<int> nbrecv2=new int[grpSize];
 118   //
 119   mcIdType *tmp=0;
 120   serializeMatrixStep0ST(nbrecv,
 121                          tmp,nbsend2,nbsend3,
 122                          nbrecv1,nbrecv2);
 123   INTERP_KERNEL::AutoPtr<mcIdType> bigArr=tmp;
 124   INTERP_KERNEL::AutoPtr<mcIdType> bigArrRecv=new mcIdType[nbrecv2[grpSize-1]+nbrecv1[grpSize-1]];
 125   commInterface.allToAllV(bigArr,nbsend2,nbsend3,MPI_ID_TYPE,
 126                           bigArrRecv,nbrecv1,nbrecv2,MPI_ID_TYPE,
 127                           *comm);// sending ids of sparse matrix (n+1 elems)
 128   //second phase echange target ids
 129   std::fill<int *>(nbsend2,nbsend2+grpSize,0);
 130   INTERP_KERNEL::AutoPtr<int> nbrecv3=new int[grpSize];
 131   INTERP_KERNEL::AutoPtr<int> nbrecv4=new int[grpSize];
 132   double *tmp2=0;
 133   mcIdType lgthOfArr=serializeMatrixStep1ST(nbrecv,bigArrRecv,nbrecv1,nbrecv2,
 134                                        tmp,tmp2,
 135                                        nbsend2,nbsend3,nbrecv3,nbrecv4);
 136   INTERP_KERNEL::AutoPtr<mcIdType> bigArr2=tmp;
 137   INTERP_KERNEL::AutoPtr<double> bigArrD2=tmp2;
 138   INTERP_KERNEL::AutoPtr<mcIdType> bigArrRecv2=new mcIdType[lgthOfArr];
 139   INTERP_KERNEL::AutoPtr<double> bigArrDRecv2=new double[lgthOfArr];
 140   commInterface.allToAllV(bigArr2,nbsend2,nbsend3,MPI_ID_TYPE,
 141                           bigArrRecv2,nbrecv3,nbrecv4,MPI_ID_TYPE,
 142                           *comm);
 143   commInterface.allToAllV(bigArrD2,nbsend2,nbsend3,MPI_DOUBLE,
 144                           bigArrDRecv2,nbrecv3,nbrecv4,MPI_DOUBLE,
 145                           *comm);
 146   //finishing
 147   unserializationST(nbOfTrgElems,nbrecv,bigArrRecv,nbrecv1,nbrecv2,
 148                     bigArrRecv2,bigArrDRecv2,nbrecv3,nbrecv4);
 149
 150   //finish to fill _the_matrix_st with already in place matrix in _matrixes_st (local computation)
 151   finishToFillFinalMatrixST();
 152
 153   //updating _src_ids_zip_st2 and _src_ids_zip_st2 with received matrix.
 154   fillSourceIdsZipReceivedForMultiply();
 155   // Prepare proc list for future field data exchange (multiply()):
 156   _proc_ids_to_send_vector_st = procsToSendField;
 157   // Make some space on local proc:
 158   _matrixes_st.clear();
 159
 160 #ifdef DEC_DEBUG
 161   printTheMatrix();
 162 #endif
 163 }
 164
 165 ///*!
 166 // * Compute denominators for ExtensiveConservation interp.
 167 // * TO BE REVISED: needs another communication since some bits are held non locally
 168 // */
 169 //void OverlapMapping::computeDenoGlobConstraint()
 170 //{
 171 //  _the_deno_st.clear();
 172 //  std::size_t sz1=_the_matrix_st.size();
 173 //  _the_deno_st.resize(sz1);
 174 //  for(std::size_t i=0;i<sz1;i++)
 175 //    {
 176 //      std::size_t sz2=_the_matrix_st[i].size();
 177 //      _the_deno_st[i].resize(sz2);
 178 //      for(std::size_t j=0;j<sz2;j++)
 179 //        {
 180 //          double sum=0;
 181 //          SparseDoubleVec& mToFill=_the_deno_st[i][j];
 182 //          const SparseDoubleVec& m=_the_matrix_st[i][j];
 183 //          for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++)
 184 //            sum+=(*it).second;
 185 //          for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++)
 186 //            mToFill[(*it).first]=sum;
 187 //        }
 188 //    }
 189 //    printDenoMatrix();
 190 //}
 191
 192 ///*! Compute integral denominators
 193 // * TO BE REVISED: needs another communication since some source areas are held non locally
 194 // */
 195 //void OverlapMapping::computeDenoIntegral()
 196 //{
 197 //  _the_deno_st.clear();
 198 //  std::size_t sz1=_the_matrix_st.size();
 199 //  _the_deno_st.resize(sz1);
 200 //  for(std::size_t i=0;i<sz1;i++)
 201 //    {
 202 //      std::size_t sz2=_the_matrix_st[i].size();
 203 //      _the_deno_st[i].resize(sz2);
 204 //      for(std::size_t j=0;j<sz2;j++)
 205 //        {
 206 //          SparseDoubleVec& mToFill=_the_deno_st[i][j];
 207 //          for(SparseDoubleVec::const_iterator it=mToFill.begin();it!=mToFill.end();it++)
 208 //            mToFill[(*it).first] = sourceAreas;
 209 //        }
 210 //    }
 211 //    printDenoMatrix();
 212 //}
 213
 214 /*! Compute rev integral denominators
 215   */
 216 void OverlapMapping::computeDenoRevIntegral(const DataArrayDouble & targetAreas)
 217 {
 218   _the_deno_st.clear();
 219   std::size_t sz1=_the_matrix_st.size();
 220   _the_deno_st.resize(sz1);
 221   const double * targetAreasP = targetAreas.getConstPointer();
 222   for(std::size_t i=0;i<sz1;i++)
 223     {
 224       std::size_t sz2=_the_matrix_st[i].size();
 225       _the_deno_st[i].resize(sz2);
 226       for(std::size_t j=0;j<sz2;j++)
 227         {
 228           SparseDoubleVec& mToFill=_the_deno_st[i][j];
 229           SparseDoubleVec& mToIterate=_the_matrix_st[i][j];
 230           for(SparseDoubleVec::const_iterator it=mToIterate.begin();it!=mToIterate.end();it++)
 231             mToFill[(*it).first] = targetAreasP[j];
 232         }
 233     }
 234 //    printDenoMatrix();
 235 }
 236
 237
 238 /*!
 239  * Compute denominators for ConvervativeVolumic interp.
 240  */
 241 void OverlapMapping::computeDenoConservativeVolumic(mcIdType nbOfTuplesTrg)
 242 {
 243   int myProcId=_group.myRank();
 244   //
 245   _the_deno_st.clear();
 246   std::size_t sz1=_the_matrix_st.size();
 247   _the_deno_st.resize(sz1);
 248   std::vector<double> deno(nbOfTuplesTrg);
 249   // Fills in the vector indexed by target cell ID:
 250   for(std::size_t i=0;i<sz1;i++)
 251     {
 252       const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i];
 253       int curSrcId=_the_matrix_st_source_proc_id[i];
 254       map < int, MCAuto<DataArrayIdType> >::const_iterator isItem1 = _sent_trg_ids.find(curSrcId);
 255       mcIdType rowId=0;
 256       if(isItem1==_sent_trg_ids.end() || curSrcId==myProcId) // Local computation: simple, because rowId of mat are directly target cell ids.
 257         {
 258           for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
 259             for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 260               deno[rowId]+=(*it2).second;
 261         }
 262       else  // matrix was received, remote computation
 263         {
 264           const DataArrayIdType *trgIds = (*isItem1).second;
 265           const mcIdType *trgIds2=trgIds->getConstPointer();
 266           for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
 267             for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 268               deno[trgIds2[rowId]]+=(*it2).second;
 269         }
 270     }
 271   // Broadcast the vector into a structure similar to the initial sparse matrix of numerators:
 272   for(std::size_t i=0;i<sz1;i++)
 273     {
 274       mcIdType rowId=0;
 275       const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i];
 276       int curSrcId=_the_matrix_st_source_proc_id[i];
 277       map < int, MCAuto<DataArrayIdType> >::const_iterator isItem1 = _sent_trg_ids.find(curSrcId);
 278       std::vector< SparseDoubleVec >& denoM=_the_deno_st[i];
 279       denoM.resize(mat.size());
 280       if(isItem1==_sent_trg_ids.end() || curSrcId==myProcId)//item1 of step2 main algo. Simple, because rowId of mat are directly target ids.
 281         {
 282           mcIdType rowId2=0;
 283           for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId2++)
 284             for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 285               denoM[rowId2][(*it2).first]=deno[rowId2];
 286         }
 287       else
 288         {
 289           const DataArrayIdType *trgIds = (*isItem1).second;
 290           const mcIdType *trgIds2=trgIds->getConstPointer();
 291           for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
 292             for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 293               denoM[rowId][(*it2).first]=deno[trgIds2[rowId]];
 294         }
 295     }
 296 //  printDenoMatrix();
 297 }
 298
 299 /*!
 300  * This method performs step #0/3 in serialization process.
 301  * \param count tells specifies nb of elems to send to corresponding proc id. size equal to _group.size().
 302  * \param offsets tells for a proc i where to start serialize#0 matrix. size equal to _group.size().
 303  * \param nbOfElemsSrc of size _group.size(). Comes from previous all2all call. tells how many srcIds per proc contains matrix for current proc.
 304  */
 305 void OverlapMapping::serializeMatrixStep0ST(const mcIdType *nbOfElemsSrc, mcIdType *&bigArr, int *count, int *offsets,
 306                                             int *countForRecv, int *offsetsForRecv) const
 307 {
 308   std::size_t grpSize=_group.size();
 309   std::fill<int *>(count,count+grpSize,0);
 310   std::size_t szz=0;
 311   int myProcId=_group.myRank();
 312   for(std::size_t i=0;i<_matrixes_st.size();i++)
 313     {
 314       if(_source_proc_id_st[i]==myProcId && _matrixes_st[i].size())// && _target_proc_id_st[i]!=myProcId
 315         {
 316           count[_target_proc_id_st[i]]=(int)_matrixes_st[i].size()+1;
 317           szz+=_matrixes_st[i].size()+1;
 318         }
 319     }
 320   bigArr=new mcIdType[szz];
 321   offsets[0]=0;
 322   for(std::size_t i=1;i<grpSize;i++)
 323     offsets[i]=offsets[i-1]+count[i-1];
 324   for(std::size_t i=0;i<_matrixes_st.size();i++)
 325     {
 326       if(_source_proc_id_st[i]==myProcId && _matrixes_st[i].size())
 327         {
 328           mcIdType start=offsets[_target_proc_id_st[i]];
 329           mcIdType *work=bigArr+start;
 330           *work=0;
 331           const std::vector< SparseDoubleVec >& mat=_matrixes_st[i];
 332           for(std::vector< SparseDoubleVec >::const_iterator it=mat.begin();it!=mat.end();it++,work++)
 333             work[1]=work[0]+ToIdType((*it).size());
 334         }
 335     }
 336   //
 337   offsetsForRecv[0]=0;
 338   for(std::size_t i=0;i<grpSize;i++)
 339     {
 340       if(nbOfElemsSrc[i]>0)
 341         countForRecv[i]=(int)nbOfElemsSrc[i]+1;
 342       else
 343         countForRecv[i]=0;
 344       if(i>0)
 345         offsetsForRecv[i]=offsetsForRecv[i-1]+countForRecv[i-1];
 346     }
 347 }
 348
 349 /*!
 350  * This method performs step#1 and step#2/3. It returns the size of expected array to get allToAllV.
 351  * It is where the locally computed matrices are serialized to be sent to adequate final proc.
 352  */
 353 mcIdType OverlapMapping::serializeMatrixStep1ST(const mcIdType *nbOfElemsSrc, const mcIdType *recvStep0, const int *countStep0, const int *offsStep0,
 354                                            mcIdType *&bigArrI, double *&bigArrD, int *count, int *offsets,
 355                                            int *countForRecv, int *offsForRecv) const
 356 {
 357   std::size_t grpSize=_group.size();
 358   int myProcId=_group.myRank();
 359   offsForRecv[0]=0;
 360   mcIdType szz=0;
 361   for(std::size_t i=0;i<grpSize;i++)
 362     {
 363       if(nbOfElemsSrc[i]!=0)
 364         countForRecv[i]=(int)recvStep0[offsStep0[i]+nbOfElemsSrc[i]];
 365       else
 366         countForRecv[i]=0;
 367       szz+=countForRecv[i];
 368       if(i>0)
 369         offsForRecv[i]=offsForRecv[i-1]+countForRecv[i-1];
 370     }
 371   //
 372   std::fill(count,count+grpSize,0);
 373   offsets[0]=0;
 374   std::size_t fullLgth=0;
 375   for(std::size_t i=0;i<_matrixes_st.size();i++)
 376     {
 377       if(_source_proc_id_st[i]==myProcId && _matrixes_st[i].size())
 378         {
 379           const std::vector< SparseDoubleVec >& mat=_matrixes_st[i];
 380           mcIdType lgthToSend=0;
 381           for(std::vector< SparseDoubleVec >::const_iterator it=mat.begin();it!=mat.end();it++)
 382             lgthToSend+=ToIdType((*it).size());
 383           count[_target_proc_id_st[i]]=(int)lgthToSend;
 384           fullLgth+=lgthToSend;
 385         }
 386     }
 387   for(std::size_t i=1;i<grpSize;i++)
 388     offsets[i]=offsets[i-1]+count[i-1];
 389   //
 390   bigArrI=new mcIdType[fullLgth];
 391   bigArrD=new double[fullLgth];
 392   // feeding arrays
 393   fullLgth=0;
 394   for(std::size_t i=0;i<_matrixes_st.size();i++)
 395     {
 396       if(_source_proc_id_st[i]==myProcId && _matrixes_st[i].size())
 397         {
 398           const std::vector< SparseDoubleVec >& mat=_matrixes_st[i];
 399           for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++)
 400             {
 401               mcIdType j=0;
 402               for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++,j++)
 403                 {
 404                   bigArrI[fullLgth+j]=(*it2).first;
 405                   bigArrD[fullLgth+j]=(*it2).second;
 406                 }
 407               fullLgth+=(*it1).size();
 408             }
 409         }
 410     }
 411   return szz;
 412 }
 413
 414 /*!
 415  * This is the last step after all2Alls for matrix exchange.
 416  * _the_matrix_st is the final matrix :
 417  *      - The first entry is srcId in current proc.
 418  *      - The second is the pseudo id of source proc (correspondence with true id is in attribute _the_matrix_st_source_proc_id and _the_matrix_st_source_ids)
 419  *      - the third is the srcId in the pseudo source proc
 420  */
 421 void OverlapMapping::unserializationST(mcIdType nbOfTrgElems,
 422                                        const mcIdType *nbOfElemsSrcPerProc,//first all2all
 423                                        const mcIdType *bigArrRecv, const int *bigArrRecvCounts, const int *bigArrRecvOffs,//2nd all2all
 424                                        const mcIdType *bigArrRecv2, const double *bigArrDRecv2, const int *bigArrRecv2Count, const int *bigArrRecv2Offs)//3rd and 4th all2alls
 425 {
 426   _the_matrix_st.clear();
 427   _the_matrix_st_source_proc_id.clear();
 428   //
 429   std::size_t grpSize=_group.size();
 430   for(unsigned int i=0;i<grpSize;i++)
 431     if(nbOfElemsSrcPerProc[i]!=0)
 432       _the_matrix_st_source_proc_id.push_back(i);
 433   std::size_t nbOfPseudoProcs=_the_matrix_st_source_proc_id.size();//_the_matrix_st_target_proc_id.size() contains number of matrix fetched remotely whose sourceProcId==myProcId
 434   _the_matrix_st.resize(nbOfPseudoProcs);
 435   //
 436   std::size_t j=0;
 437   for(std::size_t i=0;i<grpSize;i++)
 438     if(nbOfElemsSrcPerProc[i]!=0)
 439       {
 440         _the_matrix_st[j].resize(nbOfElemsSrcPerProc[i]);
 441         for(mcIdType k=0;k<nbOfElemsSrcPerProc[i];k++)
 442           {
 443             mcIdType offs=bigArrRecv[bigArrRecvOffs[i]+k];
 444             mcIdType lgthOfMap=bigArrRecv[bigArrRecvOffs[i]+k+1]-offs;
 445             for(mcIdType l=0;l<lgthOfMap;l++)
 446               _the_matrix_st[j][k][bigArrRecv2[bigArrRecv2Offs[i]+offs+l]]=bigArrDRecv2[bigArrRecv2Offs[i]+offs+l];
 447           }
 448         j++;
 449       }
 450 }
 451
 452 /*!
 453  * This method should be called when all remote matrix with sourceProcId==thisProcId have been retrieved and are
 454  * in 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id' and 'this->_the_matrix_st_target_ids'.
 455  * This method finish the job of filling 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id'
 456  * by putting candidates in 'this->_matrixes_st' into them (i.e. local computation result).
 457  */
 458 void OverlapMapping::finishToFillFinalMatrixST()
 459 {
 460   int myProcId=_group.myRank();
 461   std::size_t sz=_matrixes_st.size();
 462   int nbOfEntryToAdd=0;
 463   for(std::size_t i=0;i<sz;i++)
 464     if(_source_proc_id_st[i]!=myProcId)
 465       nbOfEntryToAdd++;
 466   if(nbOfEntryToAdd==0)
 467     return ;
 468   std::size_t oldNbOfEntry=_the_matrix_st.size();
 469   std::size_t newNbOfEntry=oldNbOfEntry+nbOfEntryToAdd;
 470   _the_matrix_st.resize(newNbOfEntry);
 471   std::size_t j=oldNbOfEntry;
 472   for(std::size_t i=0;i<sz;i++)
 473     if(_source_proc_id_st[i]!=myProcId)
 474       {
 475         const std::vector<SparseDoubleVec >& mat=_matrixes_st[i];
 476         _the_matrix_st[j]=mat;
 477         _the_matrix_st_source_proc_id.push_back(_source_proc_id_st[i]);
 478         j++;
 479       }
 480 }
 481
 482
 483 /*!
 484  * This method performs a transpose multiply of 'fieldInput' and put the result into 'fieldOutput'.
 485  * 'fieldInput' is expected to be the sourcefield and 'fieldOutput' the targetfield.
 486  */
 487 void OverlapMapping::multiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput, double default_val) const
 488 {
 489   using namespace std;
 490
 491   std::size_t nbOfCompo=fieldInput->getNumberOfComponents();//to improve same number of components to test
 492   CommInterface commInterface=_group.getCommInterface();
 493   const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 494   const MPI_Comm *comm=group->getComm();
 495   int grpSize=_group.size();
 496   int myProcID=_group.myRank();
 497   //
 498   INTERP_KERNEL::AutoPtr<int> nbsend=new int[grpSize];
 499   INTERP_KERNEL::AutoPtr<int> nbsend2=new int[grpSize];
 500   INTERP_KERNEL::AutoPtr<int> nbrecv=new int[grpSize];
 501   INTERP_KERNEL::AutoPtr<int> nbrecv2=new int[grpSize];
 502   fill<int *>(nbsend,nbsend+grpSize,0);
 503   fill<int *>(nbrecv,nbrecv+grpSize,0);
 504   nbsend2[0]=0;
 505   nbrecv2[0]=0;
 506   vector<double> valsToSend;
 507
 508   /*
 509    * FIELD VALUE XCHGE:
 510    * We call the 'BB source IDs' (bounding box source IDs) the set of source cell IDs transmitted just based on the bounding box information.
 511    * This is potentially bigger than what is finally in the interp matrix and this is stored in _sent_src_ids.
 512    * We call 'interp source IDs' the set of source cell IDs with non null entries in the interp matrix. This is a sub-set of the above.
 513    */
 514   for(int procID=0;procID<grpSize;procID++)
 515     {
 516       /* SENDING part: compute field values to be SENT (and how many of them)
 517        *   - for all proc 'procID' in group
 518        *      * if procID == myProcID, send nothing
 519        *      * elif 'procID' in _proc_ids_to_send_vector_st (computed from the BB intersection)
 520        *        % if myProcID computed the job (myProcID, procID)
 521        *           => send only 'interp source IDs' field values (i.e. IDs stored in _src_ids_zip_comp)
 522        *        % else (=we just sent mesh data to procID, but have never seen the matrix, i.e. matrix was computed remotely by procID)
 523        *           => send 'BB source IDs' set of field values (i.e. IDs stored in _sent_src_ids)
 524        */
 525       if (procID == myProcID)
 526         nbsend[procID] = 0;
 527       else
 528         if(find(_proc_ids_to_send_vector_st.begin(),_proc_ids_to_send_vector_st.end(),procID)!=_proc_ids_to_send_vector_st.end())
 529           {
 530             MCAuto<DataArrayDouble> vals;
 531             if(_locator.isInMyTodoList(myProcID, procID))
 532               {
 533                 map<int, vector<mcIdType> >::const_iterator isItem11 = _src_ids_zip_comp.find(procID);
 534                 if (isItem11 == _src_ids_zip_comp.end())
 535                   throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: SEND: unexpected end iterator in _src_ids_zip_comp!");
 536                 const vector<mcIdType> & v = (*isItem11).second;
 537                 std::size_t sz = v.size();
 538                 vals=fieldInput->getArray()->selectByTupleId(&(v[0]),&(v[0])+sz);
 539               }
 540             else
 541               {
 542                 map < int, MCAuto<DataArrayIdType> >::const_iterator isItem11 = _sent_src_ids.find( procID );
 543                 if (isItem11 == _sent_src_ids.end())
 544                   throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: SEND: unexpected end iterator in _sent_src_ids!");
 545                 vals=fieldInput->getArray()->selectByTupleId(*(*isItem11).second);
 546               }
 547             nbsend[procID] = (int)vals->getNbOfElems();  // nb of elem = nb_tuples*nb_compo
 548             // Flat version of values to send:
 549             valsToSend.insert(valsToSend.end(),vals->getConstPointer(),vals->getConstPointer()+nbsend[procID]);
 550           }
 551
 552       /* RECEIVE: compute number of field values to be RECEIVED
 553        *   - for all proc 'procID' in group
 554        *      * if procID == myProcID, rcv nothing
 555        *      * elif 'procID' in _proc_ids_to_recv_vector_st (computed from BB intersec)
 556        *        % if myProcID computed the job (procID, myProcID)
 557        *          => receive full set ('BB source IDs') of field data from proc #procID which has never seen the matrix
 558        *             i.e. prepare to receive the numb in _nb_of_rcv_src_ids
 559        *        % else (=we did NOT compute the job, hence procID has, and knows the matrix)
 560        *          => receive 'interp source IDs' set of field values
 561        */
 562       const std::vector< int > & _proc_ids_to_recv_vector_st = _the_matrix_st_source_proc_id;
 563       if (procID == myProcID)
 564         nbrecv[procID] = 0;
 565       else
 566         if(find(_proc_ids_to_recv_vector_st.begin(),_proc_ids_to_recv_vector_st.end(),procID)!=_proc_ids_to_recv_vector_st.end())
 567           {
 568             if(_locator.isInMyTodoList(procID, myProcID))
 569               {
 570                 map <int,mcIdType>::const_iterator isItem11 = _nb_of_rcv_src_ids.find(procID);
 571                 if (isItem11 == _nb_of_rcv_src_ids.end())
 572                   throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: RCV: unexpected end iterator in _nb_of_rcv_src_ids!");
 573                 nbrecv[procID] = (int)((*isItem11).second * nbOfCompo);
 574               }
 575             else
 576               {
 577                 map<int, vector<mcIdType> >::const_iterator isItem22 = _src_ids_zip_recv.find(procID);
 578                 if (isItem22 == _src_ids_zip_recv.end())
 579                   throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: RCV: unexpected end iterator in _src_ids_zip_recv!");
 580                 nbrecv[procID] = (int)((*isItem22).second.size() * nbOfCompo);
 581               }
 582           }
 583     }
 584   // Compute offsets in the sending/receiving array.
 585   for(int i=1;i<grpSize;i++)
 586     {
 587       nbsend2[i]=nbsend2[i-1]+nbsend[i-1];
 588       nbrecv2[i]=nbrecv2[i-1]+nbrecv[i-1];
 589     }
 590   INTERP_KERNEL::AutoPtr<double> bigArr=new double[nbrecv2[grpSize-1]+nbrecv[grpSize-1]];
 591
 592 #ifdef DEC_DEBUG
 593   stringstream scout;
 594   scout << "("  << myProcID << ") nbsend :" << nbsend[0] << "," << nbsend[1] << "," << nbsend[2] << "\n";
 595   scout << "("  << myProcID << ") nbrecv :" << nbrecv[0] << "," << nbrecv[1] << "," << nbrecv[2] << "\n";
 596   scout << "("  << myProcID << ") valsToSend: ";
 597   for (int iii=0; iii<valsToSend.size(); iii++)
 598     scout << ", " << valsToSend[iii];
 599   cout << scout.str() << "\n";
 600 #endif
 601
 602   /*
 603    * *********************** ALL-TO-ALL
 604    */
 605   commInterface.allToAllV(&valsToSend[0],nbsend,nbsend2,MPI_DOUBLE,
 606                           bigArr,nbrecv,nbrecv2,MPI_DOUBLE,*comm);
 607 #ifdef DEC_DEBUG
 608   MPI_Barrier(MPI_COMM_WORLD);
 609   scout << "("  << myProcID << ") bigArray: ";
 610     for (int iii=0; iii<nbrecv2[grpSize-1]+nbrecv[grpSize-1]; iii++)
 611       scout << ", " << bigArr[iii];
 612   cout << scout.str() << "\n";
 613 #endif
 614
 615   /*
 616    * TARGET FIELD COMPUTATION (matrix-vec computation)
 617    */
 618   fieldOutput->getArray()->fillWithZero();
 619   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfCompo];
 620
 621   // By default field value set to default value - so mark which cells are hit
 622   mcIdType ntup = fieldOutput->getNumberOfTuples();
 623   INTERP_KERNEL::AutoPtr<bool> hit_cells = new bool[ntup];
 624   std::fill((bool *)hit_cells, (bool *)hit_cells+ntup, false);
 625
 626   for(vector<int>::const_iterator itProc=_the_matrix_st_source_proc_id.begin(); itProc != _the_matrix_st_source_proc_id.end();itProc++)
 627   // For each source processor corresponding to a locally held matrix:
 628     {
 629       int srcProcID = *itProc;
 630       std::size_t id = std::distance(_the_matrix_st_source_proc_id.begin(),itProc);
 631       const vector< SparseDoubleVec >& mat =_the_matrix_st[id];
 632       const vector< SparseDoubleVec >& deno = _the_deno_st[id];
 633
 634       /*   FINAL MULTIPLICATION
 635        *      * if srcProcID == myProcID, local multiplication without any mapping
 636        *         => for all target cell ID 'tgtCellID'
 637        *           => for all src cell ID 'srcCellID' in the sparse vector
 638        *             => tgtFieldLocal[tgtCellID] += srcFieldLocal[srcCellID] * matrix[tgtCellID][srcCellID] / deno[tgtCellID][srcCellID]
 639        */
 640       if (srcProcID == myProcID)
 641         {
 642           std::size_t nbOfTrgTuples=mat.size();
 643           double * targetBase = fieldOutput->getArray()->getPointer();
 644           for(std::size_t j=0; j<nbOfTrgTuples; j++)
 645             {
 646               const SparseDoubleVec& mat1=mat[j];
 647               const SparseDoubleVec& deno1=deno[j];
 648               SparseDoubleVec::const_iterator it5=deno1.begin();
 649               const double * localSrcField = fieldInput->getArray()->getConstPointer();
 650               double * targetPt = targetBase+j*nbOfCompo;
 651               for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
 652                 {
 653                   // Apply the multiplication for all components:
 654                   double ratio = (*it3).second/(*it5).second;
 655                   transform(localSrcField+((*it3).first)*nbOfCompo,
 656                             localSrcField+((*it3).first+1)*nbOfCompo,
 657                             (double *)tmp,
 658                             [=](double d) { return d*ratio; });
 659                   // Accumulate with current value:
 660                   transform((double *)tmp,(double *)tmp+nbOfCompo,targetPt,targetPt,plus<double>());
 661                   hit_cells[j] = true;
 662                 }
 663             }
 664         }
 665
 666       if(nbrecv[srcProcID]<=0)  // also covers the preceding 'if'
 667         continue;
 668
 669       /*      * if something was received
 670        *         %  if received matrix (=we didn't compute the job), this means that :
 671        *            1. we sent part of our targetIDs to srcProcID before, so that srcProcId can do the computation.
 672        *            2. srcProcID has sent us only the 'interp source IDs' field values
 673        *            => invert _src_ids_zip_recv -> 'revert_zip'
 674        *            => for all target cell ID 'tgtCellID'
 675        *              => mappedTgtID = _sent_trg_ids[srcProcID][tgtCellID]
 676        *              => for all src cell ID 'srcCellID' in the sparse vector
 677        *                 => idx = revert_zip[srcCellID]
 678        *                 => tgtFieldLocal[mappedTgtID] += rcvValue[srcProcID][idx] * matrix[tgtCellID][srcCellID] / deno[tgtCellID][srcCellID]
 679        */
 680       if(!_locator.isInMyTodoList(srcProcID, myProcID))
 681         {
 682           // invert _src_ids_zip_recv
 683           map<mcIdType,int> revert_zip;
 684           map<int, vector<mcIdType> >::const_iterator it11= _src_ids_zip_recv.find(srcProcID);
 685           if (it11 == _src_ids_zip_recv.end())
 686             throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: MULTIPLY: unexpected end iterator in _src_ids_zip_recv!");
 687
 688           const vector<mcIdType> & vec = (*it11).second;
 689           int newId=0;
 690           for(vector<mcIdType>::const_iterator it=vec.begin();it!=vec.end();it++,newId++)
 691             revert_zip[*it]=newId;
 692           map < int, MCAuto<DataArrayIdType> >::const_iterator isItem24 = _sent_trg_ids.find(srcProcID);
 693           if (isItem24 == _sent_trg_ids.end())
 694             throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: MULTIPLY: unexpected end iterator in _sent_trg_ids!");
 695           const DataArrayIdType *tgrIdsDA = (*isItem24).second;
 696           const mcIdType *tgrIds = tgrIdsDA->getConstPointer();
 697
 698           std::size_t nbOfTrgTuples=mat.size();
 699           double * targetBase = fieldOutput->getArray()->getPointer();
 700           for(std::size_t j=0;j<nbOfTrgTuples;j++)
 701             {
 702               const SparseDoubleVec& mat1=mat[j];
 703               const SparseDoubleVec& deno1=deno[j];
 704               SparseDoubleVec::const_iterator it5=deno1.begin();
 705               double * targetPt = targetBase+tgrIds[j]*nbOfCompo;
 706               for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
 707                 {
 708                   map<mcIdType,int>::const_iterator it4=revert_zip.find((*it3).first);
 709                   if(it4==revert_zip.end())
 710                     throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: MULTIPLY: unexpected end iterator in revert_zip!");
 711                   double ratio = (*it3).second/(*it5).second;
 712                   transform(bigArr+nbrecv2[srcProcID]+((*it4).second)*nbOfCompo,
 713                             bigArr+nbrecv2[srcProcID]+((*it4).second+1)*nbOfCompo,
 714                             (double *)tmp,
 715                             [=](double d) { return d*ratio; } );
 716                   transform((double *)tmp,(double *)tmp+nbOfCompo,targetPt,targetPt,plus<double>());
 717                   hit_cells[tgrIds[j]] = true;
 718                 }
 719             }
 720         }
 721       else
 722         /*         % else (=we computed the job and we received the 'BB source IDs' set of source field values)
 723          *            => for all target cell ID 'tgtCellID'
 724          *              => for all src cell ID 'srcCellID' in the sparse vector
 725          *                => tgtFieldLocal[tgtCellID] += rcvValue[srcProcID][srcCellID] * matrix[tgtCellID][srcCellID] / deno[tgtCellID][srcCellID]
 726          */
 727         {
 728           // Same loop as in the case srcProcID == myProcID, except that instead of working on local field data, we work on bigArr
 729           std::size_t nbOfTrgTuples=mat.size();
 730           double * targetBase = fieldOutput->getArray()->getPointer();
 731           for(std::size_t j=0;j<nbOfTrgTuples;j++)
 732             {
 733               const SparseDoubleVec& mat1=mat[j];
 734               const SparseDoubleVec& deno1=deno[j];
 735               SparseDoubleVec::const_iterator it5=deno1.begin();
 736               double * targetPt = targetBase+j*nbOfCompo;
 737               for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
 738                 {
 739                   // Apply the multiplication for all components:
 740                   double ratio = (*it3).second/(*it5).second;
 741                   transform(bigArr+nbrecv2[srcProcID]+((*it3).first)*nbOfCompo,
 742                             bigArr+nbrecv2[srcProcID]+((*it3).first+1)*nbOfCompo,
 743                             (double *)tmp,
 744                             [=](double d) { return d*ratio; } );
 745                   // Accumulate with current value:
 746                   transform((double *)tmp,(double *)tmp+nbOfCompo,targetPt,targetPt,plus<double>());
 747                   hit_cells[j] = true;
 748                 }
 749             }
 750         }
 751     }
 752
 753   // Fill in default values for cells which haven't been hit:
 754   int i = 0;
 755   for(bool * hit_cells_ptr=hit_cells; i< ntup; hit_cells_ptr++,i++)
 756     if (!(*hit_cells_ptr))
 757       {
 758         double * targetPt=fieldOutput->getArray()->getPointer();
 759         fill(targetPt+i*nbOfCompo, targetPt+(i+1)*nbOfCompo, default_val);
 760       }
 761 }
 762
 763 /*!
 764  * This method performs a transpose multiply of 'fieldInput' and put the result into 'fieldOutput'.
 765  * 'fieldInput' is expected to be the targetfield and 'fieldOutput' the sourcefield.
 766  */
 767 void OverlapMapping::transposeMultiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput)
 768 {
 769 }
 770
 771 /*!
 772  * This method should be called immediately after _the_matrix_st has been filled with remote computed matrix
 773  * put in this proc for Matrix-Vector.
 774  * It fills _src_ids_zip_recv (see member doc)
 775  */
 776 void OverlapMapping::fillSourceIdsZipReceivedForMultiply()
 777 {
 778   /* When it is called, only the bits received from other processors (i.e. the remotely executed jobs) are in the
 779     big matrix _the_matrix_st. */
 780
 781   CommInterface commInterface=_group.getCommInterface();
 782   int myProcId=_group.myRank();
 783   std::size_t nbOfMatrixRecveived=_the_matrix_st_source_proc_id.size();
 784   for(std::size_t i=0;i<nbOfMatrixRecveived;i++)
 785     {
 786       int curSrcProcId=_the_matrix_st_source_proc_id[i];
 787       if(curSrcProcId!=myProcId)  // if =, data has been populated by addContributionST()
 788         {
 789           const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i];
 790           std::set<mcIdType> s;
 791           for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++)
 792             for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 793               s.insert((*it2).first);
 794           vector<mcIdType> vec(s.begin(),s.end());
 795           _src_ids_zip_recv[curSrcProcId] = vec;
 796         }
 797     }
 798 }
 799
 800 #ifdef DEC_DEBUG
 801  void OverlapMapping::printTheMatrix() const
 802  {
 803    CommInterface commInterface=_group.getCommInterface();
 804    const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 805    const MPI_Comm *comm=group->getComm();
 806    int grpSize=_group.size();
 807    int myProcId=_group.myRank();
 808    std::stringstream oscerr;
 809    int nbOfMat=_the_matrix_st.size();
 810    oscerr << "(" <<  myProcId <<  ") I hold " << nbOfMat << " matrix(ces) : "<< std::endl;
 811    for(int i=0;i<nbOfMat;i++)
 812      {
 813        oscerr << "   - Matrix #" << i << " coming from source proc #" << _the_matrix_st_source_proc_id[i] << ":\n ";
 814        const std::vector< SparseDoubleVec >& locMat=_the_matrix_st[i];
 815        int j = 0;
 816        for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++)
 817          {
 818            oscerr << " Target Cell #" << j;
 819            for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 820              oscerr << " (" << (*it2).first << "," << (*it2).second << "), ";
 821            oscerr << std::endl;
 822          }
 823      }
 824    oscerr << "*********" << std::endl;
 825
 826    // Hope this will be flushed in one go:
 827    std::cerr << oscerr.str() << std::endl;
 828 //   if(myProcId != 0)
 829 //     MPI_Barrier(MPI_COMM_WORLD);
 830  }
 831
 832  void OverlapMapping::printMatrixesST() const
 833   {
 834     CommInterface commInterface=_group.getCommInterface();
 835     const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 836     const MPI_Comm *comm=group->getComm();
 837     int grpSize=_group.size();
 838     int myProcId=_group.myRank();
 839     std::stringstream oscerr;
 840     int nbOfMat=_matrixes_st.size();
 841     oscerr << "(" <<  myProcId <<  ") I hold " << nbOfMat << " LOCAL matrix(ces) : "<< std::endl;
 842     for(int i=0;i<nbOfMat;i++)
 843       {
 844         oscerr << "   - Matrix #" << i << ": (source proc #" << _source_proc_id_st[i] << " / tgt proc#" << _target_proc_id_st[i] << "): \n";
 845         const std::vector< SparseDoubleVec >& locMat=_matrixes_st[i];
 846         int j = 0;
 847         for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++)
 848           {
 849             oscerr << " Target Cell #" << j;
 850             for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 851               oscerr << " (" << (*it2).first << "," << (*it2).second << "), ";
 852             oscerr << std::endl;
 853           }
 854       }
 855     oscerr << "*********" << std::endl;
 856
 857     // Hope this will be flushed in one go:
 858     std::cerr << oscerr.str() << std::endl;
 859   }
 860
 861  void OverlapMapping::printDenoMatrix() const
 862    {
 863      CommInterface commInterface=_group.getCommInterface();
 864      const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 865      const MPI_Comm *comm=group->getComm();
 866      int grpSize=_group.size();
 867      int myProcId=_group.myRank();
 868      std::stringstream oscerr;
 869      int nbOfMat=_the_deno_st.size();
 870      oscerr << "(" <<  myProcId <<  ") I hold " << nbOfMat << " DENOMINATOR matrix(ces) : "<< std::endl;
 871      for(int i=0;i<nbOfMat;i++)
 872        {
 873          oscerr << "   - Matrix #" << i << " coming from source proc #" << _the_matrix_st_source_proc_id[i] << ": \n";
 874          const std::vector< SparseDoubleVec >& locMat=_the_deno_st[i];
 875          int j = 0;
 876          for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++)
 877            {
 878              oscerr << " Target Cell #" << j;
 879              for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 880                oscerr << " (" << (*it2).first << "," << (*it2).second << "), ";
 881              oscerr << std::endl;
 882            }
 883        }
 884      oscerr << "*********" << std::endl;
 885
 886      // Hope this will be flushed in one go:
 887      std::cerr << oscerr.str() << std::endl;
 888    }
 889 #endif