src/ParaMEDMEM/BlockTopology.cxx

   1 //  Copyright (C) 2007-2008  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 #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"
  27
  28 #include <vector>
  29 #include <algorithm>
  30 #include <utility>
  31 #include <iostream>
  32
  33 using namespace std;
  34
  35 namespace ParaMEDMEM
  36 {
  37
  38   //!converts a pair <subdomainid,local> to a global number
  39   std::pair<int,int> BlockTopology::globalToLocal(const int global) const
  40   {
  41     int subdomain_id=0;
  42     int position=global;
  43     int size=_nb_elems;
  44     int size_procs=_proc_group->size();
  45     int increment=size;
  46     vector<int>axis_position(_dimension);
  47     vector<int>axis_offset(_dimension);
  48     for (int idim=0; idim<_dimension; idim++)
  49       {
  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;
  56         int iaxis=1;
  57         while (_local_array_indices[idim][iaxis]<=axis_pos)
  58           {
  59             subdomain_id+=proc_increment;
  60             iaxis++;
  61           }
  62         axis_position[idim]=axis_pos-_local_array_indices[idim][iaxis-1];
  63         axis_offset[idim]=iaxis;
  64       }
  65     int local=0;
  66     int local_increment=1;
  67     for (int idim=_dimension-1; idim>=0; idim--)
  68       {
  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];
  71       }
  72     return make_pair(subdomain_id,local);
  73   }
  74
  75   //!converts local number to a global number
  76   int BlockTopology::localToGlobal(const pair<int,int> local) const
  77   {
  78
  79     int subdomain_id=local.first;
  80     int global=0;
  81     int loc=local.second;
  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++)
  86       {
  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;
  98       }
  99     return global;
 100   }
 101
 102   //Retrieves the local number of elements
 103   int BlockTopology::getNbLocalElements()const
 104   {
 105     int position=_proc_group->myRank();
 106     int nb_elem = 1;
 107     int increment=1;
 108     for (int i=_dimension-1; i>=0; i--)
 109       {
 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);
 116       }
 117     return nb_elem;
 118   }
 119
 120   /*!
 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.
 125    */
 126   BlockTopology::BlockTopology(const ProcessorGroup& group, MEDCouplingSMesh *grid):
 127     _proc_group(&group), _dimension(grid->getSpaceDimension()), _owns_processor_group(false)
 128   {
 129     vector <int> axis_length(_dimension);
 130     _nb_elems=1;
 131     for (int idim=0; idim <_dimension; idim++)
 132       {
 133         DataArrayDouble *arr=grid->getCoordsAt(idim);
 134         axis_length[idim]=arr->getNbOfElems();
 135         _nb_elems*=axis_length[idim];
 136       }
 137     //default splitting along 1st dimension
 138     _local_array_indices.resize(_dimension);
 139     _nb_procs_per_dim.resize(_dimension);
 140
 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();
 144
 145     for (int i=1; i<=_proc_group->size(); i++)
 146       {
 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;
 151       }
 152     for (int i=1; i<_dimension; i++)
 153       {
 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;
 158       }
 159     _cycle_type.resize(_dimension);
 160     for (int i=0; i<_dimension; i++)
 161       _cycle_type[i]=ParaMEDMEM::Block;
 162   }
 163
 164   /*!
 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.
 171    *
 172    */
 173   BlockTopology::BlockTopology(const BlockTopology& geom_topo, const ComponentTopology& comp_topo):_owns_processor_group(false)
 174   {
 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"));
 179
 180     if (comp_topo.nbComponents()==1)
 181       {
 182         *this=geom_topo;
 183         return;
 184       }
 185     else
 186       {
 187         _dimension = geom_topo.getDimension()+1;
 188         if (comp_topo.nbBlocks()>1)
 189           _proc_group=comp_topo.getProcGroup();
 190         else
 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();
 200       }
 201   }
 202
 203   /*! Constructor for creating a one-dimensional
 204    * topology from a processor group and a local
 205    * number of elements on each processor
 206    *
 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.
 211    */
 212   BlockTopology::BlockTopology(const ProcessorGroup& group, int nb_elem):_proc_group(&group),_dimension(1),_owns_processor_group(false)
 213   {
 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();
 217     int nbtemp=nb_elem;
 218     mpi_group->getCommInterface().allGather(&nbtemp, 1, MPI_INTEGER,
 219                                             nbelems_per_proc, 1, MPI_INTEGER,
 220                                             *comm);
 221     _nb_elems=0;
 222
 223     //splitting along only dimension
 224     _local_array_indices.resize(1);
 225     _nb_procs_per_dim.resize(1);
 226
 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();
 230
 231     for (int i=1; i<=_proc_group->size(); i++)
 232       {
 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];
 236       }
 237     _cycle_type.resize(1);
 238     _cycle_type[0]=ParaMEDMEM::Block;
 239     delete[] nbelems_per_proc;
 240   }
 241
 242   BlockTopology::~BlockTopology()
 243   {
 244     if (_owns_processor_group)
 245       delete _proc_group;
 246   }
 247
 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.
 252    */
 253   std::vector<std::pair<int,int> > BlockTopology::getLocalArrayMinMax() const
 254   {
 255     vector<pair<int,int> > local_indices (_dimension);
 256     int myrank=_proc_group->myRank();
 257     int increment=1;
 258     for (int i=_dimension-1; i>=0; i--)
 259       {
 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;
 265       }
 266     return local_indices;
 267   }
 268
 269   /*! Serializes the data contained in the Block Topology
 270    * for communication purposes*/
 271   void BlockTopology::serialize(int* & serializer, int& size) const
 272   {
 273     vector<int> buffer;
 274
 275     buffer.push_back(_dimension);
 276     buffer.push_back(_nb_elems);
 277     for (int i=0; i<_dimension; i++)
 278       {
 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]);
 284       }
 285
 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++)
 291       {
 292         int world_rank=world_group.translateRank(_proc_group, i);
 293         buffer.push_back(world_rank);
 294       }
 295
 296     serializer=new int[buffer.size()];
 297     size=buffer.size();
 298     copy(buffer.begin(), buffer.end(), serializer);
 299   }
 300
 301   /*!
 302    *
 303    * Unserializes the data contained in the Block Topology
 304    * after communication. Uses the same structure as the one used for serialize()
 305    *
 306    */
 307   void BlockTopology::unserialize(const int* serializer,const CommInterface& comm_interface)
 308   {
 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++)
 319       {
 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++);
 325       }
 326     set<int> procs;
 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
 334   }
 335 }