Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 :
23 : // libMesh includes
24 : #include "libmesh/libmesh_config.h"
25 : #include "libmesh/point.h"
26 : #include "libmesh/tensor_value.h"
27 : #include "libmesh/vector_value.h"
28 :
29 : // TIMPI includes
30 : #include "timpi/attributes.h"
31 : #include "timpi/op_function.h"
32 : #include "timpi/standard_type.h"
33 :
34 : // C++ includes
35 : #include <cstddef>
36 : #include <memory>
37 : #include <type_traits>
38 :
39 : namespace libMesh {
40 :
41 : // OpFunction<> specializations to return an MPI_Op version of the
42 : // reduction operations on 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 : // min() and max() are applied component-wise; this makes them useful
51 : // for bounding box reduction operations.
52 :
53 : template <typename V>
54 : class TypeVectorOpFunction
55 : {
56 : public:
57 : #ifdef LIBMESH_HAVE_MPI
58 2898326 : static void vector_max (void * invec, void * inoutvec, int * len, MPI_Datatype *)
59 : {
60 465906 : V *in = static_cast<V *>(invec);
61 465906 : V *inout = static_cast<V *>(inoutvec);
62 5796652 : for (int i=0; i != *len; ++i)
63 11593304 : for (int d=0; d != LIBMESH_DIM; ++d)
64 10107525 : inout[i](d) = std::max(in[i](d), inout[i](d));
65 2898326 : }
66 :
67 2898464 : static void vector_min (void * invec, void * inoutvec, int * len, MPI_Datatype *)
68 : {
69 465906 : V *in = static_cast<V *>(invec);
70 465906 : V *inout = static_cast<V *>(inoutvec);
71 5796928 : for (int i=0; i != *len; ++i)
72 11593856 : for (int d=0; d != LIBMESH_DIM; ++d)
73 10045680 : inout[i](d) = std::min(in[i](d), inout[i](d));
74 2898464 : }
75 :
76 : static void vector_sum (void * invec, void * inoutvec, int * len, MPI_Datatype *)
77 : {
78 : V *in = static_cast<V *>(invec);
79 : V *inout = static_cast<V *>(inoutvec);
80 : for (int i=0; i != *len; ++i)
81 : for (int d=0; d != LIBMESH_DIM; ++d)
82 : inout[i](d) += in[i](d);
83 : }
84 :
85 2069326 : static MPI_Op max()
86 : {
87 : // _static_op never gets freed, but it only gets committed once
88 : // per T, so it's not a *huge* memory leak...
89 : static MPI_Op _static_op;
90 : static bool _is_initialized = false;
91 2085938 : if (!_is_initialized)
92 : {
93 13410 : timpi_call_mpi
94 : (MPI_Op_create (vector_max, /*commute=*/ true,
95 : &_static_op));
96 :
97 13410 : _is_initialized = true;
98 : }
99 :
100 2085938 : return _static_op;
101 : }
102 2069326 : static MPI_Op min()
103 : {
104 : // _static_op never gets freed, but it only gets committed once
105 : // per T, so it's not a *huge* memory leak...
106 : static MPI_Op _static_op;
107 : static bool _is_initialized = false;
108 2085938 : if (!_is_initialized)
109 : {
110 13410 : timpi_call_mpi
111 : (MPI_Op_create (vector_min, /*commute=*/ true,
112 : &_static_op));
113 :
114 13410 : _is_initialized = true;
115 : }
116 :
117 2085938 : return _static_op;
118 : }
119 : static MPI_Op sum()
120 : {
121 : // _static_op never gets freed, but it only gets committed once
122 : // per T, so it's not a *huge* memory leak...
123 : static MPI_Op _static_op;
124 : static bool _is_initialized = false;
125 : if (!_is_initialized)
126 : {
127 : timpi_call_mpi
128 : (MPI_Op_create (vector_sum, /*commute=*/ true,
129 : &_static_op));
130 :
131 : _is_initialized = true;
132 : }
133 :
134 : return _static_op;
135 : }
136 :
137 : #endif // LIBMESH_HAVE_MPI
138 : };
139 :
140 :
141 : template <typename V>
142 : struct TypeVectorAttributes
143 : {
144 : static const bool has_min_max = true;
145 0 : static void set_lowest(V & x) {
146 0 : for (int d=0; d != LIBMESH_DIM; ++d)
147 0 : TIMPI::Attributes<typename std::remove_reference<decltype(x(d))>::type>::set_lowest(x(d));
148 0 : }
149 0 : static void set_highest(V & x) {
150 0 : for (int d=0; d != LIBMESH_DIM; ++d)
151 0 : TIMPI::Attributes<typename std::remove_reference<decltype(x(d))>::type>::set_highest(x(d));
152 0 : }
153 : };
154 :
155 : } // namespace libMesh
156 :
157 :
158 : namespace TIMPI {
159 :
160 : // StandardType<> specializations to return a derived MPI datatype
161 : // to handle communication of LIBMESH_DIM-vectors.
162 : //
163 : // We use MPI_Create_struct here because our vector classes might
164 : // have vptrs, and I'd rather not have the datatype break in those
165 : // cases.
166 : template <typename T>
167 : class StandardType<libMesh::TypeVector<T>, typename std::enable_if<StandardType<T>::is_fixed_type>::type>
168 : : public DataType
169 : {
170 : public:
171 : explicit
172 : StandardType(const libMesh::TypeVector<T> * example=nullptr)
173 : {
174 : using libMesh::TypeVector;
175 :
176 : // We need an example for MPI_Address to use
177 : TypeVector<T> * ex;
178 : std::unique_ptr<TypeVector<T>> temp;
179 : if (example)
180 : ex = const_cast<TypeVector<T> *>(example);
181 : else
182 : {
183 : temp = std::make_unique<TypeVector<T>>();
184 : ex = temp.get();
185 : }
186 :
187 : #ifdef LIBMESH_HAVE_MPI
188 : StandardType<T> T_type(&((*ex)(0)));
189 :
190 : // We require MPI-2 here:
191 : int blocklength = LIBMESH_DIM;
192 : MPI_Aint displs, start;
193 : MPI_Datatype tmptype, type = T_type;
194 :
195 : timpi_call_mpi
196 : (MPI_Get_address (ex, &start));
197 : timpi_call_mpi
198 : (MPI_Get_address (&((*ex)(0)), &displs));
199 :
200 : // subtract off offset to first value from the beginning of the structure
201 : displs -= start;
202 :
203 : // create a prototype structure
204 : timpi_call_mpi
205 : (MPI_Type_create_struct (1, &blocklength, &displs, &type,
206 : &tmptype));
207 : timpi_call_mpi
208 : (MPI_Type_commit (&tmptype));
209 :
210 : // resize the structure type to account for padding, if any
211 : timpi_call_mpi
212 : (MPI_Type_create_resized (tmptype, 0, sizeof(TypeVector<T>),
213 : &_datatype));
214 :
215 : timpi_call_mpi
216 : (MPI_Type_commit (&_datatype));
217 :
218 : timpi_call_mpi
219 : (MPI_Type_free (&tmptype));
220 : #endif // #ifdef LIBMESH_HAVE_MPI
221 : }
222 :
223 : StandardType(const StandardType<libMesh::TypeVector<T>> & timpi_mpi_var(t))
224 : : DataType()
225 : {
226 : timpi_call_mpi (MPI_Type_dup (t._datatype, &_datatype));
227 : }
228 :
229 : ~StandardType() { this->free(); }
230 :
231 : static const bool is_fixed_type = true;
232 : };
233 :
234 :
235 : template <typename T>
236 : class StandardType<libMesh::VectorValue<T>, typename std::enable_if<StandardType<T>::is_fixed_type>::type>
237 : : public DataType
238 : {
239 : public:
240 : explicit
241 0 : StandardType(const libMesh::VectorValue<T> * example=nullptr)
242 0 : {
243 : using libMesh::VectorValue;
244 :
245 : // We need an example for MPI_Address to use
246 : VectorValue<T> * ex;
247 0 : std::unique_ptr<VectorValue<T>> temp;
248 0 : if (example)
249 0 : ex = const_cast<VectorValue<T> *>(example);
250 : else
251 : {
252 0 : temp = std::make_unique<VectorValue<T>>();
253 0 : ex = temp.get();
254 : }
255 :
256 : #ifdef LIBMESH_HAVE_MPI
257 0 : StandardType<T> T_type(&((*ex)(0)));
258 :
259 0 : int blocklength = LIBMESH_DIM;
260 : MPI_Aint displs, start;
261 0 : MPI_Datatype tmptype, type = T_type;
262 :
263 0 : timpi_call_mpi
264 : (MPI_Get_address (ex, &start));
265 0 : timpi_call_mpi
266 : (MPI_Get_address (&((*ex)(0)), &displs));
267 :
268 : // subtract off offset to first value from the beginning of the structure
269 0 : displs -= start;
270 :
271 : // create a prototype structure
272 0 : timpi_call_mpi
273 : (MPI_Type_create_struct (1, &blocklength, &displs, &type,
274 : &tmptype));
275 0 : timpi_call_mpi
276 : (MPI_Type_commit (&tmptype));
277 :
278 : // resize the structure type to account for padding, if any
279 0 : timpi_call_mpi
280 : (MPI_Type_create_resized (tmptype, 0,
281 : sizeof(VectorValue<T>),
282 : &_datatype));
283 :
284 0 : timpi_call_mpi
285 : (MPI_Type_commit (&_datatype));
286 :
287 0 : timpi_call_mpi
288 : (MPI_Type_free (&tmptype));
289 : #endif // #ifdef LIBMESH_HAVE_MPI
290 0 : }
291 :
292 : StandardType(const StandardType<libMesh::VectorValue<T>> & timpi_mpi_var(t))
293 : : DataType()
294 : {
295 : #ifdef LIBMESH_HAVE_MPI
296 : timpi_call_mpi (MPI_Type_dup (t._datatype, &_datatype));
297 : #endif
298 : }
299 :
300 0 : ~StandardType() { this->free(); }
301 :
302 : static const bool is_fixed_type = true;
303 : };
304 :
305 : template <>
306 : class StandardType<libMesh::Point> : public DataType
307 : {
308 : public:
309 : explicit
310 2202408 : StandardType(const libMesh::Point * example=nullptr)
311 2202408 : {
312 : using libMesh::Point;
313 :
314 : // Prevent unused variable warnings when !LIBMESH_HAVE_MPI
315 1835524 : libmesh_ignore(example);
316 :
317 : #ifdef LIBMESH_HAVE_MPI
318 :
319 : // We need an example for MPI_Address to use
320 : Point * ex;
321 :
322 2161894 : std::unique_ptr<Point> temp;
323 2202408 : if (example)
324 1828494 : ex = const_cast<Point *>(example);
325 : else
326 : {
327 12230 : temp = std::make_unique<Point>();
328 7030 : ex = temp.get();
329 : }
330 :
331 3671048 : StandardType<libMesh::Real> T_type(&((*ex)(0)));
332 :
333 2202408 : int blocklength = LIBMESH_DIM;
334 : MPI_Aint displs, start;
335 2202408 : MPI_Datatype tmptype, type = T_type;
336 :
337 2202408 : timpi_call_mpi
338 : (MPI_Get_address (ex, &start));
339 2202408 : timpi_call_mpi
340 : (MPI_Get_address (&((*ex)(0)), &displs));
341 :
342 : // subtract off offset to first value from the beginning of the structure
343 2202408 : displs -= start;
344 :
345 : // create a prototype structure
346 2202408 : timpi_call_mpi
347 : (MPI_Type_create_struct (1, &blocklength, &displs, &type,
348 : &tmptype));
349 2202408 : timpi_call_mpi
350 : (MPI_Type_commit (&tmptype));
351 :
352 : // resize the structure type to account for padding, if any
353 2202408 : timpi_call_mpi
354 : (MPI_Type_create_resized (tmptype, 0, sizeof(Point),
355 : &_datatype));
356 :
357 2202408 : timpi_call_mpi
358 : (MPI_Type_commit (&_datatype));
359 :
360 2202408 : timpi_call_mpi
361 : (MPI_Type_free (&tmptype));
362 : #endif // #ifdef LIBMESH_HAVE_MPI
363 2202408 : }
364 :
365 : StandardType(const StandardType<libMesh::Point> & timpi_mpi_var(t))
366 : : DataType()
367 : {
368 : timpi_call_mpi (MPI_Type_dup (t._datatype, &_datatype));
369 : }
370 :
371 4174950 : ~StandardType() { this->free(); }
372 :
373 : static const bool is_fixed_type = true;
374 : };
375 :
376 : template <typename T>
377 : class OpFunction<libMesh::TypeVector<T>> : public libMesh::TypeVectorOpFunction<libMesh::TypeVector<T>> {};
378 :
379 : template <typename T>
380 : class OpFunction<libMesh::VectorValue<T>> : public libMesh::TypeVectorOpFunction<libMesh::VectorValue<T>> {};
381 :
382 : template <>
383 : class OpFunction<libMesh::Point> : public libMesh::TypeVectorOpFunction<libMesh::Point> {};
384 :
385 :
386 : template <typename T>
387 : class Attributes<libMesh::TypeVector<T>> : public libMesh::TypeVectorAttributes<libMesh::TypeVector<T>> {};
388 :
389 : template <typename T>
390 : class Attributes<libMesh::VectorValue<T>> : public libMesh::TypeVectorAttributes<libMesh::VectorValue<T>> {};
391 :
392 : template <>
393 : class Attributes<libMesh::Point> : public libMesh::TypeVectorAttributes<libMesh::Point> {};
394 :
395 :
396 : // StandardType<> specializations to return a derived MPI datatype
397 : // to handle communication of LIBMESH_DIM*LIBMESH_DIM-tensors.
398 : //
399 : // We assume contiguous storage here
400 : template <typename T>
401 : class StandardType<libMesh::TypeTensor<T>, typename std::enable_if<StandardType<T>::is_fixed_type>::type>
402 : : public DataType
403 : {
404 : public:
405 : explicit
406 : StandardType(const libMesh::TypeTensor<T> * example=nullptr) :
407 : DataType(StandardType<T>(example ? &((*example)(0,0)) : nullptr), LIBMESH_DIM*LIBMESH_DIM) {}
408 :
409 : inline ~StandardType() { this->free(); }
410 :
411 : static const bool is_fixed_type = true;
412 : };
413 :
414 : template <typename T>
415 : class StandardType<libMesh::TensorValue<T>, typename std::enable_if<StandardType<T>::is_fixed_type>::type>
416 : : public DataType
417 : {
418 : public:
419 : explicit
420 0 : StandardType(const libMesh::TensorValue<T> * example=nullptr) :
421 0 : DataType(StandardType<T>(example ? &((*example)(0,0)) : nullptr), LIBMESH_DIM*LIBMESH_DIM) {}
422 :
423 0 : inline ~StandardType() { this->free(); }
424 :
425 : static const bool is_fixed_type = true;
426 : };
427 : } // namespace TIMPI
428 :
429 : #endif // LIBMESH_PARALLEL_ALGEBRA_H
|
