linear_algebra_distribution.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// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
27
28 namespace oomph
29 {
30  //============================================================================
31  /// Sets the distribution. Takes first_row, local_nrow and
32  /// global_nrow as arguments. If global_nrow is not provided or equal to
33  /// 0 then it is computed automatically
34  //============================================================================
36  const unsigned& first_row,
37  const unsigned& local_nrow,
38  const unsigned& global_nrow)
39  {
40  // copy the communicator
41  delete Comm_pt;
42  Comm_pt = new OomphCommunicator(*comm_pt);
43
44  // get the rank and the number of processors
45  int my_rank = Comm_pt->my_rank();
46  int nproc = Comm_pt->nproc();
47
48  // resize the storage
49  First_row.clear();
50  First_row.resize(nproc);
51  Nrow_local.clear();
52  Nrow_local.resize(nproc);
53
54  // set first row and local nrow for this processor
55  First_row[my_rank] = first_row;
56  Nrow_local[my_rank] = local_nrow;
57
58 #ifdef OOMPH_HAS_MPI
59  // gather the First_row vector
60  unsigned my_nrow_local = Nrow_local[my_rank];
61  MPI_Allgather(&my_nrow_local,
62  1,
63  MPI_UNSIGNED,
64  &Nrow_local[0],
65  1,
66  MPI_UNSIGNED,
67  Comm_pt->mpi_comm());
68
69  // gather the Nrow_local vector
70  unsigned my_first_row = First_row[my_rank];
71  MPI_Allgather(&my_first_row,
72  1,
73  MPI_UNSIGNED,
74  &First_row[0],
75  1,
76  MPI_UNSIGNED,
77  Comm_pt->mpi_comm());
78 #endif
79
80  // if global nrow is not provided then compute by summing local_nrow over
81  // all processors
82  if (global_nrow == 0)
83  {
84  if (nproc == 1)
85  {
86  Nrow = local_nrow;
87  }
88  else
89  {
90  Nrow = 0;
91  for (int p = 0; p < nproc; p++)
92  {
93  Nrow += Nrow_local[p];
94  }
95  }
96  }
97  else
98  {
99  Nrow = global_nrow;
100  }
101
102  // the distribution is true, unless we only have 1 processor
103  if (nproc != 1)
104  {
105  Distributed = true;
106  }
107  else
108  {
109  Distributed = false;
110  }
111
112 #ifdef OOMPH_HAS_MPI
113 #ifdef PARANOID
114  // paranoid check that the distribution works
115
116
117  // check that none of the processors partition overlap
118  for (int p = 0; p < nproc; p++)
119  {
120  if (Nrow_local[p] > 0)
121  {
122  for (int pp = p + 1; pp < nproc; pp++)
123  {
124  if (Nrow_local[pp] > 0)
125  {
126  if ((First_row[p] >= First_row[pp] &&
127  First_row[p] < First_row[pp] + Nrow_local[pp]) ||
128  (First_row[p] + Nrow_local[p] - 1 >= First_row[pp] &&
129  First_row[p] + Nrow_local[p] - 1 <
130  First_row[pp] + Nrow_local[pp]))
131  {
132  std::ostringstream error_message;
133  error_message
134  << "The distributed rows on processor " << p
135  << " and processor " << pp << " overlap.\n"
136  << "Processor " << p << " : first_row = " << First_row[p]
137  << ", nrow = " << Nrow_local[p] << ".\n"
138  << "Processor " << pp << " : first_row = " << First_row[pp]
139  << ", nrow = " << Nrow_local[pp] << ".\n";
140  throw OomphLibWarning(
141  error_message.str(),
142  "LinearAlgebraDistribution::distribute(...)",
143  OOMPH_EXCEPTION_LOCATION);
144  }
145  }
146  }
147  }
148  }
149
150  // check that no processor has a row with a global row index greater than
151  // the number of global rows
152  for (int p = 0; p < nproc; p++)
153  {
154  if (First_row[p] + Nrow_local[p] > Nrow)
155  {
156  std::ostringstream error_message;
157  error_message << "Processor " << p << " contains rows " << First_row[p]
158  << " to " << First_row[p] + Nrow_local[p] - 1
159  << " but there are only " << Nrow << " to be distributed."
160  << std::endl;
161  throw OomphLibWarning(error_message.str(),
162  "LinearAlgebraDistribution::distribute(...)",
163  OOMPH_EXCEPTION_LOCATION);
164  }
165  }
166 #endif
167 #endif
168  }
169
170  //============================================================================
171  /// Uniformly distribute global_nrow over all processors where
172  /// processors 0 holds approximately the first
173  /// global_nrow/n_proc, processor 1 holds the next
174  /// global_nrow/n_proc and so on...
175  //============================================================================
177  const unsigned& global_nrow,
178  const bool& distribute)
179  {
180  // copy the communicator
181  delete Comm_pt;
182  Comm_pt = new OomphCommunicator(*comm_pt);
183
184  // delete existing storage
185  First_row.clear();
186  Nrow_local.clear();
187
188  // set global nrow
189  Nrow = global_nrow;
190
191  // store the distributed flag
192  Distributed = distribute;
193
194 #ifdef OOMPH_HAS_MPI
195
196  // if distributed object then compute uniform distribution
197  if (distribute == true)
198  {
199  // get the number of processors
200  int nproc = Comm_pt->nproc();
201
202  // resize the storage
203  First_row.resize(nproc);
204  Nrow_local.resize(nproc);
205
206  // compute first row
207  for (int p = 0; p < nproc; p++)
208  {
209  First_row[p] =
210  unsigned((double(p) / double(nproc)) * double(global_nrow));
211  }
212
213  // compute local nrow
214  for (int p = 0; p < nproc - 1; p++)
215  {
216  Nrow_local[p] = First_row[p + 1] - First_row[p];
217  }
218  Nrow_local[nproc - 1] = global_nrow - First_row[nproc - 1];
219  }
220 #endif
221  }
222
223  //============================================================================
224  /// helper method for the =assignment operator and copy constructor
225  //============================================================================
227  const LinearAlgebraDistribution& new_dist)
228  {
229  // delete the existing storage
230  First_row.clear();
231  Nrow_local.clear();
232
233  // if new_dist is not setup
234  if (new_dist.communicator_pt() == 0)
235  {
236  delete Comm_pt;
237  Comm_pt = 0;
238  Distributed = true;
239  Nrow = 0;
240  }
241  else
242  {
243  // copy the communicator
244  delete Comm_pt;
245  Comm_pt = new OomphCommunicator(*new_dist.communicator_pt());
246
247  // the new distribution is distributed
248  if (new_dist.distributed())
249  {
250  First_row = new_dist.first_row_vector();
251  Nrow_local = new_dist.nrow_local_vector();
252
253  Distributed = true;
254  }
255  // else if the new ditribution is not distributed
256  else
257  {
258  Distributed = false;
259  }
260  Nrow = new_dist.nrow();
261  }
262  }
263
264
265  //============================================================================
266  /// operator==
267  //============================================================================
269  const LinearAlgebraDistribution& other_dist) const
270  {
271 #ifdef OOMPH_HAS_MPI
272  // compare the communicators
273  if (!((*Comm_pt) == (*other_dist.communicator_pt())))
274  {
275  return false;
276  }
277
278  // compare Distributed
279  if (Distributed != other_dist.distributed())
280  {
281  return false;
282  }
283
284  // if not distributed compare nrow
285  if (!Distributed)
286  {
287  if (other_dist.nrow() == Nrow)
288  {
289  return true;
290  }
291  return false;
292  }
293
294  // compare
295  bool flag = true;
296  int nproc = Comm_pt->nproc();
297  for (int i = 0; i < nproc && flag == true; i++)
298  {
299  if (other_dist.first_row(i) != First_row[i] ||
300  other_dist.nrow_local(i) != Nrow_local[i])
301  {
302  flag = false;
303  }
304  }
305  return flag;
306 #else
307  if (other_dist.nrow() == Nrow)
308  {
309  return true;
310  }
311  return false;
312 #endif
313  }
314
315  //=============================================================================
316  /// output operator
317  //=============================================================================
318  std::ostream& operator<<(std::ostream& stream,
320  {
321  stream << "nrow()=" << dist.nrow() << ", first_row()=" << dist.first_row()
322  << ", nrow_local()=" << dist.nrow_local()
323  << ", distributed()=" << dist.distributed() << std::endl;
324  return stream;
325  }
326
327  //=============================================================================
328  /// Namespace for helper functions for LinearAlgebraDistributions
329  //=============================================================================
330  namespace LinearAlgebraDistributionHelpers
331  {
332  //===========================================================================
333  /// Takes a vector of LinearAlgebraDistribution objects and
334  /// concatenates them such that the nrow_local of the out_distribution
335  /// is the sum of the nrow_local of all the in_distributions and the number
336  /// of global rows of the out_distribution is the sum of the number of
337  /// global rows of all the in_distributions. This results in a permutation
338  /// of the rows in the out_distribution. Think of this in terms of
339  /// DoubleVectors, if we have DoubleVectors with distributions A and B,
340  /// distributed across two processors (p0 and p1), A: [a0] (on p0) B:
341  /// [b0] (on p0)
342  /// [a1] (on p1) [b1] (on P1),
343  ///
344  /// then the out_distribution is
345  /// [a0 (on p0)
346  /// b0] (on p0)
347  /// [a1 (on p1)
348  /// b1] (on p1),
349  ///
350  /// as opposed to
351  /// [a0 (on p0)
352  /// a1] (on p0)
353  /// [b0 (on p1)
354  /// b1] (on p1).
355  ///
356  /// Note (1): The out_distribution may not be uniformly distributed even
357  /// if the in_distributions are uniform distributions.
358  /// Try this out with two distributions of global rows 3 and 5, uniformly
359  /// distributed across two processors. Compare this against a distribution
360  /// of global row 8 distributed across two processors.
361  ///
362  /// Note (2): There is no equivalent function which takes a Vector of
363  /// LinearAlgebraDistribution objects (as opposed to pointers), there should
364  /// not be one since we do not want to invoke the assignment operator when
365  /// creating the Vector of LinearAlgebraDistribution objects.
366  //===========================================================================
368  const Vector<LinearAlgebraDistribution*>& in_distribution_pt,
369  LinearAlgebraDistribution& out_distribution)
370  {
371  // How many distributions are in in_distribution?
372  unsigned ndistributions = in_distribution_pt.size();
373
374 #ifdef PARANOID
375
376  // Are any of the in_distribution pointers null?
377  // We do not want null pointers.
378  for (unsigned dist_i = 0; dist_i < ndistributions; dist_i++)
379  {
380  if (in_distribution_pt[dist_i] == 0)
381  {
382  std::ostringstream error_message;
383  error_message << "The pointer in_distribution_pt[" << dist_i
384  << "] is null.\n";
385  throw OomphLibError(error_message.str(),
386  OOMPH_CURRENT_FUNCTION,
387  OOMPH_EXCEPTION_LOCATION);
388  }
389  }
390
391  /// /////////////////////////////
392
393  // Check that all in distributions are built.
394  for (unsigned dist_i = 0; dist_i < ndistributions; dist_i++)
395  {
396  if (!in_distribution_pt[dist_i]->built())
397  {
398  std::ostringstream error_message;
399  error_message << "The in_distribution_pt[" << dist_i
400  << "] is not built.\n";
401  throw OomphLibError(error_message.str(),
402  OOMPH_CURRENT_FUNCTION,
403  OOMPH_EXCEPTION_LOCATION);
404  }
405  }
406
407  /// /////////////////////////////
408
409  // Check that all communicators to concatenate are the same
410  // by comparing all communicators against the first one.
411  const OomphCommunicator first_comm =
412  *(in_distribution_pt[0]->communicator_pt());
413
414  for (unsigned dist_i = 0; dist_i < ndistributions; dist_i++)
415  {
416  // Get the Communicator for the current vector.
417  const OomphCommunicator another_comm =
418  *(in_distribution_pt[dist_i]->communicator_pt());
419
420  if (!(first_comm == another_comm))
421  {
422  std::ostringstream error_message;
423  error_message << "The communicator in position " << dist_i << " is \n"
424  << "not the same as the first one.\n";
425  throw OomphLibError(error_message.str(),
426  OOMPH_CURRENT_FUNCTION,
427  OOMPH_EXCEPTION_LOCATION);
428  }
429  }
430
431  /// /////////////////////////////
432
433  // Ensure that all distributions are either distributed or not.
434  // This is because we use the distributed() function from the first
435  // distribution to indicate if the result distribution should be
436  // distributed or not.
437  bool first_distributed = in_distribution_pt[0]->distributed();
438  for (unsigned dist_i = 0; dist_i < ndistributions; dist_i++)
439  {
440  // Is the current distribution distributed?
441  bool another_distributed = in_distribution_pt[dist_i]->distributed();
442  if (first_distributed != another_distributed)
443  {
444  std::ostringstream error_message;
445  error_message
446  << "The distribution in position " << dist_i << " has a different\n"
447  << "different distributed boolean than the distribution.\n";
448  throw OomphLibError(error_message.str(),
449  OOMPH_CURRENT_FUNCTION,
450  OOMPH_EXCEPTION_LOCATION);
451  }
452  }
453
454  /// /////////////////////////////
455
456  // Check that the out distribution is not built.
457  if (out_distribution.built())
458  {
459  std::ostringstream error_message;
460  error_message << "The out distribution is built.\n"
461  << "Please clear it.\n";
462  throw OomphLibError(error_message.str(),
463  OOMPH_CURRENT_FUNCTION,
464  OOMPH_EXCEPTION_LOCATION);
465  }
466
467 #endif
468
469  // Get the communicator pointer
470  const OomphCommunicator* const comm_pt =
471  in_distribution_pt[0]->communicator_pt();
472
473  // Number of processors
474  unsigned nproc = comm_pt->nproc();
475
476  // First determine the out_nrow_local and the out_nrow (the global nrow)
477  unsigned out_nrow_local = 0;
478  unsigned out_nrow = 0;
479  for (unsigned in_dist_i = 0; in_dist_i < ndistributions; in_dist_i++)
480  {
481  out_nrow_local += in_distribution_pt[in_dist_i]->nrow_local();
482  out_nrow += in_distribution_pt[in_dist_i]->nrow();
483  }
484
485  // Now calculate the first row for this processor. We need to know the
486  // out_nrow_local for all the other processors, MPI_Allgather must be
487  // used.
488  unsigned out_first_row = 0;
489
490  // Distributed case: We need to gather all the out_nrow_local from all
491  // processors, the out_first_row for this processors is the partial sum up
492  // of the out_nrow_local to my_rank.
493  bool distributed = in_distribution_pt[0]->distributed();
494  if (distributed)
495  {
496 #ifdef OOMPH_HAS_MPI
497  // My rank
498  unsigned my_rank = comm_pt->my_rank();
499
500  unsigned* out_nrow_local_all = new unsigned[nproc];
501  MPI_Allgather(&out_nrow_local,
502  1,
503  MPI_UNSIGNED,
504  out_nrow_local_all,
505  1,
506  MPI_UNSIGNED,
507  comm_pt->mpi_comm());
508  for (unsigned proc_i = 0; proc_i < my_rank; proc_i++)
509  {
510  out_first_row += out_nrow_local_all[proc_i];
511  }
512  delete[] out_nrow_local_all;
513 #endif
514  }
515  // Non-distributed case
516  else
517  {
518  out_first_row = 0;
519  }
520
521  // Build the distribution.
522  if (nproc == 1 || !distributed)
523  {
524  // Some ambiguity here -- on a single processor it doesn't really
525  // matter if we think of ourselves as distributed or not but this
526  // follows the pattern used elsewhere.
527  out_distribution.build(comm_pt, out_nrow, false);
528  }
529  else
530  {
531  out_distribution.build(
532  comm_pt, out_first_row, out_nrow_local, out_nrow);
533  }
534  } // function concatenate
535
536  } // end of namespace LinearAlgebraDistributionHelpers
537
538 } // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed
Vector< unsigned > first_row_vector() const
return the first_row Vector
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero.
bool operator==(const LinearAlgebraDistribution &other_dist) const
== Operator
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
Vector< unsigned > Nrow_local
the number of local rows on the processor
void build(const OomphCommunicator *const comm_pt, const unsigned &first_row, const unsigned &nrow_local, const unsigned &nrow=0)
Sets the distribution. Takes first_row, nrow_local and nrow as arguments. If nrow is not provided or ...
OomphCommunicator * Comm_pt
the pointer to the MPI communicator object in this distribution
unsigned Nrow
the number of global rows
Vector< unsigned > First_row
the first row on this processor
Vector< unsigned > nrow_local_vector() const
return the nrow_local Vector
bool built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
bool Distributed
flag to indicate whether this distribution describes an object that is distributed over the processor...
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned.
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
void concatenate(const Vector< LinearAlgebraDistribution * > &in_distribution_pt, LinearAlgebraDistribution &out_distribution)
Takes a vector of LinearAlgebraDistribution objects and concatenates them such that the nrow_local of...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
std::ostream & operator<<(std::ostream &out, const DoubleVector &v)
output operator