src/ParaMEDMEM/OverlapMapping.cxx

   1 // Copyright (C) 2007-2012  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.
   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
  20 #include "OverlapMapping.hxx"
  21 #include "MPIProcessorGroup.hxx"
  22
  23 #include "MEDCouplingFieldDouble.hxx"
  24 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
  25
  26 #include "InterpKernelAutoPtr.hxx"
  27
  28 #include <numeric>
  29 #include <algorithm>
  30
  31 using namespace ParaMEDMEM;
  32
  33 OverlapMapping::OverlapMapping(const ProcessorGroup& group):_group(group)
  34 {
  35 }
  36
  37 /*!
  38  * This method keeps tracks of source ids to know in step 6 of main algorithm, which tuple ids to send away.
  39  * This method incarnates item#1 of step2 algorithm.
  40  */
  41 void OverlapMapping::keepTracksOfSourceIds(int procId, DataArrayInt *ids)
  42 {
  43   ids->incrRef();
  44   _src_ids_st2.push_back(ids);
  45   _src_proc_st2.push_back(procId);
  46 }
  47
  48 /*!
  49  * This method keeps tracks of target ids to know in step 6 of main algorithm.
  50  * This method incarnates item#0 of step2 algorithm.
  51  */
  52 void OverlapMapping::keepTracksOfTargetIds(int procId, DataArrayInt *ids)
  53 {
  54   ids->incrRef();
  55   _trg_ids_st2.push_back(ids);
  56   _trg_proc_st2.push_back(procId);
  57 }
  58
  59 /*!
  60  * This method stores from a matrix in format Target(rows)/Source(cols) for a source procId 'srcProcId' and for a target procId 'trgProcId'.
  61  * All ids (source and target) are in format of local ids.
  62  */
  63 void OverlapMapping::addContributionST(const std::vector< std::map<int,double> >& matrixST, const DataArrayInt *srcIds, int srcProcId, const DataArrayInt *trgIds, int trgProcId)
  64 {
  65   _matrixes_st.push_back(matrixST);
  66   _source_proc_id_st.push_back(srcProcId);
  67   _target_proc_id_st.push_back(trgProcId);
  68   if(srcIds)
  69     {//item#1 of step2 algorithm in proc m. Only to know in advanced nb of recv ids [ (0,1) computed on proc1 and Matrix-Vector on proc1 ]
  70       _nb_of_src_ids_proc_st2.push_back(srcIds->getNumberOfTuples());
  71       _src_ids_proc_st2.push_back(srcProcId);
  72     }
  73   else
  74     {//item#0 of step2 algorithm in proc k
  75       std::set<int> s;
  76       for(std::vector< std::map<int,double> >::const_iterator it1=matrixST.begin();it1!=matrixST.end();it1++)
  77         for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
  78           s.insert((*it2).first);
  79       _src_ids_zip_st2.resize(_src_ids_zip_st2.size()+1);
  80       _src_ids_zip_st2.back().insert(_src_ids_zip_st2.back().end(),s.begin(),s.end());
  81       _src_ids_zip_proc_st2.push_back(trgProcId);
  82     }
  83 }
  84
  85 /*!
  86  * 'procsInInteraction' gives the global view of interaction between procs.
  87  * In 'procsInInteraction' for a proc with id i, is in interaction with procs listed in procsInInteraction[i].
  88  *
  89  * This method is in charge to send matrixes in AlltoAll mode.
  90  * After the call of this method 'this' contains the matrixST for all source elements of the current proc
  91  */
  92 void OverlapMapping::prepare(const std::vector< std::vector<int> >& procsInInteraction, int nbOfTrgElems)
  93 {
  94   CommInterface commInterface=_group.getCommInterface();
  95   const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
  96   const MPI_Comm *comm=group->getComm();
  97   int grpSize=_group.size();
  98   INTERP_KERNEL::AutoPtr<int> nbsend=new int[grpSize];
  99   INTERP_KERNEL::AutoPtr<int> nbsend2=new int[grpSize];
 100   INTERP_KERNEL::AutoPtr<int> nbsend3=new int[grpSize];
 101   std::fill<int *>(nbsend,nbsend+grpSize,0);
 102   int myProcId=_group.myRank();
 103   _proc_ids_to_recv_vector_st.clear();
 104   int curProc=0;
 105   for(std::vector< std::vector<int> >::const_iterator it1=procsInInteraction.begin();it1!=procsInInteraction.end();it1++,curProc++)
 106     if(std::find((*it1).begin(),(*it1).end(),myProcId)!=(*it1).end())
 107       _proc_ids_to_recv_vector_st.push_back(curProc);
 108   _proc_ids_to_send_vector_st=procsInInteraction[myProcId];
 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]]=_matrixes_st[i].size();
 112   INTERP_KERNEL::AutoPtr<int> nbrecv=new int[grpSize];
 113   commInterface.allToAll(nbsend,1,MPI_INT,nbrecv,1,MPI_INT,*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   int *tmp=0;
 120   serializeMatrixStep0ST(nbrecv,
 121                          tmp,nbsend2,nbsend3,
 122                          nbrecv1,nbrecv2);
 123   INTERP_KERNEL::AutoPtr<int> bigArr=tmp;
 124   INTERP_KERNEL::AutoPtr<int> bigArrRecv=new int[nbrecv2[grpSize-1]+nbrecv1[grpSize-1]];
 125   commInterface.allToAllV(bigArr,nbsend2,nbsend3,MPI_INT,
 126                           bigArrRecv,nbrecv1,nbrecv2,MPI_INT,
 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   int lgthOfArr=serializeMatrixStep1ST(nbrecv,bigArrRecv,nbrecv1,nbrecv2,
 134                                        tmp,tmp2,
 135                                        nbsend2,nbsend3,nbrecv3,nbrecv4);
 136   INTERP_KERNEL::AutoPtr<int> bigArr2=tmp;
 137   INTERP_KERNEL::AutoPtr<double> bigArrD2=tmp2;
 138   INTERP_KERNEL::AutoPtr<int> bigArrRecv2=new int[lgthOfArr];
 139   INTERP_KERNEL::AutoPtr<double> bigArrDRecv2=new double[lgthOfArr];
 140   commInterface.allToAllV(bigArr2,nbsend2,nbsend3,MPI_INT,
 141                           bigArrRecv2,nbrecv3,nbrecv4,MPI_INT,
 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   //updating _src_ids_zip_st2 and _src_ids_zip_st2 with received matrix.
 150   updateZipSourceIdsForFuture();
 151   //finish to fill _the_matrix_st with already in place matrix in _matrixes_st
 152   finishToFillFinalMatrixST();
 153   //printTheMatrix();
 154 }
 155
 156 /*!
 157  * Compute denominators.
 158  */
 159 void OverlapMapping::computeDenoGlobConstraint()
 160 {
 161   _the_deno_st.clear();
 162   std::size_t sz1=_the_matrix_st.size();
 163   _the_deno_st.resize(sz1);
 164   for(std::size_t i=0;i<sz1;i++)
 165     {
 166       std::size_t sz2=_the_matrix_st[i].size();
 167       _the_deno_st[i].resize(sz2);
 168       for(std::size_t j=0;j<sz2;j++)
 169         {
 170           double sum=0;
 171           std::map<int,double>& mToFill=_the_deno_st[i][j];
 172           const std::map<int,double>& m=_the_matrix_st[i][j];
 173           for(std::map<int,double>::const_iterator it=m.begin();it!=m.end();it++)
 174             sum+=(*it).second;
 175           for(std::map<int,double>::const_iterator it=m.begin();it!=m.end();it++)
 176             mToFill[(*it).first]=sum;
 177         }
 178     }
 179 }
 180
 181 /*!
 182  * Compute denominators.
 183  */
 184 void OverlapMapping::computeDenoConservativeVolumic(int nbOfTuplesTrg)
 185 {
 186   CommInterface commInterface=_group.getCommInterface();
 187   int myProcId=_group.myRank();
 188   //
 189   _the_deno_st.clear();
 190   std::size_t sz1=_the_matrix_st.size();
 191   _the_deno_st.resize(sz1);
 192   std::vector<double> deno(nbOfTuplesTrg);
 193   for(std::size_t i=0;i<sz1;i++)
 194     {
 195       const std::vector< std::map<int,double> >& mat=_the_matrix_st[i];
 196       int curSrcId=_the_matrix_st_source_proc_id[i];
 197       std::vector<int>::iterator isItem1=std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),curSrcId);
 198       int rowId=0;
 199       if(isItem1==_trg_proc_st2.end() || curSrcId==myProcId)//item1 of step2 main algo. Simple, because rowId of mat are directly target ids.
 200         {
 201           for(std::vector< std::map<int,double> >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
 202             for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 203               deno[rowId]+=(*it2).second;
 204         }
 205       else
 206         {//item0 of step2 main algo. More complicated.
 207           std::vector<int>::iterator fnd=isItem1;//std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),curSrcId);
 208           int locId=std::distance(_trg_proc_st2.begin(),fnd);
 209           const DataArrayInt *trgIds=_trg_ids_st2[locId];
 210           const int *trgIds2=trgIds->getConstPointer();
 211           for(std::vector< std::map<int,double> >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
 212             for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 213               deno[trgIds2[rowId]]+=(*it2).second;
 214         }
 215     }
 216   //
 217   for(std::size_t i=0;i<sz1;i++)
 218     {
 219       int rowId=0;
 220       const std::vector< std::map<int,double> >& mat=_the_matrix_st[i];
 221       int curSrcId=_the_matrix_st_source_proc_id[i];
 222       std::vector<int>::iterator isItem1=std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),curSrcId);
 223       std::vector< std::map<int,double> >& denoM=_the_deno_st[i];
 224       denoM.resize(mat.size());
 225       if(isItem1==_trg_proc_st2.end() || curSrcId==myProcId)//item1 of step2 main algo. Simple, because rowId of mat are directly target ids.
 226         {
 227           int rowId=0;
 228           for(std::vector< std::map<int,double> >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
 229             for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 230               denoM[rowId][(*it2).first]=deno[rowId];
 231         }
 232       else
 233         {
 234           std::vector<int>::iterator fnd=isItem1;
 235           int locId=std::distance(_trg_proc_st2.begin(),fnd);
 236           const DataArrayInt *trgIds=_trg_ids_st2[locId];
 237           const int *trgIds2=trgIds->getConstPointer();
 238           for(std::vector< std::map<int,double> >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
 239             for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 240               denoM[rowId][(*it2).first]=deno[trgIds2[rowId]];
 241         }
 242     }
 243 }
 244
 245 /*!
 246  * This method performs step #0/3 in serialization process.
 247  * \param count tells specifies nb of elems to send to corresponding proc id. size equal to _group.size().
 248  * \param offsets tells for a proc i where to start serialize#0 matrix. size equal to _group.size().
 249  * \param nbOfElemsSrc of size _group.size(). Comes from previous all2all call. tells how many srcIds per proc contains matrix for current proc.
 250  */
 251 void OverlapMapping::serializeMatrixStep0ST(const int *nbOfElemsSrc, int *&bigArr, int *count, int *offsets,
 252                                             int *countForRecv, int *offsetsForRecv) const
 253 {
 254   int grpSize=_group.size();
 255   std::fill<int *>(count,count+grpSize,0);
 256   int szz=0;
 257   int myProcId=_group.myRank();
 258   for(std::size_t i=0;i<_matrixes_st.size();i++)
 259     {
 260       if(_source_proc_id_st[i]==myProcId)// && _target_proc_id_st[i]!=myProcId
 261         {
 262           count[_target_proc_id_st[i]]=_matrixes_st[i].size()+1;
 263           szz+=_matrixes_st[i].size()+1;
 264         }
 265     }
 266   bigArr=new int[szz];
 267   offsets[0]=0;
 268   for(int i=1;i<grpSize;i++)
 269     offsets[i]=offsets[i-1]+count[i-1];
 270   for(std::size_t i=0;i<_matrixes_st.size();i++)
 271     {
 272       if(_source_proc_id_st[i]==myProcId)
 273         {
 274           int start=offsets[_target_proc_id_st[i]];
 275           int *work=bigArr+start;
 276           *work=0;
 277           const std::vector< std::map<int,double> >& mat=_matrixes_st[i];
 278           for(std::vector< std::map<int,double> >::const_iterator it=mat.begin();it!=mat.end();it++,work++)
 279             work[1]=work[0]+(*it).size();
 280         }
 281     }
 282   //
 283   offsetsForRecv[0]=0;
 284   for(int i=0;i<grpSize;i++)
 285     {
 286       if(nbOfElemsSrc[i]>0)
 287         countForRecv[i]=nbOfElemsSrc[i]+1;
 288       else
 289         countForRecv[i]=0;
 290       if(i>0)
 291         offsetsForRecv[i]=offsetsForRecv[i-1]+countForRecv[i-1];
 292     }
 293 }
 294
 295 /*!
 296  * This method performs step#1 and step#2/3. It returns the size of expected array to get allToAllV.
 297  */
 298 int OverlapMapping::serializeMatrixStep1ST(const int *nbOfElemsSrc, const int *recvStep0, const int *countStep0, const int *offsStep0,
 299                                            int *&bigArrI, double *&bigArrD, int *count, int *offsets,
 300                                            int *countForRecv, int *offsForRecv) const
 301 {
 302   int grpSize=_group.size();
 303   int myProcId=_group.myRank();
 304   offsForRecv[0]=0;
 305   int szz=0;
 306   for(int i=0;i<grpSize;i++)
 307     {
 308       if(nbOfElemsSrc[i]!=0)
 309         countForRecv[i]=recvStep0[offsStep0[i]+nbOfElemsSrc[i]];
 310       else
 311         countForRecv[i]=0;
 312       szz+=countForRecv[i];
 313       if(i>0)
 314         offsForRecv[i]=offsForRecv[i-1]+countForRecv[i-1];
 315     }
 316   //
 317   std::fill(count,count+grpSize,0);
 318   offsets[0]=0;
 319   int fullLgth=0;
 320   for(std::size_t i=0;i<_matrixes_st.size();i++)
 321     {
 322       if(_source_proc_id_st[i]==myProcId)
 323         {
 324           const std::vector< std::map<int,double> >& mat=_matrixes_st[i];
 325           int lgthToSend=0;
 326           for(std::vector< std::map<int,double> >::const_iterator it=mat.begin();it!=mat.end();it++)
 327             lgthToSend+=(*it).size();
 328           count[_target_proc_id_st[i]]=lgthToSend;
 329           fullLgth+=lgthToSend;
 330         }
 331     }
 332   for(int i=1;i<grpSize;i++)
 333     offsets[i]=offsets[i-1]+count[i-1];
 334   //
 335   bigArrI=new int[fullLgth];
 336   bigArrD=new double[fullLgth];
 337   // feeding arrays
 338   fullLgth=0;
 339   for(std::size_t i=0;i<_matrixes_st.size();i++)
 340     {
 341       if(_source_proc_id_st[i]==myProcId)
 342         {
 343           const std::vector< std::map<int,double> >& mat=_matrixes_st[i];
 344           for(std::vector< std::map<int,double> >::const_iterator it1=mat.begin();it1!=mat.end();it1++)
 345             {
 346               int j=0;
 347               for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++,j++)
 348                 {
 349                   bigArrI[fullLgth+j]=(*it2).first;
 350                   bigArrD[fullLgth+j]=(*it2).second;
 351                 }
 352               fullLgth+=(*it1).size();
 353             }
 354         }
 355     }
 356   return szz;
 357 }
 358
 359 /*!
 360  * This is the last step after all2Alls for matrix exchange.
 361  * _the_matrix_st is the final matrix :
 362  *      - The first entry is srcId in current proc.
 363  *      - The second is the pseudo id of source proc (correspondance with true id is in attribute _the_matrix_st_source_proc_id and _the_matrix_st_source_ids)
 364  *      - the third is the srcId in the pseudo source proc
 365  */
 366 void OverlapMapping::unserializationST(int nbOfTrgElems,
 367                                        const int *nbOfElemsSrcPerProc,//first all2all
 368                                        const int *bigArrRecv, const int *bigArrRecvCounts, const int *bigArrRecvOffs,//2nd all2all
 369                                        const int *bigArrRecv2, const double *bigArrDRecv2, const int *bigArrRecv2Count, const int *bigArrRecv2Offs)//3rd and 4th all2alls
 370 {
 371   _the_matrix_st.clear();
 372   _the_matrix_st_source_proc_id.clear();
 373   //
 374   int grpSize=_group.size();
 375   for(int i=0;i<grpSize;i++)
 376     if(nbOfElemsSrcPerProc[i]!=0)
 377       _the_matrix_st_source_proc_id.push_back(i);
 378   int nbOfPseudoProcs=_the_matrix_st_source_proc_id.size();//_the_matrix_st_target_proc_id.size() contains number of matrix fetched remotely whose sourceProcId==myProcId
 379   _the_matrix_st.resize(nbOfPseudoProcs);
 380   //
 381   int j=0;
 382   for(int i=0;i<grpSize;i++)
 383     if(nbOfElemsSrcPerProc[i]!=0)
 384       {
 385         _the_matrix_st[j].resize(nbOfElemsSrcPerProc[i]);
 386         for(int k=0;k<nbOfElemsSrcPerProc[i];k++)
 387           {
 388             int offs=bigArrRecv[bigArrRecvOffs[i]+k];
 389             int lgthOfMap=bigArrRecv[bigArrRecvOffs[i]+k+1]-offs;
 390             for(int l=0;l<lgthOfMap;l++)
 391               _the_matrix_st[j][k][bigArrRecv2[bigArrRecv2Offs[i]+offs+l]]=bigArrDRecv2[bigArrRecv2Offs[i]+offs+l];
 392           }
 393         j++;
 394       }
 395 }
 396
 397 /*!
 398  * This method should be called when all remote matrix with sourceProcId==thisProcId have been retrieved and are in 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id'
 399  * and 'this->_the_matrix_st_target_ids'.
 400  * This method finish the job of filling 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id' by putting candidates in 'this->_matrixes_st' into them.
 401  */
 402 void OverlapMapping::finishToFillFinalMatrixST()
 403 {
 404   int myProcId=_group.myRank();
 405   int sz=_matrixes_st.size();
 406   int nbOfEntryToAdd=0;
 407   for(int i=0;i<sz;i++)
 408     if(_source_proc_id_st[i]!=myProcId)
 409       nbOfEntryToAdd++;
 410   if(nbOfEntryToAdd==0)
 411     return ;
 412   int oldNbOfEntry=_the_matrix_st.size();
 413   int newNbOfEntry=oldNbOfEntry+nbOfEntryToAdd;
 414   _the_matrix_st.resize(newNbOfEntry);
 415   int j=oldNbOfEntry;
 416   for(int i=0;i<sz;i++)
 417     if(_source_proc_id_st[i]!=myProcId)
 418       {
 419         const std::vector<std::map<int,double> >& mat=_matrixes_st[i];
 420         _the_matrix_st[j]=mat;
 421         _the_matrix_st_source_proc_id.push_back(_source_proc_id_st[i]);
 422         j++;
 423       }
 424   _matrixes_st.clear();
 425 }
 426
 427 /*!
 428  * This method performs the operation of target ids broadcasting.
 429  */
 430 void OverlapMapping::prepareIdsToSendST()
 431 {
 432   CommInterface commInterface=_group.getCommInterface();
 433   const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 434   const MPI_Comm *comm=group->getComm();
 435   int grpSize=_group.size();
 436   _source_ids_to_send_st.clear();
 437   _source_ids_to_send_st.resize(grpSize);
 438   INTERP_KERNEL::AutoPtr<int> nbsend=new int[grpSize];
 439   std::fill<int *>(nbsend,nbsend+grpSize,0);
 440   for(std::size_t i=0;i<_the_matrix_st_source_proc_id.size();i++)
 441     nbsend[_the_matrix_st_source_proc_id[i]]=_the_matrix_st_source_ids[i].size();
 442   INTERP_KERNEL::AutoPtr<int> nbrecv=new int[grpSize];
 443   commInterface.allToAll(nbsend,1,MPI_INT,nbrecv,1,MPI_INT,*comm);
 444   //
 445   INTERP_KERNEL::AutoPtr<int> nbsend2=new int[grpSize];
 446   std::copy((int *)nbsend,((int *)nbsend)+grpSize,(int *)nbsend2);
 447   INTERP_KERNEL::AutoPtr<int> nbsend3=new int[grpSize];
 448   nbsend3[0]=0;
 449   for(int i=1;i<grpSize;i++)
 450     nbsend3[i]=nbsend3[i-1]+nbsend2[i-1];
 451   int sendSz=nbsend3[grpSize-1]+nbsend2[grpSize-1];
 452   INTERP_KERNEL::AutoPtr<int> bigDataSend=new int[sendSz];
 453   for(std::size_t i=0;i<_the_matrix_st_source_proc_id.size();i++)
 454     {
 455       int offset=nbsend3[_the_matrix_st_source_proc_id[i]];
 456       std::copy(_the_matrix_st_source_ids[i].begin(),_the_matrix_st_source_ids[i].end(),((int *)nbsend3)+offset);
 457     }
 458   INTERP_KERNEL::AutoPtr<int> nbrecv2=new int[grpSize];
 459   INTERP_KERNEL::AutoPtr<int> nbrecv3=new int[grpSize];
 460   std::copy((int *)nbrecv,((int *)nbrecv)+grpSize,(int *)nbrecv2);
 461   nbrecv3[0]=0;
 462   for(int i=1;i<grpSize;i++)
 463     nbrecv3[i]=nbrecv3[i-1]+nbrecv2[i-1];
 464   int recvSz=nbrecv3[grpSize-1]+nbrecv2[grpSize-1];
 465   INTERP_KERNEL::AutoPtr<int> bigDataRecv=new int[recvSz];
 466   //
 467   commInterface.allToAllV(bigDataSend,nbsend2,nbsend3,MPI_INT,
 468                           bigDataRecv,nbrecv2,nbrecv3,MPI_INT,
 469                           *comm);
 470   for(int i=0;i<grpSize;i++)
 471     {
 472       if(nbrecv2[i]>0)
 473         {
 474           _source_ids_to_send_st[i].insert(_source_ids_to_send_st[i].end(),((int *)bigDataRecv)+nbrecv3[i],((int *)bigDataRecv)+nbrecv3[i]+nbrecv2[i]);
 475         }
 476     }
 477 }
 478
 479 /*!
 480  * This method performs a transpose multiply of 'fieldInput' and put the result into 'fieldOutput'.
 481  * 'fieldInput' is expected to be the sourcefield and 'fieldOutput' the targetfield.
 482  */
 483 void OverlapMapping::multiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput) const
 484 {
 485   int nbOfCompo=fieldInput->getNumberOfComponents();//to improve same number of components to test
 486   CommInterface commInterface=_group.getCommInterface();
 487   const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 488   const MPI_Comm *comm=group->getComm();
 489   int grpSize=_group.size();
 490   int myProcId=_group.myRank();
 491   //
 492   INTERP_KERNEL::AutoPtr<int> nbsend=new int[grpSize];
 493   INTERP_KERNEL::AutoPtr<int> nbsend2=new int[grpSize];
 494   INTERP_KERNEL::AutoPtr<int> nbrecv=new int[grpSize];
 495   INTERP_KERNEL::AutoPtr<int> nbrecv2=new int[grpSize];
 496   std::fill<int *>(nbsend,nbsend+grpSize,0);
 497   std::fill<int *>(nbrecv,nbrecv+grpSize,0);
 498   nbsend2[0]=0;
 499   nbrecv2[0]=0;
 500   std::vector<double> valsToSend;
 501   for(int i=0;i<grpSize;i++)
 502     {
 503       if(std::find(_proc_ids_to_send_vector_st.begin(),_proc_ids_to_send_vector_st.end(),i)!=_proc_ids_to_send_vector_st.end())
 504         {
 505           std::vector<int>::const_iterator isItem1=std::find(_src_proc_st2.begin(),_src_proc_st2.end(),i);
 506           MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> vals;
 507           if(isItem1!=_src_proc_st2.end())//item1 of step2 main algo
 508             {
 509               int id=std::distance(_src_proc_st2.begin(),isItem1);
 510               vals=fieldInput->getArray()->selectByTupleId(_src_ids_st2[id]->getConstPointer(),_src_ids_st2[id]->getConstPointer()+_src_ids_st2[id]->getNumberOfTuples());
 511             }
 512           else
 513             {//item0 of step2 main algo
 514               int id=std::distance(_src_ids_zip_proc_st2.begin(),std::find(_src_ids_zip_proc_st2.begin(),_src_ids_zip_proc_st2.end(),i));
 515               vals=fieldInput->getArray()->selectByTupleId(&(_src_ids_zip_st2[id])[0],&(_src_ids_zip_st2[id])[0]+_src_ids_zip_st2[id].size());
 516             }
 517           nbsend[i]=vals->getNbOfElems();
 518           valsToSend.insert(valsToSend.end(),vals->getConstPointer(),vals->getConstPointer()+nbsend[i]);
 519         }
 520       if(std::find(_proc_ids_to_recv_vector_st.begin(),_proc_ids_to_recv_vector_st.end(),i)!=_proc_ids_to_recv_vector_st.end())
 521         {
 522           std::vector<int>::const_iterator isItem0=std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),i);
 523           if(isItem0==_trg_proc_st2.end())//item1 of step2 main algo [ (0,1) computed on proc1 and Matrix-Vector on proc1 ]
 524             {
 525               std::vector<int>::const_iterator it1=std::find(_src_ids_proc_st2.begin(),_src_ids_proc_st2.end(),i);
 526               if(it1!=_src_ids_proc_st2.end())
 527                 {
 528                   int id=std::distance(_src_ids_proc_st2.begin(),it1);
 529                   nbrecv[i]=_nb_of_src_ids_proc_st2[id]*nbOfCompo;
 530                 }
 531               else if(i==myProcId)
 532                 {
 533                   nbrecv[i]=fieldInput->getNumberOfTuplesExpected()*nbOfCompo;
 534                 }
 535               else
 536                 throw INTERP_KERNEL::Exception("Plouff ! send email to anthony.geay@cea.fr ! ");
 537             }
 538           else
 539             {//item0 of step2 main algo [ (2,1) computed on proc2 but Matrix-Vector on proc1 ] [(1,0) computed on proc1 but Matrix-Vector on proc0]
 540               int id=std::distance(_src_ids_zip_proc_st2.begin(),std::find(_src_ids_zip_proc_st2.begin(),_src_ids_zip_proc_st2.end(),i));
 541               nbrecv[i]=_src_ids_zip_st2[id].size()*nbOfCompo;
 542             }
 543         }
 544     }
 545   for(int i=1;i<grpSize;i++)
 546     {
 547       nbsend2[i]=nbsend2[i-1]+nbsend[i-1];
 548       nbrecv2[i]=nbrecv2[i-1]+nbrecv[i-1];
 549     }
 550   INTERP_KERNEL::AutoPtr<double> bigArr=new double[nbrecv2[grpSize-1]+nbrecv[grpSize-1]];
 551   commInterface.allToAllV(&valsToSend[0],nbsend,nbsend2,MPI_DOUBLE,
 552                           bigArr,nbrecv,nbrecv2,MPI_DOUBLE,*comm);
 553   fieldOutput->getArray()->fillWithZero();
 554   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfCompo];
 555   for(int i=0;i<grpSize;i++)
 556     {
 557       if(nbrecv[i]>0)
 558         {
 559           double *pt=fieldOutput->getArray()->getPointer();
 560           std::vector<int>::const_iterator it=std::find(_the_matrix_st_source_proc_id.begin(),_the_matrix_st_source_proc_id.end(),i);
 561           if(it==_the_matrix_st_source_proc_id.end())
 562             throw INTERP_KERNEL::Exception("Big problem !");
 563           int id=std::distance(_the_matrix_st_source_proc_id.begin(),it);
 564           const std::vector< std::map<int,double> >& mat=_the_matrix_st[id];
 565           const std::vector< std::map<int,double> >& deno=_the_deno_st[id];
 566           std::vector<int>::const_iterator isItem0=std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),i);
 567           if(isItem0==_trg_proc_st2.end())//item1 of step2 main algo [ (0,1) computed on proc1 and Matrix-Vector on proc1 ]
 568             {
 569               int nbOfTrgTuples=mat.size();
 570               for(int j=0;j<nbOfTrgTuples;j++,pt+=nbOfCompo)
 571                 {
 572                   const std::map<int,double>& mat1=mat[j];
 573                   const std::map<int,double>& deno1=deno[j];
 574                   std::map<int,double>::const_iterator it4=deno1.begin();
 575                   for(std::map<int,double>::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it4++)
 576                     {
 577                       std::transform(bigArr+nbrecv2[i]+((*it3).first)*nbOfCompo,bigArr+nbrecv2[i]+((*it3).first+1)*(nbOfCompo),(double *)tmp,std::bind2nd(std::multiplies<double>(),(*it3).second/(*it4).second));
 578                       std::transform((double *)tmp,(double *)tmp+nbOfCompo,pt,pt,std::plus<double>());
 579                     }
 580                 }
 581             }
 582           else
 583             {//item0 of step2 main algo [ (2,1) computed on proc2 but Matrix-Vector on proc1 ]
 584               double *pt=fieldOutput->getArray()->getPointer();
 585               std::map<int,int> zipCor;
 586               int id=std::distance(_src_ids_zip_proc_st2.begin(),std::find(_src_ids_zip_proc_st2.begin(),_src_ids_zip_proc_st2.end(),i));
 587               const std::vector<int> zipIds=_src_ids_zip_st2[id];
 588               int newId=0;
 589               for(std::vector<int>::const_iterator it=zipIds.begin();it!=zipIds.end();it++,newId++)
 590                 zipCor[*it]=newId;
 591               int id2=std::distance(_trg_proc_st2.begin(),std::find(_trg_proc_st2.begin(),_trg_proc_st2.end(),i));
 592               const DataArrayInt *tgrIds=_trg_ids_st2[id2];
 593               const int *tgrIds2=tgrIds->getConstPointer();
 594               int nbOfTrgTuples=mat.size();
 595               for(int j=0;j<nbOfTrgTuples;j++)
 596                 {
 597                   const std::map<int,double>& mat1=mat[j];
 598                   const std::map<int,double>& deno1=deno[j];
 599                   std::map<int,double>::const_iterator it5=deno1.begin();
 600                   for(std::map<int,double>::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
 601                     {
 602                       std::map<int,int>::const_iterator it4=zipCor.find((*it3).first);
 603                       if(it4==zipCor.end())
 604                         throw INTERP_KERNEL::Exception("Hmmmmm send e mail to anthony.geay@cea.fr !");
 605                       std::transform(bigArr+nbrecv2[i]+((*it4).second)*nbOfCompo,bigArr+nbrecv2[i]+((*it4).second+1)*(nbOfCompo),(double *)tmp,std::bind2nd(std::multiplies<double>(),(*it3).second/(*it5).second));
 606                       std::transform((double *)tmp,(double *)tmp+nbOfCompo,pt+tgrIds2[j]*nbOfCompo,pt+tgrIds2[j]*nbOfCompo,std::plus<double>());
 607                     }
 608                 }
 609             }
 610         }
 611     }
 612 }
 613
 614 /*!
 615  * This method performs a transpose multiply of 'fieldInput' and put the result into 'fieldOutput'.
 616  * 'fieldInput' is expected to be the targetfield and 'fieldOutput' the sourcefield.
 617  */
 618 void OverlapMapping::transposeMultiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput)
 619 {
 620 }
 621
 622 /*!
 623  * This method should be called immediately after _the_matrix_st has been filled with remote computed matrix put in this proc for Matrix-Vector.
 624  * This method computes for these matrix the minimal set of source ids corresponding to the source proc id.
 625  */
 626 void OverlapMapping::updateZipSourceIdsForFuture()
 627 {
 628   CommInterface commInterface=_group.getCommInterface();
 629   int myProcId=_group.myRank();
 630   int nbOfMatrixRecveived=_the_matrix_st_source_proc_id.size();
 631   for(int i=0;i<nbOfMatrixRecveived;i++)
 632     {
 633       int curSrcProcId=_the_matrix_st_source_proc_id[i];
 634       if(curSrcProcId!=myProcId)
 635         {
 636           const std::vector< std::map<int,double> >& mat=_the_matrix_st[i];
 637           _src_ids_zip_proc_st2.push_back(curSrcProcId);
 638           _src_ids_zip_st2.resize(_src_ids_zip_st2.size()+1);
 639           std::set<int> s;
 640           for(std::vector< std::map<int,double> >::const_iterator it1=mat.begin();it1!=mat.end();it1++)
 641             for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 642               s.insert((*it2).first);
 643           _src_ids_zip_st2.back().insert(_src_ids_zip_st2.back().end(),s.begin(),s.end());
 644         }
 645     }
 646 }
 647
 648 // #include <iostream>
 649
 650 // void OverlapMapping::printTheMatrix() const
 651 // {
 652 //   CommInterface commInterface=_group.getCommInterface();
 653 //   const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
 654 //   const MPI_Comm *comm=group->getComm();
 655 //   int grpSize=_group.size();
 656 //   int myProcId=_group.myRank();
 657 //   std::cerr << "I am proc #" << myProcId << std::endl;
 658 //   int nbOfMat=_the_matrix_st.size();
 659 //   std::cerr << "I do manage " << nbOfMat << "matrix : "<< std::endl;
 660 //   for(int i=0;i<nbOfMat;i++)
 661 //     {
 662 //       std::cerr << "   - Matrix #" << i << " on source proc #" << _the_matrix_st_source_proc_id[i];
 663 //       const std::vector< std::map<int,double> >& locMat=_the_matrix_st[i];
 664 //       for(std::vector< std::map<int,double> >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++)
 665 //         {
 666 //           for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 667 //             std::cerr << "(" << (*it2).first << "," << (*it2).second << "), ";
 668 //           std::cerr << std::endl;
 669 //         }
 670 //     }
 671 //   std::cerr << "*********" << std::endl;
 672 // }