]> SALOME platform Git repositories - modules/med.git/blob - src/ParaMEDMEM/OverlapMapping.cxx
Salome HOME
Merge from V6_main 01/04/2013
[modules/med.git] / src / ParaMEDMEM / OverlapMapping.cxx
1 // Copyright (C) 2007-2013  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 // 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 // }