src/ParaMEDMEM/OverlapMapping.cxx

   1 // Copyright (C) 2007-2021  CEA/DEN, EDF R&D
   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)// && _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)
 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)
 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)
 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();
 548             valsToSend.insert(valsToSend.end(),vals->getConstPointer(),vals->getConstPointer()+nbsend[procID]);
 549           }
 550
 551       /* RECEIVE: compute number of field values to be RECEIVED
 552        *   - for all proc 'procID' in group
 553        *      * if procID == myProcID, rcv nothing
 554        *      * elif 'procID' in _proc_ids_to_recv_vector_st (computed from BB intersec)
 555        *        % if myProcID computed the job (procID, myProcID)
 556        *          => receive full set ('BB source IDs') of field data from proc #procID which has never seen the matrix
 557        *             i.e. prepare to receive the numb in _nb_of_rcv_src_ids
 558        *        % else (=we did NOT compute the job, hence procID has, and knows the matrix)
 559        *          => receive 'interp source IDs' set of field values
 560        */
 561       const std::vector< int > & _proc_ids_to_recv_vector_st = _the_matrix_st_source_proc_id;
 562       if (procID == myProcID)
 563         nbrecv[procID] = 0;
 564       else
 565         if(find(_proc_ids_to_recv_vector_st.begin(),_proc_ids_to_recv_vector_st.end(),procID)!=_proc_ids_to_recv_vector_st.end())
 566           {
 567             if(_locator.isInMyTodoList(procID, myProcID))
 568               {
 569                 map <int,mcIdType>::const_iterator isItem11 = _nb_of_rcv_src_ids.find(procID);
 570                 if (isItem11 == _nb_of_rcv_src_ids.end())
 571                   throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: RCV: unexpected end iterator in _nb_of_rcv_src_ids!");
 572                 nbrecv[procID] = (int)((*isItem11).second);
 573               }
 574             else
 575               {
 576                 map<int, vector<mcIdType> >::const_iterator isItem11 = _src_ids_zip_recv.find(procID);
 577                 if (isItem11 == _src_ids_zip_recv.end())
 578                   throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: RCV: unexpected end iterator in _src_ids_zip_recv!");
 579                 nbrecv[procID] = (int)((*isItem11).second.size()*nbOfCompo);
 580               }
 581           }
 582     }
 583   // Compute offsets in the sending/receiving array.
 584   for(int i=1;i<grpSize;i++)
 585     {
 586       nbsend2[i]=nbsend2[i-1]+nbsend[i-1];
 587       nbrecv2[i]=nbrecv2[i-1]+nbrecv[i-1];
 588     }
 589   INTERP_KERNEL::AutoPtr<double> bigArr=new double[nbrecv2[grpSize-1]+nbrecv[grpSize-1]];
 590
 591 #ifdef DEC_DEBUG
 592   stringstream scout;
 593   scout << "("  << myProcID << ") nbsend :" << nbsend[0] << "," << nbsend[1] << "," << nbsend[2] << "\n";
 594   scout << "("  << myProcID << ") nbrecv :" << nbrecv[0] << "," << nbrecv[1] << "," << nbrecv[2] << "\n";
 595   scout << "("  << myProcID << ") valsToSend: ";
 596   for (int iii=0; iii<valsToSend.size(); iii++)
 597     scout << ", " << valsToSend[iii];
 598   cout << scout.str() << "\n";
 599 #endif
 600
 601   /*
 602    * *********************** ALL-TO-ALL
 603    */
 604   commInterface.allToAllV(&valsToSend[0],nbsend,nbsend2,MPI_DOUBLE,
 605                           bigArr,nbrecv,nbrecv2,MPI_DOUBLE,*comm);
 606 #ifdef DEC_DEBUG
 607   MPI_Barrier(MPI_COMM_WORLD);
 608   scout << "("  << myProcID << ") bigArray: ";
 609     for (int iii=0; iii<nbrecv2[grpSize-1]+nbrecv[grpSize-1]; iii++)
 610       scout << ", " << bigArr[iii];
 611   cout << scout.str() << "\n";
 612 #endif
 613
 614   /*
 615    * TARGET FIELD COMPUTATION (matrix-vec computation)
 616    */
 617   fieldOutput->getArray()->fillWithZero();
 618   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfCompo];
 619
 620   // By default field value set to default value - so mark which cells are hit
 621   INTERP_KERNEL::AutoPtr<bool> hit_cells = new bool[fieldOutput->getNumberOfTuples()];
 622
 623   for(vector<int>::const_iterator itProc=_the_matrix_st_source_proc_id.begin(); itProc != _the_matrix_st_source_proc_id.end();itProc++)
 624   // For each source processor corresponding to a locally held matrix:
 625     {
 626       int srcProcID = *itProc;
 627       std::size_t id = std::distance(_the_matrix_st_source_proc_id.begin(),itProc);
 628       const vector< SparseDoubleVec >& mat =_the_matrix_st[id];
 629       const vector< SparseDoubleVec >& deno = _the_deno_st[id];
 630
 631       /*   FINAL MULTIPLICATION
 632        *      * if srcProcID == myProcID, local multiplication without any mapping
 633        *         => for all target cell ID 'tgtCellID'
 634        *           => for all src cell ID 'srcCellID' in the sparse vector
 635        *             => tgtFieldLocal[tgtCellID] += srcFieldLocal[srcCellID] * matrix[tgtCellID][srcCellID] / deno[tgtCellID][srcCellID]
 636        */
 637       if (srcProcID == myProcID)
 638         {
 639           std::size_t nbOfTrgTuples=mat.size();
 640           double * targetBase = fieldOutput->getArray()->getPointer();
 641           for(std::size_t j=0; j<nbOfTrgTuples; j++)
 642             {
 643               const SparseDoubleVec& mat1=mat[j];
 644               const SparseDoubleVec& deno1=deno[j];
 645               SparseDoubleVec::const_iterator it5=deno1.begin();
 646               const double * localSrcField = fieldInput->getArray()->getConstPointer();
 647               double * targetPt = targetBase+j*nbOfCompo;
 648               for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
 649                 {
 650                   // Apply the multiplication for all components:
 651                   double ratio = (*it3).second/(*it5).second;
 652                   transform(localSrcField+((*it3).first)*nbOfCompo,
 653                             localSrcField+((*it3).first+1)*nbOfCompo,
 654                             (double *)tmp,
 655                             bind2nd(multiplies<double>(),ratio) );
 656                   // Accumulate with current value:
 657                   transform((double *)tmp,(double *)tmp+nbOfCompo,targetPt,targetPt,plus<double>());
 658                   hit_cells[j] = true;
 659                 }
 660             }
 661         }
 662
 663       if(nbrecv[srcProcID]<=0)  // also covers the preceding 'if'
 664         continue;
 665
 666       /*      * if something was received
 667        *         %  if received matrix (=we didn't compute the job), this means that :
 668        *            1. we sent part of our targetIDs to srcProcID before, so that srcProcId can do the computation.
 669        *            2. srcProcID has sent us only the 'interp source IDs' field values
 670        *            => invert _src_ids_zip_recv -> 'revert_zip'
 671        *            => for all target cell ID 'tgtCellID'
 672        *              => mappedTgtID = _sent_trg_ids[srcProcID][tgtCellID]
 673        *              => for all src cell ID 'srcCellID' in the sparse vector
 674        *                 => idx = revert_zip[srcCellID]
 675        *                 => tgtFieldLocal[mappedTgtID] += rcvValue[srcProcID][idx] * matrix[tgtCellID][srcCellID] / deno[tgtCellID][srcCellID]
 676        */
 677       if(!_locator.isInMyTodoList(srcProcID, myProcID))
 678         {
 679           // invert _src_ids_zip_recv
 680           map<mcIdType,int> revert_zip;
 681           map<int, vector<mcIdType> >::const_iterator it11= _src_ids_zip_recv.find(srcProcID);
 682           if (it11 == _src_ids_zip_recv.end())
 683             throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: MULTIPLY: unexpected end iterator in _src_ids_zip_recv!");
 684
 685           const vector<mcIdType> & vec = (*it11).second;
 686           int newId=0;
 687           for(vector<mcIdType>::const_iterator it=vec.begin();it!=vec.end();it++,newId++)
 688             revert_zip[*it]=newId;
 689           map < int, MCAuto<DataArrayIdType> >::const_iterator isItem24 = _sent_trg_ids.find(srcProcID);
 690           if (isItem24 == _sent_trg_ids.end())
 691             throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: MULTIPLY: unexpected end iterator in _sent_trg_ids!");
 692           const DataArrayIdType *tgrIdsDA = (*isItem24).second;
 693           const mcIdType *tgrIds = tgrIdsDA->getConstPointer();
 694
 695           std::size_t nbOfTrgTuples=mat.size();
 696           double * targetBase = fieldOutput->getArray()->getPointer();
 697           for(std::size_t j=0;j<nbOfTrgTuples;j++)
 698             {
 699               const SparseDoubleVec& mat1=mat[j];
 700               const SparseDoubleVec& deno1=deno[j];
 701               SparseDoubleVec::const_iterator it5=deno1.begin();
 702               double * targetPt = targetBase+tgrIds[j]*nbOfCompo;
 703               for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
 704                 {
 705                   map<mcIdType,int>::const_iterator it4=revert_zip.find((*it3).first);
 706                   if(it4==revert_zip.end())
 707                     throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: MULTIPLY: unexpected end iterator in revert_zip!");
 708                   double ratio = (*it3).second/(*it5).second;
 709                   transform(bigArr+nbrecv2[srcProcID]+((*it4).second)*nbOfCompo,
 710                             bigArr+nbrecv2[srcProcID]+((*it4).second+1)*nbOfCompo,
 711                             (double *)tmp,
 712                             bind2nd(multiplies<double>(),ratio) );
 713                   transform((double *)tmp,(double *)tmp+nbOfCompo,targetPt,targetPt,plus<double>());
 714                   hit_cells[tgrIds[j]] = true;
 715                 }
 716             }
 717         }
 718       else
 719         /*         % else (=we computed the job and we received the 'BB source IDs' set of source field values)
 720          *            => for all target cell ID 'tgtCellID'
 721          *              => for all src cell ID 'srcCellID' in the sparse vector
 722          *                => tgtFieldLocal[tgtCellID] += rcvValue[srcProcID][srcCellID] * matrix[tgtCellID][srcCellID] / deno[tgtCellID][srcCellID]
 723          */
 724         {
 725           // Same loop as in the case srcProcID == myProcID, except that instead of working on local field data, we work on bigArr
 726           std::size_t nbOfTrgTuples=mat.size();
 727           double * targetBase = fieldOutput->getArray()->getPointer();
 728           for(std::size_t j=0;j<nbOfTrgTuples;j++)
 729             {
 730               const SparseDoubleVec& mat1=mat[j];
 731               const SparseDoubleVec& deno1=deno[j];
 732               SparseDoubleVec::const_iterator it5=deno1.begin();
 733               double * targetPt = targetBase+j*nbOfCompo;
 734               for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
 735                 {
 736                   // Apply the multiplication for all components:
 737                   double ratio = (*it3).second/(*it5).second;
 738                   transform(bigArr+nbrecv2[srcProcID]+((*it3).first)*nbOfCompo,
 739                             bigArr+nbrecv2[srcProcID]+((*it3).first+1)*nbOfCompo,
 740                             (double *)tmp,
 741                             bind2nd(multiplies<double>(),ratio));
 742                   // Accumulate with current value:
 743                   transform((double *)tmp,(double *)tmp+nbOfCompo,targetPt,targetPt,plus<double>());
 744                   hit_cells[j] = true;
 745                 }
 746             }
 747         }
 748     }
 749
 750   // Fill in default values for cells which haven't been hit:
 751   int i = 0;
 752   for(bool * hit_cells_ptr=hit_cells; i< fieldOutput->getNumberOfTuples(); hit_cells_ptr++,i++)
 753     if (!(*hit_cells_ptr))
 754       {
 755         double * targetPt=fieldOutput->getArray()->getPointer();
 756         fill(targetPt+i*nbOfCompo, targetPt+(i+1)*nbOfCompo, default_val);
 757       }
 758 }
 759
 760 /*!
 761  * This method performs a transpose multiply of 'fieldInput' and put the result into 'fieldOutput'.
 762  * 'fieldInput' is expected to be the targetfield and 'fieldOutput' the sourcefield.
 763  */
 764 void OverlapMapping::transposeMultiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput)
 765 {
 766 }
 767
 768 /*!
 769  * This method should be called immediately after _the_matrix_st has been filled with remote computed matrix
 770  * put in this proc for Matrix-Vector.
 771  * It fills _src_ids_zip_recv (see member doc)
 772  */
 773 void OverlapMapping::fillSourceIdsZipReceivedForMultiply()
 774 {
 775   /* When it is called, only the bits received from other processors (i.e. the remotely executed jobs) are in the
 776     big matrix _the_matrix_st. */
 777
 778   CommInterface commInterface=_group.getCommInterface();
 779   int myProcId=_group.myRank();
 780   std::size_t nbOfMatrixRecveived=_the_matrix_st_source_proc_id.size();
 781   for(std::size_t i=0;i<nbOfMatrixRecveived;i++)
 782     {
 783       int curSrcProcId=_the_matrix_st_source_proc_id[i];
 784       if(curSrcProcId!=myProcId)  // if =, data has been populated by addContributionST()
 785         {
 786           const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i];
 787           std::set<mcIdType> s;
 788           for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++)
 789             for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 790               s.insert((*it2).first);
 791           vector<mcIdType> vec(s.begin(),s.end());
 792           _src_ids_zip_recv[curSrcProcId] = vec;
 793         }
 794     }
 795 }
 796
 797 #ifdef DEC_DEBUG
 798  void OverlapMapping::printTheMatrix() const
 799  {
 800    CommInterface commInterface=_group.getCommInterface();
 801    const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 802    const MPI_Comm *comm=group->getComm();
 803    int grpSize=_group.size();
 804    int myProcId=_group.myRank();
 805    std::stringstream oscerr;
 806    int nbOfMat=_the_matrix_st.size();
 807    oscerr << "(" <<  myProcId <<  ") I hold " << nbOfMat << " matrix(ces) : "<< std::endl;
 808    for(int i=0;i<nbOfMat;i++)
 809      {
 810        oscerr << "   - Matrix #" << i << " coming from source proc #" << _the_matrix_st_source_proc_id[i] << ":\n ";
 811        const std::vector< SparseDoubleVec >& locMat=_the_matrix_st[i];
 812        int j = 0;
 813        for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++)
 814          {
 815            oscerr << " Target Cell #" << j;
 816            for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 817              oscerr << " (" << (*it2).first << "," << (*it2).second << "), ";
 818            oscerr << std::endl;
 819          }
 820      }
 821    oscerr << "*********" << std::endl;
 822
 823    // Hope this will be flushed in one go:
 824    std::cerr << oscerr.str() << std::endl;
 825 //   if(myProcId != 0)
 826 //     MPI_Barrier(MPI_COMM_WORLD);
 827  }
 828
 829  void OverlapMapping::printMatrixesST() const
 830   {
 831     CommInterface commInterface=_group.getCommInterface();
 832     const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 833     const MPI_Comm *comm=group->getComm();
 834     int grpSize=_group.size();
 835     int myProcId=_group.myRank();
 836     std::stringstream oscerr;
 837     int nbOfMat=_matrixes_st.size();
 838     oscerr << "(" <<  myProcId <<  ") I hold " << nbOfMat << " LOCAL matrix(ces) : "<< std::endl;
 839     for(int i=0;i<nbOfMat;i++)
 840       {
 841         oscerr << "   - Matrix #" << i << ": (source proc #" << _source_proc_id_st[i] << " / tgt proc#" << _target_proc_id_st[i] << "): \n";
 842         const std::vector< SparseDoubleVec >& locMat=_matrixes_st[i];
 843         int j = 0;
 844         for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++)
 845           {
 846             oscerr << " Target Cell #" << j;
 847             for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 848               oscerr << " (" << (*it2).first << "," << (*it2).second << "), ";
 849             oscerr << std::endl;
 850           }
 851       }
 852     oscerr << "*********" << std::endl;
 853
 854     // Hope this will be flushed in one go:
 855     std::cerr << oscerr.str() << std::endl;
 856   }
 857
 858  void OverlapMapping::printDenoMatrix() const
 859    {
 860      CommInterface commInterface=_group.getCommInterface();
 861      const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 862      const MPI_Comm *comm=group->getComm();
 863      int grpSize=_group.size();
 864      int myProcId=_group.myRank();
 865      std::stringstream oscerr;
 866      int nbOfMat=_the_deno_st.size();
 867      oscerr << "(" <<  myProcId <<  ") I hold " << nbOfMat << " DENOMINATOR matrix(ces) : "<< std::endl;
 868      for(int i=0;i<nbOfMat;i++)
 869        {
 870          oscerr << "   - Matrix #" << i << " coming from source proc #" << _the_matrix_st_source_proc_id[i] << ": \n";
 871          const std::vector< SparseDoubleVec >& locMat=_the_deno_st[i];
 872          int j = 0;
 873          for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++)
 874            {
 875              oscerr << " Target Cell #" << j;
 876              for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 877                oscerr << " (" << (*it2).first << "," << (*it2).second << "), ";
 878              oscerr << std::endl;
 879            }
 880        }
 881      oscerr << "*********" << std::endl;
 882
 883      // Hope this will be flushed in one go:
 884      std::cerr << oscerr.str() << std::endl;
 885    }
 886 #endif