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