1 // Copyright (C) 2007-2008 CEA/DEN, EDF R&D
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.
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 #include "BlockTopology.hxx"
20 #include "MemArray.hxx"
21 #include "MEDCouplingSMesh.hxx"
22 #include "CommInterface.hxx"
23 #include "ProcessorGroup.hxx"
24 #include "MPIProcessorGroup.hxx"
25 #include "ComponentTopology.hxx"
26 #include "InterpKernelUtilities.hxx"
38 //!converts a pair <subdomainid,local> to a global number
39 std::pair<int,int> BlockTopology::globalToLocal(const int global) const
44 int size_procs=_proc_group->size();
46 vector<int>axis_position(_dimension);
47 vector<int>axis_offset(_dimension);
48 for (int idim=0; idim<_dimension; idim++)
50 int axis_size=_local_array_indices[idim].size()-1;
51 int axis_nb_elem=_local_array_indices[idim][axis_size];
52 increment=increment/axis_nb_elem;
53 int proc_increment = size_procs/(axis_size);
54 int axis_pos=position/increment;
55 position=position%increment;
57 while (_local_array_indices[idim][iaxis]<=axis_pos)
59 subdomain_id+=proc_increment;
62 axis_position[idim]=axis_pos-_local_array_indices[idim][iaxis-1];
63 axis_offset[idim]=iaxis;
66 int local_increment=1;
67 for (int idim=_dimension-1; idim>=0; idim--)
69 local+=axis_position[idim]*local_increment;
70 local_increment*=_local_array_indices[idim][axis_offset[idim]]-_local_array_indices[idim][axis_offset[idim]-1];
72 return make_pair(subdomain_id,local);
75 //!converts local number to a global number
76 int BlockTopology::localToGlobal(const pair<int,int> local) const
79 int subdomain_id=local.first;
82 int increment=_nb_elems;
83 int proc_increment=_proc_group->size();
84 int local_increment=getNbLocalElements();
85 for (int idim=0; idim < _dimension; idim++)
87 int axis_size=_local_array_indices[idim].size()-1;
88 int axis_nb_elem=_local_array_indices[idim][axis_size];
89 increment=increment/axis_nb_elem;
90 proc_increment = proc_increment/(axis_size);
91 int proc_axis=subdomain_id/proc_increment;
92 subdomain_id=subdomain_id%proc_increment;
93 int local_axis_nb_elem=_local_array_indices[idim][proc_axis+1]-_local_array_indices[idim][proc_axis];
94 local_increment = local_increment/local_axis_nb_elem;
95 int iaxis=loc/local_increment+_local_array_indices[idim][proc_axis];
96 global+=increment*iaxis;
97 loc = loc%local_increment;
102 //Retrieves the local number of elements
103 int BlockTopology::getNbLocalElements()const
105 int position=_proc_group->myRank();
108 for (int i=_dimension-1; i>=0; i--)
110 increment *=_nb_procs_per_dim[i];
111 int idim=position%increment;
112 position=position/increment;
113 int imin=_local_array_indices[i][idim];
114 int imax=_local_array_indices[i][idim+1];
115 nb_elem*=(imax-imin);
121 * Constructor of a block topology from a grid.
122 * This preliminary version simply splits along the first axis
123 * instead of making the best choice with respect to the
124 * values of the different axes.
126 BlockTopology::BlockTopology(const ProcessorGroup& group, MEDCouplingSMesh *grid):
127 _proc_group(&group), _dimension(grid->getSpaceDimension()), _owns_processor_group(false)
129 vector <int> axis_length(_dimension);
131 for (int idim=0; idim <_dimension; idim++)
133 DataArrayDouble *arr=grid->getCoordsAt(idim);
134 axis_length[idim]=arr->getNbOfElems();
135 _nb_elems*=axis_length[idim];
137 //default splitting along 1st dimension
138 _local_array_indices.resize(_dimension);
139 _nb_procs_per_dim.resize(_dimension);
141 _local_array_indices[0].resize(_proc_group->size()+1);
142 _local_array_indices[0][0]=0;
143 _nb_procs_per_dim[0]=_proc_group->size();
145 for (int i=1; i<=_proc_group->size(); i++)
147 _local_array_indices[0][i]=_local_array_indices[0][i-1]+
148 axis_length[0]/_proc_group->size();
149 if (i<= axis_length[0]%_proc_group->size())
150 _local_array_indices[0][i]+=1;
152 for (int i=1; i<_dimension; i++)
154 _local_array_indices[i].resize(2);
155 _local_array_indices[i][0]=0;
156 _local_array_indices[i][1]=axis_length[i];
157 _nb_procs_per_dim[i]=1;
159 _cycle_type.resize(_dimension);
160 for (int i=0; i<_dimension; i++)
161 _cycle_type[i]=ParaMEDMEM::Block;
165 * Creation of a block topology by composing
166 * a geometrical topology and a component topology.
167 * This constructor is intended for creating fields
168 * for which the parallel distribution is made on the
169 * components of the field rather than on the geometrical
170 * partitioning of the underlying mesh.
173 BlockTopology::BlockTopology(const BlockTopology& geom_topo, const ComponentTopology& comp_topo):_owns_processor_group(false)
175 // so far, the block topology can only be created if the proc group
176 // is either on geom_topo or on comp_topo
177 if (geom_topo.getProcGroup()->size()>1 && comp_topo.nbBlocks()>1)
178 throw INTERP_KERNEL::Exception(LOCALIZED("BlockTopology cannot yet be constructed with both complex geo and components topology"));
180 if (comp_topo.nbComponents()==1)
187 _dimension = geom_topo.getDimension()+1;
188 if (comp_topo.nbBlocks()>1)
189 _proc_group=comp_topo.getProcGroup();
191 _proc_group=geom_topo.getProcGroup();
192 _local_array_indices=geom_topo._local_array_indices;
193 vector<int> comp_indices = *(comp_topo.getBlockIndices());
194 _local_array_indices.push_back(comp_indices);
195 _nb_procs_per_dim=geom_topo._nb_procs_per_dim;
196 _nb_procs_per_dim.push_back(comp_topo.nbBlocks());
197 _cycle_type=geom_topo._cycle_type;
198 _cycle_type.push_back(Block);
199 _nb_elems=geom_topo.getNbElements()*comp_topo.nbComponents();
203 /*! Constructor for creating a one-dimensional
204 * topology from a processor group and a local
205 * number of elements on each processor
207 * The function must be called only by the processors belonging
208 * to group \a group. Calling it from a processor not belonging
209 * to \a group will cause an MPI error, while calling from a subset
210 * of \a group will result in a deadlock.
212 BlockTopology::BlockTopology(const ProcessorGroup& group, int nb_elem):_proc_group(&group),_dimension(1),_owns_processor_group(false)
214 int* nbelems_per_proc = new int[group.size()];
215 const MPIProcessorGroup* mpi_group=dynamic_cast<const MPIProcessorGroup*>(_proc_group);
216 const MPI_Comm* comm=mpi_group->getComm();
218 mpi_group->getCommInterface().allGather(&nbtemp, 1, MPI_INTEGER,
219 nbelems_per_proc, 1, MPI_INTEGER,
223 //splitting along only dimension
224 _local_array_indices.resize(1);
225 _nb_procs_per_dim.resize(1);
227 _local_array_indices[0].resize(_proc_group->size()+1);
228 _local_array_indices[0][0]=0;
229 _nb_procs_per_dim[0]=_proc_group->size();
231 for (int i=1; i<=_proc_group->size(); i++)
233 _local_array_indices[0][i]=_local_array_indices[0][i-1]+
234 nbelems_per_proc[i-1];
235 _nb_elems+=nbelems_per_proc[i-1];
237 _cycle_type.resize(1);
238 _cycle_type[0]=ParaMEDMEM::Block;
239 delete[] nbelems_per_proc;
242 BlockTopology::~BlockTopology()
244 if (_owns_processor_group)
248 /*! Retrieves the min and max indices of the domain stored locally
249 * for each dimension. The output vector has the topology dimension
250 * as a size and each pair <int,int> contains min and max. Indices
251 * range from min to max-1.
253 std::vector<std::pair<int,int> > BlockTopology::getLocalArrayMinMax() const
255 vector<pair<int,int> > local_indices (_dimension);
256 int myrank=_proc_group->myRank();
258 for (int i=_dimension-1; i>=0; i--)
260 increment *=_nb_procs_per_dim[i];
261 int idim=myrank%increment;
262 local_indices[i].first=_local_array_indices[i][idim];
263 local_indices[i].second=_local_array_indices[i][idim+1];
264 cout << local_indices[i].first << " "<< local_indices[i].second<<endl;
266 return local_indices;
269 /*! Serializes the data contained in the Block Topology
270 * for communication purposes*/
271 void BlockTopology::serialize(int* & serializer, int& size) const
275 buffer.push_back(_dimension);
276 buffer.push_back(_nb_elems);
277 for (int i=0; i<_dimension; i++)
279 buffer.push_back(_nb_procs_per_dim[i]);
280 buffer.push_back(_cycle_type[i]);
281 buffer.push_back(_local_array_indices[i].size());
282 for (int j=0; j<_local_array_indices[i].size(); j++)
283 buffer.push_back(_local_array_indices[i][j]);
286 //serializing the comm group
287 int size_comm=_proc_group->size();
288 buffer.push_back(size_comm);
289 MPIProcessorGroup world_group(_proc_group->getCommInterface());
290 for (int i=0; i<size_comm;i++)
292 int world_rank=world_group.translateRank(_proc_group, i);
293 buffer.push_back(world_rank);
296 serializer=new int[buffer.size()];
298 copy(buffer.begin(), buffer.end(), serializer);
303 * Unserializes the data contained in the Block Topology
304 * after communication. Uses the same structure as the one used for serialize()
307 void BlockTopology::unserialize(const int* serializer,const CommInterface& comm_interface)
309 const int* ptr_serializer=serializer;
310 cout << "unserialize..."<<endl;
311 _dimension=*(ptr_serializer++);
312 cout << "dimension "<<_dimension<<endl;
313 _nb_elems=*(ptr_serializer++);
314 cout << "nbelems "<<_nb_elems<<endl;
315 _nb_procs_per_dim.resize(_dimension);
316 _cycle_type.resize(_dimension);
317 _local_array_indices.resize(_dimension);
318 for (int i=0; i<_dimension; i++)
320 _nb_procs_per_dim[i]=*(ptr_serializer++);
321 _cycle_type[i]=(CYCLE_TYPE)*(ptr_serializer++);
322 _local_array_indices[i].resize(*(ptr_serializer++));
323 for (int j=0; j<_local_array_indices[i].size(); j++)
324 _local_array_indices[i][j]=*(ptr_serializer++);
327 int size_comm=*(ptr_serializer++);
328 for (int i=0; i<size_comm; i++)
329 procs.insert(*(ptr_serializer++));
330 cout << "unserialize..."<<procs.size()<<endl;
331 _proc_group=new MPIProcessorGroup(comm_interface,procs);
332 _owns_processor_group=true;
333 //TODO manage memory ownership of _proc_group