1 // Copyright (C) 2007-2024 CEA, EDF
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "OverlapMapping.hxx"
22 #include "MPIProcessorGroup.hxx"
24 #include "MEDCouplingFieldDouble.hxx"
27 #include "InterpKernelAutoPtr.hxx"
32 using namespace MEDCoupling;
34 OverlapMapping::OverlapMapping(const ProcessorGroup& group, const OverlapElementLocator & loc):
35 _group(group),_locator(loc)
40 * Keeps the link between a given a proc holding source mesh data, and the corresponding cell IDs.
42 void OverlapMapping::keepTracksOfSourceIds(int procId, DataArrayIdType *ids)
45 _sent_src_ids[procId] = ids;
49 * Same as keepTracksOfSourceIds() but for target mesh data.
51 void OverlapMapping::keepTracksOfTargetIds(int procId, DataArrayIdType *ids)
54 _sent_trg_ids[procId] = ids;
58 * This method stores in the local members the contribution coming from a matrix in format
59 * Target(rows)/Source(cols) for a source procId 'srcProcId' and for a target procId 'trgProcId'.
60 * All IDs received here (source and target) are in the format of local IDs.
62 * @param srcIds is null if the source mesh is on the local proc
63 * @param trgIds is null if the source mesh is on the local proc
65 * One of the 2 is necessarily null (the two can be null together)
67 void OverlapMapping::addContributionST(const std::vector< SparseDoubleVec >& matrixST, const DataArrayIdType *srcIds, int srcProcId, const DataArrayIdType *trgIds, int trgProcId)
69 _matrixes_st.push_back(matrixST);
70 _source_proc_id_st.push_back(srcProcId);
71 _target_proc_id_st.push_back(trgProcId);
72 if(srcIds) // source mesh part is remote <=> srcProcId != myRank
73 _nb_of_rcv_src_ids[srcProcId] = srcIds->getNumberOfTuples();
74 else // source mesh part is local
77 // For all source IDs (=col indices) in the sparse matrix:
78 for(std::vector< SparseDoubleVec >::const_iterator it1=matrixST.begin();it1!=matrixST.end();it1++)
79 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
80 s.insert((*it2).first);
81 vector<mcIdType> v(s.begin(), s.end()); // turn set into vector
82 _src_ids_zip_comp[trgProcId] = v;
87 * This method is in charge to send matrices in AlltoAll mode.
89 * 'procsToSendField' gives the list of procs field data has to be sent to.
90 * See OverlapElementLocator::computeBoundingBoxesAndTodoList()
92 * After the call of this method, 'this' contains the matrixST for all source cells of the current proc
94 void OverlapMapping::prepare(const std::vector< int >& procsToSendField, mcIdType nbOfTrgElems)
100 CommInterface commInterface=_group.getCommInterface();
101 const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
102 const MPI_Comm *comm=group->getComm();
103 std::size_t grpSize=_group.size();
104 INTERP_KERNEL::AutoPtr<mcIdType> nbsend=new mcIdType[grpSize];
105 INTERP_KERNEL::AutoPtr<int> nbsend2=new int[grpSize];
106 INTERP_KERNEL::AutoPtr<int> nbsend3=new int[grpSize];
107 std::fill<mcIdType *>(nbsend,nbsend+grpSize,0);
108 int myProcId=_group.myRank();
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]]=(int)_matrixes_st[i].size();
112 INTERP_KERNEL::AutoPtr<mcIdType> nbrecv=new mcIdType[grpSize];
113 commInterface.allToAll(nbsend,1,MPI_ID_TYPE,nbrecv,1,MPI_ID_TYPE,*comm);
115 //first exchanging offsets+ids_source
116 INTERP_KERNEL::AutoPtr<int> nbrecv1=new int[grpSize];
117 INTERP_KERNEL::AutoPtr<int> nbrecv2=new int[grpSize];
120 serializeMatrixStep0ST(nbrecv,
123 INTERP_KERNEL::AutoPtr<mcIdType> bigArr=tmp;
124 INTERP_KERNEL::AutoPtr<mcIdType> bigArrRecv=new mcIdType[nbrecv2[grpSize-1]+nbrecv1[grpSize-1]];
125 commInterface.allToAllV(bigArr,nbsend2,nbsend3,MPI_ID_TYPE,
126 bigArrRecv,nbrecv1,nbrecv2,MPI_ID_TYPE,
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];
133 mcIdType lgthOfArr=serializeMatrixStep1ST(nbrecv,bigArrRecv,nbrecv1,nbrecv2,
135 nbsend2,nbsend3,nbrecv3,nbrecv4);
136 INTERP_KERNEL::AutoPtr<mcIdType> bigArr2=tmp;
137 INTERP_KERNEL::AutoPtr<double> bigArrD2=tmp2;
138 INTERP_KERNEL::AutoPtr<mcIdType> bigArrRecv2=new mcIdType[lgthOfArr];
139 INTERP_KERNEL::AutoPtr<double> bigArrDRecv2=new double[lgthOfArr];
140 commInterface.allToAllV(bigArr2,nbsend2,nbsend3,MPI_ID_TYPE,
141 bigArrRecv2,nbrecv3,nbrecv4,MPI_ID_TYPE,
143 commInterface.allToAllV(bigArrD2,nbsend2,nbsend3,MPI_DOUBLE,
144 bigArrDRecv2,nbrecv3,nbrecv4,MPI_DOUBLE,
147 unserializationST(nbOfTrgElems,nbrecv,bigArrRecv,nbrecv1,nbrecv2,
148 bigArrRecv2,bigArrDRecv2,nbrecv3,nbrecv4);
150 //finish to fill _the_matrix_st with already in place matrix in _matrixes_st (local computation)
151 finishToFillFinalMatrixST();
153 //updating _src_ids_zip_st2 and _src_ids_zip_st2 with received matrix.
154 fillSourceIdsZipReceivedForMultiply();
155 // Prepare proc list for future field data exchange (multiply()):
156 _proc_ids_to_send_vector_st = procsToSendField;
157 // Make some space on local proc:
158 _matrixes_st.clear();
166 // * Compute denominators for ExtensiveConservation interp.
167 // * TO BE REVISED: needs another communication since some bits are held non locally
169 //void OverlapMapping::computeDenoGlobConstraint()
171 // _the_deno_st.clear();
172 // std::size_t sz1=_the_matrix_st.size();
173 // _the_deno_st.resize(sz1);
174 // for(std::size_t i=0;i<sz1;i++)
176 // std::size_t sz2=_the_matrix_st[i].size();
177 // _the_deno_st[i].resize(sz2);
178 // for(std::size_t j=0;j<sz2;j++)
181 // SparseDoubleVec& mToFill=_the_deno_st[i][j];
182 // const SparseDoubleVec& m=_the_matrix_st[i][j];
183 // for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++)
184 // sum+=(*it).second;
185 // for(SparseDoubleVec::const_iterator it=m.begin();it!=m.end();it++)
186 // mToFill[(*it).first]=sum;
189 // printDenoMatrix();
192 ///*! Compute integral denominators
193 // * TO BE REVISED: needs another communication since some source areas are held non locally
195 //void OverlapMapping::computeDenoIntegral()
197 // _the_deno_st.clear();
198 // std::size_t sz1=_the_matrix_st.size();
199 // _the_deno_st.resize(sz1);
200 // for(std::size_t i=0;i<sz1;i++)
202 // std::size_t sz2=_the_matrix_st[i].size();
203 // _the_deno_st[i].resize(sz2);
204 // for(std::size_t j=0;j<sz2;j++)
206 // SparseDoubleVec& mToFill=_the_deno_st[i][j];
207 // for(SparseDoubleVec::const_iterator it=mToFill.begin();it!=mToFill.end();it++)
208 // mToFill[(*it).first] = sourceAreas;
211 // printDenoMatrix();
214 /*! Compute rev integral denominators
216 void OverlapMapping::computeDenoRevIntegral(const DataArrayDouble & targetAreas)
218 _the_deno_st.clear();
219 std::size_t sz1=_the_matrix_st.size();
220 _the_deno_st.resize(sz1);
221 const double * targetAreasP = targetAreas.getConstPointer();
222 for(std::size_t i=0;i<sz1;i++)
224 std::size_t sz2=_the_matrix_st[i].size();
225 _the_deno_st[i].resize(sz2);
226 for(std::size_t j=0;j<sz2;j++)
228 SparseDoubleVec& mToFill=_the_deno_st[i][j];
229 SparseDoubleVec& mToIterate=_the_matrix_st[i][j];
230 for(SparseDoubleVec::const_iterator it=mToIterate.begin();it!=mToIterate.end();it++)
231 mToFill[(*it).first] = targetAreasP[j];
234 // printDenoMatrix();
239 * Compute denominators for ConvervativeVolumic interp.
241 void OverlapMapping::computeDenoConservativeVolumic(mcIdType nbOfTuplesTrg)
243 int myProcId=_group.myRank();
245 _the_deno_st.clear();
246 std::size_t sz1=_the_matrix_st.size();
247 _the_deno_st.resize(sz1);
248 std::vector<double> deno(nbOfTuplesTrg);
249 // Fills in the vector indexed by target cell ID:
250 for(std::size_t i=0;i<sz1;i++)
252 const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i];
253 int curSrcId=_the_matrix_st_source_proc_id[i];
254 map < int, MCAuto<DataArrayIdType> >::const_iterator isItem1 = _sent_trg_ids.find(curSrcId);
256 if(isItem1==_sent_trg_ids.end() || curSrcId==myProcId) // Local computation: simple, because rowId of mat are directly target cell ids.
258 for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
259 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
260 deno[rowId]+=(*it2).second;
262 else // matrix was received, remote computation
264 const DataArrayIdType *trgIds = (*isItem1).second;
265 const mcIdType *trgIds2=trgIds->getConstPointer();
266 for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
267 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
268 deno[trgIds2[rowId]]+=(*it2).second;
271 // Broadcast the vector into a structure similar to the initial sparse matrix of numerators:
272 for(std::size_t i=0;i<sz1;i++)
275 const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i];
276 int curSrcId=_the_matrix_st_source_proc_id[i];
277 map < int, MCAuto<DataArrayIdType> >::const_iterator isItem1 = _sent_trg_ids.find(curSrcId);
278 std::vector< SparseDoubleVec >& denoM=_the_deno_st[i];
279 denoM.resize(mat.size());
280 if(isItem1==_sent_trg_ids.end() || curSrcId==myProcId)//item1 of step2 main algo. Simple, because rowId of mat are directly target ids.
283 for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId2++)
284 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
285 denoM[rowId2][(*it2).first]=deno[rowId2];
289 const DataArrayIdType *trgIds = (*isItem1).second;
290 const mcIdType *trgIds2=trgIds->getConstPointer();
291 for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++,rowId++)
292 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
293 denoM[rowId][(*it2).first]=deno[trgIds2[rowId]];
296 // printDenoMatrix();
300 * This method performs step #0/3 in serialization process.
301 * \param count tells specifies nb of elems to send to corresponding proc id. size equal to _group.size().
302 * \param offsets tells for a proc i where to start serialize#0 matrix. size equal to _group.size().
303 * \param nbOfElemsSrc of size _group.size(). Comes from previous all2all call. tells how many srcIds per proc contains matrix for current proc.
305 void OverlapMapping::serializeMatrixStep0ST(const mcIdType *nbOfElemsSrc, mcIdType *&bigArr, int *count, int *offsets,
306 int *countForRecv, int *offsetsForRecv) const
308 std::size_t grpSize=_group.size();
309 std::fill<int *>(count,count+grpSize,0);
311 int myProcId=_group.myRank();
312 for(std::size_t i=0;i<_matrixes_st.size();i++)
314 if(_source_proc_id_st[i]==myProcId)// && _target_proc_id_st[i]!=myProcId
316 count[_target_proc_id_st[i]]=(int)_matrixes_st[i].size()+1;
317 szz+=_matrixes_st[i].size()+1;
320 bigArr=new mcIdType[szz];
322 for(std::size_t i=1;i<grpSize;i++)
323 offsets[i]=offsets[i-1]+count[i-1];
324 for(std::size_t i=0;i<_matrixes_st.size();i++)
326 if(_source_proc_id_st[i]==myProcId)
328 mcIdType start=offsets[_target_proc_id_st[i]];
329 mcIdType *work=bigArr+start;
331 const std::vector< SparseDoubleVec >& mat=_matrixes_st[i];
332 for(std::vector< SparseDoubleVec >::const_iterator it=mat.begin();it!=mat.end();it++,work++)
333 work[1]=work[0]+ToIdType((*it).size());
338 for(std::size_t i=0;i<grpSize;i++)
340 if(nbOfElemsSrc[i]>0)
341 countForRecv[i]=(int)nbOfElemsSrc[i]+1;
345 offsetsForRecv[i]=offsetsForRecv[i-1]+countForRecv[i-1];
350 * This method performs step#1 and step#2/3. It returns the size of expected array to get allToAllV.
351 * It is where the locally computed matrices are serialized to be sent to adequate final proc.
353 mcIdType OverlapMapping::serializeMatrixStep1ST(const mcIdType *nbOfElemsSrc, const mcIdType *recvStep0, const int *countStep0, const int *offsStep0,
354 mcIdType *&bigArrI, double *&bigArrD, int *count, int *offsets,
355 int *countForRecv, int *offsForRecv) const
357 std::size_t grpSize=_group.size();
358 int myProcId=_group.myRank();
361 for(std::size_t i=0;i<grpSize;i++)
363 if(nbOfElemsSrc[i]!=0)
364 countForRecv[i]=(int)recvStep0[offsStep0[i]+nbOfElemsSrc[i]];
367 szz+=countForRecv[i];
369 offsForRecv[i]=offsForRecv[i-1]+countForRecv[i-1];
372 std::fill(count,count+grpSize,0);
374 std::size_t fullLgth=0;
375 for(std::size_t i=0;i<_matrixes_st.size();i++)
377 if(_source_proc_id_st[i]==myProcId)
379 const std::vector< SparseDoubleVec >& mat=_matrixes_st[i];
380 mcIdType lgthToSend=0;
381 for(std::vector< SparseDoubleVec >::const_iterator it=mat.begin();it!=mat.end();it++)
382 lgthToSend+=ToIdType((*it).size());
383 count[_target_proc_id_st[i]]=(int)lgthToSend;
384 fullLgth+=lgthToSend;
387 for(std::size_t i=1;i<grpSize;i++)
388 offsets[i]=offsets[i-1]+count[i-1];
390 bigArrI=new mcIdType[fullLgth];
391 bigArrD=new double[fullLgth];
394 for(std::size_t i=0;i<_matrixes_st.size();i++)
396 if(_source_proc_id_st[i]==myProcId)
398 const std::vector< SparseDoubleVec >& mat=_matrixes_st[i];
399 for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++)
402 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++,j++)
404 bigArrI[fullLgth+j]=(*it2).first;
405 bigArrD[fullLgth+j]=(*it2).second;
407 fullLgth+=(*it1).size();
415 * This is the last step after all2Alls for matrix exchange.
416 * _the_matrix_st is the final matrix :
417 * - The first entry is srcId in current proc.
418 * - The second is the pseudo id of source proc (correspondence with true id is in attribute _the_matrix_st_source_proc_id and _the_matrix_st_source_ids)
419 * - the third is the srcId in the pseudo source proc
421 void OverlapMapping::unserializationST(mcIdType nbOfTrgElems,
422 const mcIdType *nbOfElemsSrcPerProc,//first all2all
423 const mcIdType *bigArrRecv, const int *bigArrRecvCounts, const int *bigArrRecvOffs,//2nd all2all
424 const mcIdType *bigArrRecv2, const double *bigArrDRecv2, const int *bigArrRecv2Count, const int *bigArrRecv2Offs)//3rd and 4th all2alls
426 _the_matrix_st.clear();
427 _the_matrix_st_source_proc_id.clear();
429 std::size_t grpSize=_group.size();
430 for(unsigned int i=0;i<grpSize;i++)
431 if(nbOfElemsSrcPerProc[i]!=0)
432 _the_matrix_st_source_proc_id.push_back(i);
433 std::size_t nbOfPseudoProcs=_the_matrix_st_source_proc_id.size();//_the_matrix_st_target_proc_id.size() contains number of matrix fetched remotely whose sourceProcId==myProcId
434 _the_matrix_st.resize(nbOfPseudoProcs);
437 for(std::size_t i=0;i<grpSize;i++)
438 if(nbOfElemsSrcPerProc[i]!=0)
440 _the_matrix_st[j].resize(nbOfElemsSrcPerProc[i]);
441 for(mcIdType k=0;k<nbOfElemsSrcPerProc[i];k++)
443 mcIdType offs=bigArrRecv[bigArrRecvOffs[i]+k];
444 mcIdType lgthOfMap=bigArrRecv[bigArrRecvOffs[i]+k+1]-offs;
445 for(mcIdType l=0;l<lgthOfMap;l++)
446 _the_matrix_st[j][k][bigArrRecv2[bigArrRecv2Offs[i]+offs+l]]=bigArrDRecv2[bigArrRecv2Offs[i]+offs+l];
453 * This method should be called when all remote matrix with sourceProcId==thisProcId have been retrieved and are
454 * in 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id' and 'this->_the_matrix_st_target_ids'.
455 * This method finish the job of filling 'this->_the_matrix_st' and 'this->_the_matrix_st_target_proc_id'
456 * by putting candidates in 'this->_matrixes_st' into them (i.e. local computation result).
458 void OverlapMapping::finishToFillFinalMatrixST()
460 int myProcId=_group.myRank();
461 std::size_t sz=_matrixes_st.size();
462 int nbOfEntryToAdd=0;
463 for(std::size_t i=0;i<sz;i++)
464 if(_source_proc_id_st[i]!=myProcId)
466 if(nbOfEntryToAdd==0)
468 std::size_t oldNbOfEntry=_the_matrix_st.size();
469 std::size_t newNbOfEntry=oldNbOfEntry+nbOfEntryToAdd;
470 _the_matrix_st.resize(newNbOfEntry);
471 std::size_t j=oldNbOfEntry;
472 for(std::size_t i=0;i<sz;i++)
473 if(_source_proc_id_st[i]!=myProcId)
475 const std::vector<SparseDoubleVec >& mat=_matrixes_st[i];
476 _the_matrix_st[j]=mat;
477 _the_matrix_st_source_proc_id.push_back(_source_proc_id_st[i]);
484 * This method performs a transpose multiply of 'fieldInput' and put the result into 'fieldOutput'.
485 * 'fieldInput' is expected to be the sourcefield and 'fieldOutput' the targetfield.
487 void OverlapMapping::multiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput, double default_val) const
491 std::size_t nbOfCompo=fieldInput->getNumberOfComponents();//to improve same number of components to test
492 CommInterface commInterface=_group.getCommInterface();
493 const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
494 const MPI_Comm *comm=group->getComm();
495 int grpSize=_group.size();
496 int myProcID=_group.myRank();
498 INTERP_KERNEL::AutoPtr<int> nbsend=new int[grpSize];
499 INTERP_KERNEL::AutoPtr<int> nbsend2=new int[grpSize];
500 INTERP_KERNEL::AutoPtr<int> nbrecv=new int[grpSize];
501 INTERP_KERNEL::AutoPtr<int> nbrecv2=new int[grpSize];
502 fill<int *>(nbsend,nbsend+grpSize,0);
503 fill<int *>(nbrecv,nbrecv+grpSize,0);
506 vector<double> valsToSend;
510 * We call the 'BB source IDs' (bounding box source IDs) the set of source cell IDs transmitted just based on the bounding box information.
511 * This is potentially bigger than what is finally in the interp matrix and this is stored in _sent_src_ids.
512 * 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.
514 for(int procID=0;procID<grpSize;procID++)
516 /* SENDING part: compute field values to be SENT (and how many of them)
517 * - for all proc 'procID' in group
518 * * if procID == myProcID, send nothing
519 * * elif 'procID' in _proc_ids_to_send_vector_st (computed from the BB intersection)
520 * % if myProcID computed the job (myProcID, procID)
521 * => send only 'interp source IDs' field values (i.e. IDs stored in _src_ids_zip_comp)
522 * % else (=we just sent mesh data to procID, but have never seen the matrix, i.e. matrix was computed remotely by procID)
523 * => send 'BB source IDs' set of field values (i.e. IDs stored in _sent_src_ids)
525 if (procID == myProcID)
528 if(find(_proc_ids_to_send_vector_st.begin(),_proc_ids_to_send_vector_st.end(),procID)!=_proc_ids_to_send_vector_st.end())
530 MCAuto<DataArrayDouble> vals;
531 if(_locator.isInMyTodoList(myProcID, procID))
533 map<int, vector<mcIdType> >::const_iterator isItem11 = _src_ids_zip_comp.find(procID);
534 if (isItem11 == _src_ids_zip_comp.end())
535 throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: SEND: unexpected end iterator in _src_ids_zip_comp!");
536 const vector<mcIdType> & v = (*isItem11).second;
537 std::size_t sz = v.size();
538 vals=fieldInput->getArray()->selectByTupleId(&(v[0]),&(v[0])+sz);
542 map < int, MCAuto<DataArrayIdType> >::const_iterator isItem11 = _sent_src_ids.find( procID );
543 if (isItem11 == _sent_src_ids.end())
544 throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: SEND: unexpected end iterator in _sent_src_ids!");
545 vals=fieldInput->getArray()->selectByTupleId(*(*isItem11).second);
547 nbsend[procID] = (int)vals->getNbOfElems(); // nb of elem = nb_tuples*nb_compo
548 // Flat version of values to send:
549 valsToSend.insert(valsToSend.end(),vals->getConstPointer(),vals->getConstPointer()+nbsend[procID]);
552 /* RECEIVE: compute number of field values to be RECEIVED
553 * - for all proc 'procID' in group
554 * * if procID == myProcID, rcv nothing
555 * * elif 'procID' in _proc_ids_to_recv_vector_st (computed from BB intersec)
556 * % if myProcID computed the job (procID, myProcID)
557 * => receive full set ('BB source IDs') of field data from proc #procID which has never seen the matrix
558 * i.e. prepare to receive the numb in _nb_of_rcv_src_ids
559 * % else (=we did NOT compute the job, hence procID has, and knows the matrix)
560 * => receive 'interp source IDs' set of field values
562 const std::vector< int > & _proc_ids_to_recv_vector_st = _the_matrix_st_source_proc_id;
563 if (procID == myProcID)
566 if(find(_proc_ids_to_recv_vector_st.begin(),_proc_ids_to_recv_vector_st.end(),procID)!=_proc_ids_to_recv_vector_st.end())
568 if(_locator.isInMyTodoList(procID, myProcID))
570 map <int,mcIdType>::const_iterator isItem11 = _nb_of_rcv_src_ids.find(procID);
571 if (isItem11 == _nb_of_rcv_src_ids.end())
572 throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: RCV: unexpected end iterator in _nb_of_rcv_src_ids!");
573 nbrecv[procID] = (int)((*isItem11).second * nbOfCompo);
577 map<int, vector<mcIdType> >::const_iterator isItem22 = _src_ids_zip_recv.find(procID);
578 if (isItem22 == _src_ids_zip_recv.end())
579 throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: RCV: unexpected end iterator in _src_ids_zip_recv!");
580 nbrecv[procID] = (int)((*isItem22).second.size() * nbOfCompo);
584 // Compute offsets in the sending/receiving array.
585 for(int i=1;i<grpSize;i++)
587 nbsend2[i]=nbsend2[i-1]+nbsend[i-1];
588 nbrecv2[i]=nbrecv2[i-1]+nbrecv[i-1];
590 INTERP_KERNEL::AutoPtr<double> bigArr=new double[nbrecv2[grpSize-1]+nbrecv[grpSize-1]];
594 scout << "(" << myProcID << ") nbsend :" << nbsend[0] << "," << nbsend[1] << "," << nbsend[2] << "\n";
595 scout << "(" << myProcID << ") nbrecv :" << nbrecv[0] << "," << nbrecv[1] << "," << nbrecv[2] << "\n";
596 scout << "(" << myProcID << ") valsToSend: ";
597 for (int iii=0; iii<valsToSend.size(); iii++)
598 scout << ", " << valsToSend[iii];
599 cout << scout.str() << "\n";
603 * *********************** ALL-TO-ALL
605 commInterface.allToAllV(&valsToSend[0],nbsend,nbsend2,MPI_DOUBLE,
606 bigArr,nbrecv,nbrecv2,MPI_DOUBLE,*comm);
608 MPI_Barrier(MPI_COMM_WORLD);
609 scout << "(" << myProcID << ") bigArray: ";
610 for (int iii=0; iii<nbrecv2[grpSize-1]+nbrecv[grpSize-1]; iii++)
611 scout << ", " << bigArr[iii];
612 cout << scout.str() << "\n";
616 * TARGET FIELD COMPUTATION (matrix-vec computation)
618 fieldOutput->getArray()->fillWithZero();
619 INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfCompo];
621 // By default field value set to default value - so mark which cells are hit
622 mcIdType ntup = fieldOutput->getNumberOfTuples();
623 INTERP_KERNEL::AutoPtr<bool> hit_cells = new bool[ntup];
624 std::fill((bool *)hit_cells, (bool *)hit_cells+ntup, false);
626 for(vector<int>::const_iterator itProc=_the_matrix_st_source_proc_id.begin(); itProc != _the_matrix_st_source_proc_id.end();itProc++)
627 // For each source processor corresponding to a locally held matrix:
629 int srcProcID = *itProc;
630 std::size_t id = std::distance(_the_matrix_st_source_proc_id.begin(),itProc);
631 const vector< SparseDoubleVec >& mat =_the_matrix_st[id];
632 const vector< SparseDoubleVec >& deno = _the_deno_st[id];
634 /* FINAL MULTIPLICATION
635 * * if srcProcID == myProcID, local multiplication without any mapping
636 * => for all target cell ID 'tgtCellID'
637 * => for all src cell ID 'srcCellID' in the sparse vector
638 * => tgtFieldLocal[tgtCellID] += srcFieldLocal[srcCellID] * matrix[tgtCellID][srcCellID] / deno[tgtCellID][srcCellID]
640 if (srcProcID == myProcID)
642 std::size_t nbOfTrgTuples=mat.size();
643 double * targetBase = fieldOutput->getArray()->getPointer();
644 for(std::size_t j=0; j<nbOfTrgTuples; j++)
646 const SparseDoubleVec& mat1=mat[j];
647 const SparseDoubleVec& deno1=deno[j];
648 SparseDoubleVec::const_iterator it5=deno1.begin();
649 const double * localSrcField = fieldInput->getArray()->getConstPointer();
650 double * targetPt = targetBase+j*nbOfCompo;
651 for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
653 // Apply the multiplication for all components:
654 double ratio = (*it3).second/(*it5).second;
655 transform(localSrcField+((*it3).first)*nbOfCompo,
656 localSrcField+((*it3).first+1)*nbOfCompo,
658 [=](double d) { return d*ratio; });
659 // Accumulate with current value:
660 transform((double *)tmp,(double *)tmp+nbOfCompo,targetPt,targetPt,plus<double>());
666 if(nbrecv[srcProcID]<=0) // also covers the preceding 'if'
669 /* * if something was received
670 * % if received matrix (=we didn't compute the job), this means that :
671 * 1. we sent part of our targetIDs to srcProcID before, so that srcProcId can do the computation.
672 * 2. srcProcID has sent us only the 'interp source IDs' field values
673 * => invert _src_ids_zip_recv -> 'revert_zip'
674 * => for all target cell ID 'tgtCellID'
675 * => mappedTgtID = _sent_trg_ids[srcProcID][tgtCellID]
676 * => for all src cell ID 'srcCellID' in the sparse vector
677 * => idx = revert_zip[srcCellID]
678 * => tgtFieldLocal[mappedTgtID] += rcvValue[srcProcID][idx] * matrix[tgtCellID][srcCellID] / deno[tgtCellID][srcCellID]
680 if(!_locator.isInMyTodoList(srcProcID, myProcID))
682 // invert _src_ids_zip_recv
683 map<mcIdType,int> revert_zip;
684 map<int, vector<mcIdType> >::const_iterator it11= _src_ids_zip_recv.find(srcProcID);
685 if (it11 == _src_ids_zip_recv.end())
686 throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: MULTIPLY: unexpected end iterator in _src_ids_zip_recv!");
688 const vector<mcIdType> & vec = (*it11).second;
690 for(vector<mcIdType>::const_iterator it=vec.begin();it!=vec.end();it++,newId++)
691 revert_zip[*it]=newId;
692 map < int, MCAuto<DataArrayIdType> >::const_iterator isItem24 = _sent_trg_ids.find(srcProcID);
693 if (isItem24 == _sent_trg_ids.end())
694 throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: MULTIPLY: unexpected end iterator in _sent_trg_ids!");
695 const DataArrayIdType *tgrIdsDA = (*isItem24).second;
696 const mcIdType *tgrIds = tgrIdsDA->getConstPointer();
698 std::size_t nbOfTrgTuples=mat.size();
699 double * targetBase = fieldOutput->getArray()->getPointer();
700 for(std::size_t j=0;j<nbOfTrgTuples;j++)
702 const SparseDoubleVec& mat1=mat[j];
703 const SparseDoubleVec& deno1=deno[j];
704 SparseDoubleVec::const_iterator it5=deno1.begin();
705 double * targetPt = targetBase+tgrIds[j]*nbOfCompo;
706 for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
708 map<mcIdType,int>::const_iterator it4=revert_zip.find((*it3).first);
709 if(it4==revert_zip.end())
710 throw INTERP_KERNEL::Exception("OverlapMapping::multiply(): internal error: MULTIPLY: unexpected end iterator in revert_zip!");
711 double ratio = (*it3).second/(*it5).second;
712 transform(bigArr+nbrecv2[srcProcID]+((*it4).second)*nbOfCompo,
713 bigArr+nbrecv2[srcProcID]+((*it4).second+1)*nbOfCompo,
715 [=](double d) { return d*ratio; } );
716 transform((double *)tmp,(double *)tmp+nbOfCompo,targetPt,targetPt,plus<double>());
717 hit_cells[tgrIds[j]] = true;
722 /* % else (=we computed the job and we received the 'BB source IDs' set of source field values)
723 * => for all target cell ID 'tgtCellID'
724 * => for all src cell ID 'srcCellID' in the sparse vector
725 * => tgtFieldLocal[tgtCellID] += rcvValue[srcProcID][srcCellID] * matrix[tgtCellID][srcCellID] / deno[tgtCellID][srcCellID]
728 // Same loop as in the case srcProcID == myProcID, except that instead of working on local field data, we work on bigArr
729 std::size_t nbOfTrgTuples=mat.size();
730 double * targetBase = fieldOutput->getArray()->getPointer();
731 for(std::size_t j=0;j<nbOfTrgTuples;j++)
733 const SparseDoubleVec& mat1=mat[j];
734 const SparseDoubleVec& deno1=deno[j];
735 SparseDoubleVec::const_iterator it5=deno1.begin();
736 double * targetPt = targetBase+j*nbOfCompo;
737 for(SparseDoubleVec::const_iterator it3=mat1.begin();it3!=mat1.end();it3++,it5++)
739 // Apply the multiplication for all components:
740 double ratio = (*it3).second/(*it5).second;
741 transform(bigArr+nbrecv2[srcProcID]+((*it3).first)*nbOfCompo,
742 bigArr+nbrecv2[srcProcID]+((*it3).first+1)*nbOfCompo,
744 [=](double d) { return d*ratio; } );
745 // Accumulate with current value:
746 transform((double *)tmp,(double *)tmp+nbOfCompo,targetPt,targetPt,plus<double>());
753 // Fill in default values for cells which haven't been hit:
755 for(bool * hit_cells_ptr=hit_cells; i< ntup; hit_cells_ptr++,i++)
756 if (!(*hit_cells_ptr))
758 double * targetPt=fieldOutput->getArray()->getPointer();
759 fill(targetPt+i*nbOfCompo, targetPt+(i+1)*nbOfCompo, default_val);
764 * This method performs a transpose multiply of 'fieldInput' and put the result into 'fieldOutput'.
765 * 'fieldInput' is expected to be the targetfield and 'fieldOutput' the sourcefield.
767 void OverlapMapping::transposeMultiply(const MEDCouplingFieldDouble *fieldInput, MEDCouplingFieldDouble *fieldOutput)
772 * This method should be called immediately after _the_matrix_st has been filled with remote computed matrix
773 * put in this proc for Matrix-Vector.
774 * It fills _src_ids_zip_recv (see member doc)
776 void OverlapMapping::fillSourceIdsZipReceivedForMultiply()
778 /* When it is called, only the bits received from other processors (i.e. the remotely executed jobs) are in the
779 big matrix _the_matrix_st. */
781 CommInterface commInterface=_group.getCommInterface();
782 int myProcId=_group.myRank();
783 std::size_t nbOfMatrixRecveived=_the_matrix_st_source_proc_id.size();
784 for(std::size_t i=0;i<nbOfMatrixRecveived;i++)
786 int curSrcProcId=_the_matrix_st_source_proc_id[i];
787 if(curSrcProcId!=myProcId) // if =, data has been populated by addContributionST()
789 const std::vector< SparseDoubleVec >& mat=_the_matrix_st[i];
790 std::set<mcIdType> s;
791 for(std::vector< SparseDoubleVec >::const_iterator it1=mat.begin();it1!=mat.end();it1++)
792 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
793 s.insert((*it2).first);
794 vector<mcIdType> vec(s.begin(),s.end());
795 _src_ids_zip_recv[curSrcProcId] = vec;
801 void OverlapMapping::printTheMatrix() const
803 CommInterface commInterface=_group.getCommInterface();
804 const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
805 const MPI_Comm *comm=group->getComm();
806 int grpSize=_group.size();
807 int myProcId=_group.myRank();
808 std::stringstream oscerr;
809 int nbOfMat=_the_matrix_st.size();
810 oscerr << "(" << myProcId << ") I hold " << nbOfMat << " matrix(ces) : "<< std::endl;
811 for(int i=0;i<nbOfMat;i++)
813 oscerr << " - Matrix #" << i << " coming from source proc #" << _the_matrix_st_source_proc_id[i] << ":\n ";
814 const std::vector< SparseDoubleVec >& locMat=_the_matrix_st[i];
816 for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++)
818 oscerr << " Target Cell #" << j;
819 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
820 oscerr << " (" << (*it2).first << "," << (*it2).second << "), ";
824 oscerr << "*********" << std::endl;
826 // Hope this will be flushed in one go:
827 std::cerr << oscerr.str() << std::endl;
829 // MPI_Barrier(MPI_COMM_WORLD);
832 void OverlapMapping::printMatrixesST() const
834 CommInterface commInterface=_group.getCommInterface();
835 const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
836 const MPI_Comm *comm=group->getComm();
837 int grpSize=_group.size();
838 int myProcId=_group.myRank();
839 std::stringstream oscerr;
840 int nbOfMat=_matrixes_st.size();
841 oscerr << "(" << myProcId << ") I hold " << nbOfMat << " LOCAL matrix(ces) : "<< std::endl;
842 for(int i=0;i<nbOfMat;i++)
844 oscerr << " - Matrix #" << i << ": (source proc #" << _source_proc_id_st[i] << " / tgt proc#" << _target_proc_id_st[i] << "): \n";
845 const std::vector< SparseDoubleVec >& locMat=_matrixes_st[i];
847 for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++)
849 oscerr << " Target Cell #" << j;
850 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
851 oscerr << " (" << (*it2).first << "," << (*it2).second << "), ";
855 oscerr << "*********" << std::endl;
857 // Hope this will be flushed in one go:
858 std::cerr << oscerr.str() << std::endl;
861 void OverlapMapping::printDenoMatrix() const
863 CommInterface commInterface=_group.getCommInterface();
864 const MPIProcessorGroup *group=static_cast<const MPIProcessorGroup*>(&_group);
865 const MPI_Comm *comm=group->getComm();
866 int grpSize=_group.size();
867 int myProcId=_group.myRank();
868 std::stringstream oscerr;
869 int nbOfMat=_the_deno_st.size();
870 oscerr << "(" << myProcId << ") I hold " << nbOfMat << " DENOMINATOR matrix(ces) : "<< std::endl;
871 for(int i=0;i<nbOfMat;i++)
873 oscerr << " - Matrix #" << i << " coming from source proc #" << _the_matrix_st_source_proc_id[i] << ": \n";
874 const std::vector< SparseDoubleVec >& locMat=_the_deno_st[i];
876 for(std::vector< SparseDoubleVec >::const_iterator it1=locMat.begin();it1!=locMat.end();it1++, j++)
878 oscerr << " Target Cell #" << j;
879 for(SparseDoubleVec::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
880 oscerr << " (" << (*it2).first << "," << (*it2).second << "), ";
884 oscerr << "*********" << std::endl;
886 // Hope this will be flushed in one go:
887 std::cerr << oscerr.str() << std::endl;