libMesh
parallel_algebra.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 #ifndef LIBMESH_PARALLEL_ALGEBRA_H
20 #define LIBMESH_PARALLEL_ALGEBRA_H
21 
22 // This class contains all the functionality for bin sorting
23 // Templated on the type of keys you will be sorting and the
24 // type of iterator you will be using.
25 
26 
27 // Local Includes
28 #include "libmesh/libmesh_config.h"
29 
30 #include "libmesh/auto_ptr.h"
31 #include "libmesh/parallel.h"
32 #include "libmesh/point.h"
33 #include "libmesh/tensor_value.h"
34 #include "libmesh/vector_value.h"
35 
36 // C++ includes
37 #include <cstddef>
38 
39 namespace libMesh {
40 namespace Parallel {
41 // StandardType<> specializations to return a derived MPI datatype
42 // to handle communication of LIBMESH_DIM-vectors.
43 //
44 // We use static variables to minimize the number of MPI datatype
45 // construction calls executed over the course of the program.
46 //
47 // We use a singleton pattern because a global variable would
48 // have tried to call MPI functions before MPI got initialized.
49 //
50 // We use MPI_Create_struct here because our vector classes might
51 // have vptrs, and I'd rather not have the datatype break in those
52 // cases.
53 template <typename T>
54 class StandardType<TypeVector<T>> : public DataType
55 {
56 public:
57  explicit
59  // We need an example for MPI_Address to use
60  TypeVector<T> * ex;
62  if (example)
63  ex = const_cast<TypeVector<T> *>(example);
64  else
65  {
66  temp.reset(new TypeVector<T>());
67  ex = temp.get();
68  }
69 
70  // _static_type never gets freed, but it only gets committed once
71  // per T, so it's not a *huge* memory leak...
72  static data_type _static_type;
73  static bool _is_initialized = false;
74  if (!_is_initialized)
75  {
76 #ifdef LIBMESH_HAVE_MPI
77  StandardType<T> T_type(&((*ex)(0)));
78 
79 #if MPI_VERSION == 1
80 
81  int blocklengths[3] = {1, LIBMESH_DIM, 1};
82  MPI_Aint displs[3];
83  MPI_Datatype types[3] = {MPI_LB, T_type, MPI_UB};
84  MPI_Aint start, later;
85 
86  libmesh_call_mpi
87  (MPI_Address(ex, &start));
88  displs[0] = 0;
89  libmesh_call_mpi
90  (MPI_Address(&((*ex)(0)), &later));
91  displs[1] = later - start;
92  libmesh_call_mpi
93  (MPI_Address((ex+1), &later));
94  displs[2] = later - start;
95 
96  libmesh_call_mpi
97  (MPI_Type_struct (3, blocklengths, displs, types,
98  &_static_type));
99 
100 #else // MPI_VERSION >= 2
101 
102  int blocklength = LIBMESH_DIM;
103  MPI_Aint displs, start;
104  MPI_Datatype tmptype, type = T_type;
105 
106  libmesh_call_mpi
107  (MPI_Get_address (ex, &start));
108  libmesh_call_mpi
109  (MPI_Get_address (&((*ex)(0)), &displs));
110 
111  // subtract off offset to first value from the beginning of the structure
112  displs -= start;
113 
114  // create a prototype structure
115  libmesh_call_mpi
116  (MPI_Type_create_struct (1, &blocklength, &displs, &type,
117  &tmptype));
118  libmesh_call_mpi
119  (MPI_Type_commit (&tmptype));
120 
121  // resize the structure type to account for padding, if any
122  libmesh_call_mpi
123  (MPI_Type_create_resized (tmptype, 0, sizeof(TypeVector<T>),
124  &_static_type));
125 #endif
126 
127  libmesh_call_mpi
128  (MPI_Type_commit (&_static_type));
129 #endif // #ifdef LIBMESH_HAVE_MPI
130 
131  _is_initialized = true;
132  }
133  _datatype = _static_type;
134  }
135 };
136 
137 template <typename T>
139 {
140 public:
141  explicit
143  // We need an example for MPI_Address to use
144  VectorValue<T> * ex;
146  if (example)
147  ex = const_cast<VectorValue<T> *>(example);
148  else
149  {
150  temp.reset(new VectorValue<T>());
151  ex = temp.get();
152  }
153 
154  // _static_type never gets freed, but it only gets committed once
155  // per T, so it's not a *huge* memory leak...
156  static data_type _static_type;
157  static bool _is_initialized = false;
158  if (!_is_initialized)
159  {
160 #ifdef LIBMESH_HAVE_MPI
161  StandardType<T> T_type(&((*ex)(0)));
162 
163 #if MPI_VERSION == 1
164 
165  int blocklengths[3] = {1, LIBMESH_DIM, 1};
166  MPI_Aint displs[3];
167  MPI_Datatype types[3] = {MPI_LB, T_type, MPI_UB};
168  MPI_Aint start, later;
169 
170  libmesh_call_mpi
171  (MPI_Address(ex, &start));
172  displs[0] = 0;
173  libmesh_call_mpi
174  (MPI_Address(&((*ex)(0)), &later));
175  displs[1] = later - start;
176  libmesh_call_mpi
177  (MPI_Address((ex+1), &later));
178  displs[2] = later - start;
179 
180  libmesh_call_mpi
181  (MPI_Type_struct (3, blocklengths, displs, types,
182  &_static_type));
183 
184 #else // MPI_VERSION >= 2
185 
186  int blocklength = LIBMESH_DIM;
187  MPI_Aint displs, start;
188  MPI_Datatype tmptype, type = T_type;
189 
190  libmesh_call_mpi
191  (MPI_Get_address (ex, &start));
192  libmesh_call_mpi
193  (MPI_Get_address (&((*ex)(0)), &displs));
194 
195  // subtract off offset to first value from the beginning of the structure
196  displs -= start;
197 
198  // create a prototype structure
199  libmesh_call_mpi
200  (MPI_Type_create_struct (1, &blocklength, &displs, &type,
201  &tmptype));
202  libmesh_call_mpi
203  (MPI_Type_commit (&tmptype));
204 
205  // resize the structure type to account for padding, if any
206  libmesh_call_mpi
207  (MPI_Type_create_resized (tmptype, 0,
208  sizeof(VectorValue<T>),
209  &_static_type));
210 #endif
211 
212  libmesh_call_mpi
213  (MPI_Type_commit (&_static_type));
214 #endif // #ifdef LIBMESH_HAVE_MPI
215 
216  _is_initialized = true;
217  }
218  _datatype = _static_type;
219  }
220 };
221 
222 template <>
223 class StandardType<Point> : public DataType
224 {
225 public:
226  explicit
228  {
229  // Prevent unused variable warnings when !LIBMESH_HAVE_MPI
230  libmesh_ignore(example);
231 
232  // _static_type never gets freed, but it only gets committed once
233  // per T, so it's not a *huge* memory leak...
234  static data_type _static_type;
235  static bool _is_initialized = false;
236  if (!_is_initialized)
237  {
238 #ifdef LIBMESH_HAVE_MPI
239 
240  // We need an example for MPI_Address to use
241  Point * ex;
242 
243  UniquePtr<Point> temp;
244  if (example)
245  ex = const_cast<Point *>(example);
246  else
247  {
248  temp.reset(new Point());
249  ex = temp.get();
250  }
251 
252  StandardType<Real> T_type(&((*ex)(0)));
253 
254 #if MPI_VERSION == 1
255 
256  int blocklengths[3] = {1, LIBMESH_DIM, 1};
257  MPI_Aint displs[3];
258  MPI_Datatype types[3] = {MPI_LB, T_type, MPI_UB};
259  MPI_Aint start, later;
260 
261  libmesh_call_mpi
262  (MPI_Address(ex, &start));
263  displs[0] = 0;
264  libmesh_call_mpi
265  (MPI_Address(&((*ex)(0)), &later));
266  displs[1] = later - start;
267  libmesh_call_mpi
268  (MPI_Address((ex+1), &later));
269  displs[2] = later - start;
270 
271  libmesh_call_mpi
272  (MPI_Type_struct (3, blocklengths, displs, types,
273  &_static_type));
274 
275 #else // MPI_VERSION >= 2
276 
277  int blocklength = LIBMESH_DIM;
278  MPI_Aint displs, start;
279  MPI_Datatype tmptype, type = T_type;
280 
281  libmesh_call_mpi
282  (MPI_Get_address (ex, &start));
283  libmesh_call_mpi
284  (MPI_Get_address (&((*ex)(0)), &displs));
285 
286  // subtract off offset to first value from the beginning of the structure
287  displs -= start;
288 
289  // create a prototype structure
290  libmesh_call_mpi
291  (MPI_Type_create_struct (1, &blocklength, &displs, &type,
292  &tmptype));
293  libmesh_call_mpi
294  (MPI_Type_commit (&tmptype));
295 
296  // resize the structure type to account for padding, if any
297  libmesh_call_mpi
298  (MPI_Type_create_resized (tmptype, 0, sizeof(Point),
299  &_static_type));
300 #endif
301 
302  libmesh_call_mpi
303  (MPI_Type_commit (&_static_type));
304 #endif // #ifdef LIBMESH_HAVE_MPI
305 
306  _is_initialized = true;
307  }
308  _datatype = _static_type;
309  }
310 };
311 
312 // OpFunction<> specializations to return an MPI_Op version of the
313 // reduction operations on LIBMESH_DIM-vectors.
314 //
315 // We use static variables to minimize the number of MPI datatype
316 // construction calls executed over the course of the program.
317 //
318 // We use a singleton pattern because a global variable would
319 // have tried to call MPI functions before MPI got initialized.
320 //
321 // min() and max() are applied component-wise; this makes them useful
322 // for bounding box reduction operations.
323 template <typename V>
325 {
326 public:
327 #ifdef LIBMESH_HAVE_MPI
328  static void vector_max (void * invec, void * inoutvec, int * len, MPI_Datatype *)
329  {
330  V *in = static_cast<V *>(invec);
331  V *inout = static_cast<V *>(inoutvec);
332  for (int i=0; i != *len; ++i)
333  for (int d=0; d != LIBMESH_DIM; ++d)
334  inout[i](d) = std::max(in[i](d), inout[i](d));
335  }
336 
337  static void vector_min (void * invec, void * inoutvec, int * len, MPI_Datatype *)
338  {
339  V *in = static_cast<V *>(invec);
340  V *inout = static_cast<V *>(inoutvec);
341  for (int i=0; i != *len; ++i)
342  for (int d=0; d != LIBMESH_DIM; ++d)
343  inout[i](d) = std::min(in[i](d), inout[i](d));
344  }
345 
346  static void vector_sum (void * invec, void * inoutvec, int * len, MPI_Datatype *)
347  {
348  V *in = static_cast<V *>(invec);
349  V *inout = static_cast<V *>(inoutvec);
350  for (int i=0; i != *len; ++i)
351  for (int d=0; d != LIBMESH_DIM; ++d)
352  inout[i](d) += in[i](d);
353  }
354 
355  static MPI_Op max()
356  {
357  // _static_op never gets freed, but it only gets committed once
358  // per T, so it's not a *huge* memory leak...
359  static MPI_Op _static_op;
360  static bool _is_initialized = false;
361  if (!_is_initialized)
362  {
363  libmesh_call_mpi
364  (MPI_Op_create (vector_max, /*commute=*/ true,
365  &_static_op));
366 
367  _is_initialized = true;
368  }
369 
370  return _static_op;
371  }
372  static MPI_Op min()
373  {
374  // _static_op never gets freed, but it only gets committed once
375  // per T, so it's not a *huge* memory leak...
376  static MPI_Op _static_op;
377  static bool _is_initialized = false;
378  if (!_is_initialized)
379  {
380  libmesh_call_mpi
381  (MPI_Op_create (vector_min, /*commute=*/ true,
382  &_static_op));
383 
384  _is_initialized = true;
385  }
386 
387  return _static_op;
388  }
389  static MPI_Op sum()
390  {
391  // _static_op never gets freed, but it only gets committed once
392  // per T, so it's not a *huge* memory leak...
393  static MPI_Op _static_op;
394  static bool _is_initialized = false;
395  if (!_is_initialized)
396  {
397  libmesh_call_mpi
398  (MPI_Op_create (vector_sum, /*commute=*/ true,
399  &_static_op));
400 
401  _is_initialized = true;
402  }
403 
404  return _static_op;
405  }
406 
407 #endif // LIBMESH_HAVE_MPI
408 };
409 
410 template <typename T>
411 class OpFunction<TypeVector<T>> : public TypeVectorOpFunction<TypeVector<T>> {};
412 
413 template <typename T>
414 class OpFunction<VectorValue<T>> : public TypeVectorOpFunction<VectorValue<T>> {};
415 
416 template <>
417 class OpFunction<Point> : public TypeVectorOpFunction<Point> {};
418 
419 // StandardType<> specializations to return a derived MPI datatype
420 // to handle communication of LIBMESH_DIM*LIBMESH_DIM-tensors.
421 //
422 // We assume contiguous storage here
423 template <typename T>
424 class StandardType<TypeTensor<T>> : public DataType
425 {
426 public:
427  explicit
429  DataType(StandardType<T>(example ? &((*example)(0,0)) : libmesh_nullptr), LIBMESH_DIM*LIBMESH_DIM) {}
430 
431  inline ~StandardType() { this->free(); }
432 };
433 
434 template <typename T>
436 {
437 public:
438  explicit
440  DataType(StandardType<T>(example ? &((*example)(0,0)) : libmesh_nullptr), LIBMESH_DIM*LIBMESH_DIM) {}
441 
442  inline ~StandardType() { this->free(); }
443 };
444 } // namespace Parallel
445 } // namespace libMesh
446 
447 #endif // LIBMESH_PARALLEL_ALGEBRA_H
static void vector_max(void *invec, void *inoutvec, int *len, MPI_Datatype *)
StandardType(const TensorValue< T > *example=libmesh_nullptr)
StandardType(const VectorValue< T > *example=libmesh_nullptr)
const class libmesh_nullptr_t libmesh_nullptr
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
static void vector_min(void *invec, void *inoutvec, int *len, MPI_Datatype *)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
This class defines a tensor in LIBMESH_DIM dimensional space of type T.
Definition: tensor_tools.h:32
StandardType(const Point *example=libmesh_nullptr)
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:255
StandardType(const TypeTensor< T > *example=libmesh_nullptr)
Templated class to provide the appropriate MPI datatype for use with built-in C types or simple C++ c...
Definition: parallel.h:380
This class defines a vector in LIBMESH_DIM dimensional space of type T.
Definition: tensor_tools.h:30
Encapsulates the MPI_Datatype.
Definition: parallel.h:289
void libmesh_ignore(const T &)
StandardType(const TypeVector< T > *example=libmesh_nullptr)
Templated class to provide the appropriate MPI reduction operations for use with built-in C types or ...
Definition: parallel.h:406
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static void vector_sum(void *invec, void *inoutvec, int *len, MPI_Datatype *)
This class defines a tensor in LIBMESH_DIM dimensional Real or Complex space.