src/ParaMEDMEM/OverlapMapping.cxx

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