double_vector_with_halo.cc
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
31 
32 namespace oomph
33 {
34 
35  //============================================================================
36  ///\short Constructor that sets up the required information communicating
37  ///between all processors. Requires two "all to all" communications.
38  ///Arguments are the distribution of the DoubleVector and a
39  ///Vector of global unknowns required on this processor.
40  //===========================================================================
42  LinearAlgebraDistribution* const &dist_pt,
43  const Vector<unsigned> &required_global_eqn) : Distribution_pt(dist_pt)
44  {
45 #ifdef OOMPH_HAS_MPI
46  //Only bother to do anything if the vector is distributed
47  if(dist_pt->distributed())
48  {
49  //First create temporary maps for halo requests.
50  //Using the map structure ensures that the data will be sorted
51  //into processor order: the order of the unsigned integer key
52 
53  //These are the halo requests indexed by master processor and the
54  //local equation number that is to be haloed on that processor
55  std::map<unsigned,Vector<unsigned> >to_be_haloed;
56  //These are the halo requests indexed by master processor and the
57  //index in the additional storage
58  //that corresponds to the halo data on this processor
59  std::map<unsigned,Vector<unsigned> > halo_entries;
60 
61  //Find rank of current processor
62  const unsigned my_rank =
63  static_cast<int>(dist_pt->communicator_pt()->my_rank());
64  //Loop over the desired equations (on this processor)
65  //and work out whether we need to halo them according to the
66  //given distribution
67  const unsigned n_global_eqn = required_global_eqn.size();
68  //Index for the locally stored halo values
69  unsigned index = 0;
70  for(unsigned n=0;n<n_global_eqn;n++)
71  {
72  //Cache the required GLOBAL equation number
73  const unsigned i_global = required_global_eqn[n];
74  //Where is it stored?
75  unsigned rank_of_global = dist_pt->rank_of_global_row(i_global);
76  //If the equation is not stored locally then
77  //populate the two maps
78  if(my_rank!=rank_of_global)
79  {
80  //Work out the local entry on the appropriate processor
81  unsigned i_local = i_global - dist_pt->first_row(rank_of_global);
82  //Mark the local storage index as halo with rank_of_global as master
83  halo_entries[rank_of_global].push_back(index);
84  //Mark the local equation of the rank_of_global as to be
85  //haloed
86  to_be_haloed[rank_of_global].push_back(i_local);
87  //Store the local index corresponding to the global equation
88  Local_index[i_global] = index;
89  //Increment the index
90  ++index;
91  }
92  }
93 
94  //We now need to tell the other processors which of their data are
95  //haloed on this processor
96 
97  //First find out how many processors there are!
98  const int n_proc = dist_pt->communicator_pt()->nproc();
99 
100  //Setup storage for number of data to be sent to each processor
101  Vector<int> send_n(n_proc,0);
102  Vector<int> send_displacement(n_proc,0);
103  int send_data_count=0;
104 
105  //Iterate over the entries in the map
106  //This will be in rank order because of the ordering of the map
107  for(std::map<unsigned,Vector<unsigned> >::iterator it
108  = to_be_haloed.begin();it!=to_be_haloed.end();++it)
109  {
110  const unsigned rank = it->first;
111  const unsigned size_ = it->second.size();
112  //The displacement is the current number of data
113  send_displacement[rank] = send_data_count;
114  //The number to send is just the size of the array
115  send_n[rank] = static_cast<int> (size_);
116  send_data_count += size_;
117  }
118 
119  //Now send the number of haloed entries from every processor
120  //to every processor
121 
122  //Receive the data directly into the storage for haloed daa
123  Haloed_n.resize(n_proc,0);
124  MPI_Alltoall(&send_n[0],1,MPI_INT,&Haloed_n[0],1,MPI_INT,
125  dist_pt->communicator_pt()->mpi_comm());
126 
127  //Prepare the data to be sent
128  //Always resize to at least one
129  if(send_data_count==0) {++send_data_count;}
130  Vector<unsigned> send_data(send_data_count);
131  //Iterate over the entries in the map
132  unsigned count=0;
133  for(std::map<unsigned,Vector<unsigned> >::iterator it
134  = to_be_haloed.begin();it!=to_be_haloed.end();++it)
135  {
136  //Iterate over the vector
137  for(Vector<unsigned>::iterator it2 = it->second.begin();
138  it2 != it->second.end();++it2)
139  {
140  send_data[count] = (*it2);
141  ++count;
142  }
143  }
144 
145  //Prepare the data to be received,
146  //Again this can go directly into Haloed storage
147  int receive_data_count=0;
148  Haloed_displacement.resize(n_proc);
149  for(int d=0;d<n_proc;d++)
150  {
151  //The displacement is the amount of data received so far
152  Haloed_displacement[d]=receive_data_count;
153  receive_data_count += Haloed_n[d];
154  }
155 
156  //Now resize the receive buffer
157  //Always make sure that it has size of at least one
158  if(receive_data_count==0) {++receive_data_count;}
159  Haloed_eqns.resize(receive_data_count);
160  //Send the data between all the processors
161  MPI_Alltoallv(&send_data[0],&send_n[0],&send_displacement[0],
162  MPI_UNSIGNED,
163  &Haloed_eqns[0],&Haloed_n[0],
165  MPI_UNSIGNED,
166  dist_pt->communicator_pt()->mpi_comm());
167 
168  //Finally, we translate the map of halo entries into the permanent
169  //storage
170  Halo_n.resize(n_proc,0);
171  Halo_displacement.resize(n_proc,0);
172 
173  //Loop over all the entries in the map
174  unsigned receive_haloed_count=0;
175  for(int d=0;d<n_proc;d++)
176  {
177  //Pointer to the map entry
178  std::map<unsigned,Vector<unsigned> >::iterator it
179  = halo_entries.find(d);
180  //If we don't have it in the map, skip
181  if(it==halo_entries.end())
182  {
183  Halo_displacement[d] = receive_haloed_count;
184  Halo_n[d] = 0;
185  }
186  else
187  {
188  Halo_displacement[d] = receive_haloed_count;
189  const int size_ = it->second.size();
190  Halo_n[d] = size_;
191  //Resize the equations to be sent
192  Halo_eqns.resize(receive_haloed_count+size_);
193  for(int i=0;i<size_;i++)
194  {
195  Halo_eqns[receive_haloed_count+i] = it->second[i];
196  }
197  receive_haloed_count += size_;
198  }
199  }
200  }
201 #endif
202  }
203 
204  //=====================================================================
205  ///\short Function that sets up a vector of pointers to halo
206  /// data, index using the scheme in Local_index. The first arguement
207  /// is a map of pointers to all halo data index by the global equation
208  /// number
209  //====================================================================
211  const std::map<unsigned,double*> &halo_data_pt, Vector<double*> &halo_dof_pt)
212  {
213  //How many entries are there in the map
214  unsigned n_halo = Local_index.size();
215  //Resize the vector
216  halo_dof_pt.resize(n_halo);
217 
218  //Loop over all the entries in the map
219  for(std::map<unsigned,unsigned>::iterator it=Local_index.begin();
220  it!=Local_index.end();++it)
221  {
222  //Find the pointer in the halo_data_pt map
223  std::map<unsigned,double*>::const_iterator it2 =
224  halo_data_pt.find(it->first);
225  //Did we find it
226  if(it2!=halo_data_pt.end())
227  {
228  //Now set the entry
229  halo_dof_pt[it->second] = it2->second;
230  }
231  else
232  {
233  std::ostringstream error_stream;
234  error_stream << "Global equation " << it->first
235  << " reqired as halo is not stored in halo_data_pt\n";
236 
237  throw OomphLibError(error_stream.str(),
238  OOMPH_CURRENT_FUNCTION,
239  OOMPH_EXCEPTION_LOCATION);
240  }
241  }
242  }
243 
244  //------------------------------------------------------------------
245  //Member functions for the DoubleVectorWithHaloEntries
246  //-------------------------------------------------------------------
247 
248 
249  //=========================================================================
250  ///Synchronise the halo data within the vector. This requires one
251  ///"all to all" communnication.
252  //====================================================================
254  {
255 #ifdef OOMPH_HAS_MPI
256  //Only need to do anything if the DoubleVector is distributed
257  if(this->distributed())
258  {
259  //Read out the number of entries to send
260  const unsigned n_send = Halo_scheme_pt->Haloed_eqns.size();
261  Vector<double> send_data(n_send);
262  //Read out the data values
263  for(unsigned i=0;i<n_send;i++)
264  {
265  send_data[i] = (*this)[Halo_scheme_pt->Haloed_eqns[i]];
266  }
267 
268  //Read out the number of entries to receive
269  const unsigned n_receive = Halo_scheme_pt->Halo_eqns.size();
270  Vector<double> receive_data(n_receive);
271 
272  //Make sure that the send and receive data have size at least one
273  if(n_send==0) {send_data.resize(1);}
274  if(n_receive==0) {receive_data.resize(1);}
275  //Communicate
276  MPI_Alltoallv(&send_data[0],&Halo_scheme_pt->Haloed_n[0],
277  &Halo_scheme_pt->Haloed_displacement[0],MPI_DOUBLE,
278  &receive_data[0],&Halo_scheme_pt->Halo_n[0],
279  &Halo_scheme_pt->Halo_displacement[0],MPI_DOUBLE,
280  this->distribution_pt()->communicator_pt()->mpi_comm());
281 
282 
283  //Now I need simply to update my local values
284  for(unsigned i=0;i<n_receive;i++)
285  {
286  Halo_value[Halo_scheme_pt->Halo_eqns[i]] =receive_data[i];
287  }
288  }
289 #endif
290  }
291 
292  //=========================================================================
293  ///Gather all ther data from multiple processors and sum the result
294  /// which will be stored in the master copy and then synchronised to
295  /// all copies. This requires two "all to all" communications
296  //====================================================================
298  {
299 #ifdef OOMPH_HAS_MPI
300  //Only need to do anything if the DoubleVector is distributed
301  if(this->distributed())
302  {
303  //Send the Halo entries to the master processor
304  const unsigned n_send = Halo_scheme_pt->Halo_eqns.size();
305  Vector<double> send_data(n_send);
306  //Read out the data values
307  for(unsigned i=0;i<n_send;i++)
308  {
309  send_data[i] = Halo_value[Halo_scheme_pt->Halo_eqns[i]];
310  }
311 
312  //Read out the number of entries to receive
313  const unsigned n_receive = Halo_scheme_pt->Haloed_eqns.size();
314  Vector<double> receive_data(n_receive);
315 
316  //Make sure that the send and receive data have size at least one
317  if(n_send==0) {send_data.resize(1);}
318  if(n_receive==0) {receive_data.resize(1);}
319  //Communicate
320  MPI_Alltoallv(&send_data[0],&Halo_scheme_pt->Halo_n[0],
321  &Halo_scheme_pt->Halo_displacement[0],MPI_DOUBLE,
322  &receive_data[0],&Halo_scheme_pt->Haloed_n[0],
323  &Halo_scheme_pt->Haloed_displacement[0],MPI_DOUBLE,
324  this->distribution_pt()->communicator_pt()->mpi_comm());
325 
326 
327  //Now I need simply to update and sum my local values
328  for(unsigned i=0;i<n_receive;i++)
329  {
330  (*this)[Halo_scheme_pt->Haloed_eqns[i]] += receive_data[i];
331  }
332 
333  //Then synchronise
334  this->synchronise();
335  }
336 #endif
337  }
338 
339 
340  //===================================================================
341  ///Construct the halo scheme and storage for the halo data
342 //=====================================================================
344  DoubleVectorHaloScheme* const &halo_scheme_pt)
345  {
347 
348  if(Halo_scheme_pt!=0)
349  {
350  //Need to set up the halo data
351  unsigned n_halo_data = halo_scheme_pt->Local_index.size();
352 
353  //Resize the halo storage
354  Halo_value.resize(n_halo_data);
355 
356  //Now let's get the initial values from the other processors
357  this->synchronise();
358  }
359  }
360 
361 
362 }
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero...
cstr elem_len * i
Definition: cfortran.h:607
Vector< int > Halo_n
Storage for the number of entries to be received from each other processor.
DoubleVectorHaloScheme * Halo_scheme_pt
Pointer to the lookup scheme that stores information about on which processor the required informatio...
bool distributed() const
distribution is serial or distributed
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed ...
Vector< int > Halo_displacement
Storage for the offsets of the processor data in the receive buffer.
unsigned rank_of_global_row(const unsigned i) const
return the processor rank of the global row number i
Vector< unsigned > Halo_eqns
Storage for all the entries that are to be received from other processors (received_from_proc0,received_from_proc1,...received_from_procn)
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
DoubleVectorHaloScheme *& halo_scheme_pt()
Access function for halo scheme.
Vector< double > Halo_value
Vector of the halo values.
void build_halo_scheme(DoubleVectorHaloScheme *const &halo_scheme_pt)
Construct the halo scheme and storage for the halo data.
Vector< int > Haloed_displacement
Storage for the offsets of the haloed entries for each processor in the packed Haloed_eqns array...
Vector< unsigned > Haloed_eqns
The haloed entries that will be sent in a format compatible with MPI_Alltoallv i.e. (send_to_proc0,send_to_proc1 ... send_to_procn)
void synchronise()
Synchronise the halo data.
void setup_halo_dofs(const std::map< unsigned, double * > &halo_data_pt, Vector< double * > &halo_dof_pt)
Function that sets up a vector of pointers to halo data, index using the scheme in Local_index...
Vector< int > Haloed_n
Storage for the number of haloed entries to be sent to each processor.
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
std::map< unsigned, unsigned > Local_index
Storage for the translation scheme from global unknown to local index in the additional storage vecto...
DoubleVectorHaloScheme(LinearAlgebraDistribution *const &dist_pt, const Vector< unsigned > &required_global_eqn)
Constructor that sets up the required information communicating between all processors. Requires two "all to all" communications. Arguments are the distribution of the DoubleVector and a Vector of global unknowns required on this processor.