parallel_algebra.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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
58  StandardType(const TypeVector<T> *example=NULL) {
59  // We need an example for MPI_Address to use
60  TypeVector<T> *ex;
61  AutoPtr<TypeVector<T> > temp;
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[LIBMESH_DIM+2];
82  MPI_Aint displs[LIBMESH_DIM+2];
83  MPI_Datatype types[LIBMESH_DIM+2];
84  MPI_Aint start, later;
85 
86  MPI_Address(ex, &start);
87  blocklengths[0] = 1;
88  displs[0] = 0;
89  types[0] = MPI_LB;
90  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
91  {
92  MPI_Address(&((*ex)(i)), &later);
93  blocklengths[i+1] = 1;
94  displs[i+1] = later - start;
95  types[i+1] = T_type;
96  }
97  MPI_Address((ex+1), &later);
98  blocklengths[LIBMESH_DIM+1] = 1;
99  displs[LIBMESH_DIM+1] = later - start;
100  types[LIBMESH_DIM+1] = MPI_UB;
101 
102  MPI_Type_struct (LIBMESH_DIM+2, blocklengths, displs, types, &_static_type);
103 
104 #else // MPI_VERSION >= 2
105 
106  int blocklength = LIBMESH_DIM;
107  MPI_Aint displs, start;
108  MPI_Datatype tmptype, type = T_type;
109 
110  MPI_Get_address (ex, &start);
111  MPI_Get_address (&((*ex)(0)), &displs);
112 
113  // subtract off offset to first value from the beginning of the structure
114  displs -= start;
115 
116  // create a prototype structure
117  MPI_Type_create_struct (1, &blocklength, &displs, &type, &tmptype);
118 
119  // resize the structure type to account for padding, if any
120  MPI_Type_create_resized (tmptype, 0, sizeof(TypeVector<T>), &_static_type);
121 #endif
122 
123  MPI_Type_commit (&_static_type);
124 #endif // #ifdef LIBMESH_HAVE_MPI
125 
126  _is_initialized = true;
127  }
128  _datatype = _static_type;
129  }
130  };
131 
132  template <typename T>
133  class StandardType<VectorValue<T> > : public DataType
134  {
135  public:
136  explicit
137  StandardType(const VectorValue<T> *example=NULL) {
138  // We need an example for MPI_Address to use
139  VectorValue<T> *ex;
140  AutoPtr<VectorValue<T> > temp;
141  if (example)
142  ex = const_cast<VectorValue<T> *>(example);
143  else
144  {
145  temp.reset(new VectorValue<T>());
146  ex = temp.get();
147  }
148 
149  // _static_type never gets freed, but it only gets committed once
150  // per T, so it's not a *huge* memory leak...
151  static data_type _static_type;
152  static bool _is_initialized = false;
153  if (!_is_initialized)
154  {
155 #ifdef LIBMESH_HAVE_MPI
156  StandardType<T> T_type(&((*ex)(0)));
157 
158 #if MPI_VERSION == 1
159 
160  int blocklengths[LIBMESH_DIM+2];
161  MPI_Aint displs[LIBMESH_DIM+2];
162  MPI_Datatype types[LIBMESH_DIM+2];
163  MPI_Aint start, later;
164 
165  MPI_Address(ex, &start);
166  blocklengths[0] = 1;
167  displs[0] = 0;
168  types[0] = MPI_LB;
169  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
170  {
171  MPI_Address(&((*ex)(i)), &later);
172  blocklengths[i+1] = 1;
173  displs[i+1] = later - start;
174  types[i+1] = T_type;
175  }
176  MPI_Address((ex+1), &later);
177  blocklengths[LIBMESH_DIM+1] = 1;
178  displs[LIBMESH_DIM+1] = later - start;
179  types[LIBMESH_DIM+1] = MPI_UB;
180 
181  MPI_Type_struct (LIBMESH_DIM+2, blocklengths, displs, types, &_static_type);
182 
183 #else // MPI_VERSION >= 2
184 
185  int blocklength = LIBMESH_DIM;
186  MPI_Aint displs, start;
187  MPI_Datatype tmptype, type = T_type;
188 
189  MPI_Get_address (ex, &start);
190  MPI_Get_address (&((*ex)(0)), &displs);
191 
192  // subtract off offset to first value from the beginning of the structure
193  displs -= start;
194 
195  // create a prototype structure
196  MPI_Type_create_struct (1, &blocklength, &displs, &type, &tmptype);
197 
198  // resize the structure type to account for padding, if any
199  MPI_Type_create_resized (tmptype, 0, sizeof(VectorValue<T>), &_static_type);
200 #endif
201 
202  MPI_Type_commit (&_static_type);
203 #endif // #ifdef LIBMESH_HAVE_MPI
204 
205  _is_initialized = true;
206  }
207  _datatype = _static_type;
208  }
209  };
210 
211  template <>
212  class StandardType<Point> : public DataType
213  {
214  public:
215  explicit
216  StandardType(const Point *example=NULL) {
217  // We need an example for MPI_Address to use
218  Point *ex;
219  AutoPtr<Point> temp;
220  if (example)
221  ex = const_cast<Point *>(example);
222  else
223  {
224  temp.reset(new Point());
225  ex = temp.get();
226  }
227 
228  // _static_type never gets freed, but it only gets committed once
229  // per T, so it's not a *huge* memory leak...
230  static data_type _static_type;
231  static bool _is_initialized = false;
232  if (!_is_initialized)
233  {
234 #ifdef LIBMESH_HAVE_MPI
235  StandardType<Real> T_type(&((*ex)(0)));
236 
237 #if MPI_VERSION == 1
238 
239  int blocklengths[LIBMESH_DIM+2];
240  MPI_Aint displs[LIBMESH_DIM+2];
241  MPI_Datatype types[LIBMESH_DIM+2];
242  MPI_Aint start, later;
243 
244  MPI_Address(ex, &start);
245  blocklengths[0] = 1;
246  displs[0] = 0;
247  types[0] = MPI_LB;
248  for (unsigned int i=0; i != LIBMESH_DIM; ++i)
249  {
250  MPI_Address(&((*ex)(i)), &later);
251  blocklengths[i+1] = 1;
252  displs[i+1] = later - start;
253  types[i+1] = T_type;
254  }
255  MPI_Address((ex+1), &later);
256  blocklengths[LIBMESH_DIM+1] = 1;
257  displs[LIBMESH_DIM+1] = later - start;
258  types[LIBMESH_DIM+1] = MPI_UB;
259 
260  MPI_Type_struct (LIBMESH_DIM+2, blocklengths, displs, types, &_static_type);
261 
262 #else // MPI_VERSION >= 2
263 
264  int blocklength = LIBMESH_DIM;
265  MPI_Aint displs, start;
266  MPI_Datatype tmptype, type = T_type;
267 
268  MPI_Get_address (ex, &start);
269  MPI_Get_address (&((*ex)(0)), &displs);
270 
271  // subtract off offset to first value from the beginning of the structure
272  displs -= start;
273 
274  // create a prototype structure
275  MPI_Type_create_struct (1, &blocklength, &displs, &type, &tmptype);
276 
277  // resize the structure type to account for padding, if any
278  MPI_Type_create_resized (tmptype, 0, sizeof(Point), &_static_type);
279 #endif
280 
281  MPI_Type_commit (&_static_type);
282 #endif // #ifdef LIBMESH_HAVE_MPI
283 
284  _is_initialized = true;
285  }
286  _datatype = _static_type;
287  }
288  };
289 
290  // StandardType<> specializations to return a derived MPI datatype
291  // to handle communication of LIBMESH_DIM*LIBMESH_DIM-tensors.
292  //
293  // We use a singleton pattern here because a global variable would
294  // have tried to call MPI functions before MPI got initialized.
295  //
296  // We assume contiguous storage here
297  template <typename T>
298  class StandardType<TypeTensor<T> > : public DataType
299  {
300  public:
301  explicit
302  StandardType(const TypeTensor<T> *example=NULL) :
303  DataType(StandardType<T>(example ? &((*example)(0,0)) : NULL), LIBMESH_DIM*LIBMESH_DIM) {}
304 
305  inline ~StandardType() { this->free(); }
306  };
307 
308  template <typename T>
309  class StandardType<TensorValue<T> > : public DataType
310  {
311  public:
312  explicit
313  StandardType(const TensorValue<T> *example=NULL) :
314  DataType(StandardType<T>(example ? &((*example)(0,0)) : NULL), LIBMESH_DIM*LIBMESH_DIM) {}
315 
316  inline ~StandardType() { this->free(); }
317  };
318 } // namespace Parallel
319 } // namespace libMesh
320 
321 #endif // LIBMESH_PARALLEL_ALGEBRA_H

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:06 UTC

Hosted By:
SourceForge.net Logo