src/ParaMEDMEM/OverlapMapping.cxx

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