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 :
20 : #include "libmesh/libmesh_common.h"
21 :
22 :
23 :
24 : #ifndef LIBMESH_DISTRIBUTED_VECTOR_H
25 : #define LIBMESH_DISTRIBUTED_VECTOR_H
26 :
27 : // Local includes
28 : #include "libmesh/int_range.h"
29 : #include "libmesh/numeric_vector.h"
30 : #include "libmesh/parallel.h"
31 :
32 : // C++ includes
33 : #include <vector>
34 : #include <algorithm>
35 : #include <limits>
36 : #include <mutex>
37 :
38 : namespace libMesh
39 : {
40 :
41 : /**
42 : * This class provides a simple parallel, distributed vector datatype
43 : * which is specific to libmesh. Offers some collective communication
44 : * capabilities.
45 : *
46 : * \note The class will sill function without MPI, but only on one
47 : * processor. All overridden virtual functions are documented in
48 : * numeric_vector.h.
49 : *
50 : * \author Benjamin S. Kirk
51 : * \date 2003
52 : */
53 : template <typename T>
54 : class DistributedVector final : public NumericVector<T>
55 : {
56 : public:
57 :
58 : /**
59 : * Dummy-Constructor. Dimension=0
60 : */
61 : explicit
62 : DistributedVector (const Parallel::Communicator & comm,
63 : const ParallelType = AUTOMATIC);
64 :
65 : /**
66 : * Constructor. Set dimension to \p n and initialize all elements with zero.
67 : */
68 : explicit
69 : DistributedVector (const Parallel::Communicator & comm,
70 : const numeric_index_type n,
71 : const ParallelType ptype = AUTOMATIC);
72 :
73 : /**
74 : * Constructor. Set local dimension to \p n_local, the global dimension
75 : * to \p n, and initialize all elements with zero.
76 : */
77 : DistributedVector (const Parallel::Communicator & comm,
78 : const numeric_index_type n,
79 : const numeric_index_type n_local,
80 : const ParallelType ptype = AUTOMATIC);
81 :
82 : /**
83 : * Constructor. Set local dimension to \p n_local, the global
84 : * dimension to \p n, but additionally reserve memory for the
85 : * indices specified by the \p ghost argument.
86 : */
87 : DistributedVector (const Parallel::Communicator & comm,
88 : const numeric_index_type N,
89 : const numeric_index_type n_local,
90 : const std::vector<numeric_index_type> & ghost,
91 : const ParallelType ptype = AUTOMATIC);
92 :
93 : /**
94 : * Copy assignment operator. We cannot default this (although it
95 : * essentially implements the default behavior) because the
96 : * compiler-generated default attempts to automatically call the
97 : * base class (NumericVector) copy assignment operator, which we
98 : * have chosen to make pure virtual for other design reasons.
99 : */
100 : DistributedVector & operator= (const DistributedVector &);
101 :
102 : /**
103 : * The 5 special functions can be defaulted for this class, as it
104 : * does not manage any memory itself.
105 : */
106 : DistributedVector (DistributedVector &&) = default;
107 : DistributedVector (const DistributedVector &) = default;
108 : DistributedVector & operator= (DistributedVector &&) = default;
109 0 : virtual ~DistributedVector () = default;
110 :
111 : virtual void close () override;
112 :
113 : virtual void clear () override;
114 :
115 : virtual void zero () override;
116 :
117 : virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
118 :
119 : virtual std::unique_ptr<NumericVector<T>> clone () const override;
120 :
121 : virtual void init (const numeric_index_type N,
122 : const numeric_index_type n_local,
123 : const bool fast=false,
124 : const ParallelType ptype=AUTOMATIC) override;
125 :
126 : virtual void init (const numeric_index_type N,
127 : const bool fast=false,
128 : const ParallelType ptype=AUTOMATIC) override;
129 :
130 : virtual void init (const numeric_index_type N,
131 : const numeric_index_type n_local,
132 : const std::vector<numeric_index_type> & ghost,
133 : const bool fast = false,
134 : const ParallelType = AUTOMATIC) override;
135 :
136 : virtual void init (const NumericVector<T> & other,
137 : const bool fast = false) override;
138 :
139 : virtual NumericVector<T> & operator= (const T s) override;
140 :
141 : virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
142 :
143 : virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
144 :
145 : virtual Real min () const override;
146 :
147 : virtual Real max () const override;
148 :
149 : virtual T sum() const override;
150 :
151 : virtual Real l1_norm () const override;
152 :
153 : virtual Real l2_norm () const override;
154 :
155 : virtual Real linfty_norm () const override;
156 :
157 : virtual numeric_index_type size () const override;
158 :
159 : virtual numeric_index_type local_size() const override;
160 :
161 : virtual numeric_index_type first_local_index() const override;
162 :
163 : virtual numeric_index_type last_local_index() const override;
164 :
165 : virtual T operator() (const numeric_index_type i) const override;
166 :
167 : virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
168 :
169 : virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
170 :
171 : virtual NumericVector<T> & operator *= (const NumericVector<T> & v) override;
172 :
173 : virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
174 :
175 : virtual void reciprocal() override;
176 :
177 : virtual void conjugate() override;
178 :
179 : virtual void set (const numeric_index_type i, const T value) override;
180 :
181 : virtual void add (const numeric_index_type i, const T value) override;
182 :
183 : virtual void add (const T s) override;
184 :
185 : virtual void add (const NumericVector<T> & V) override;
186 :
187 : virtual void add (const T a, const NumericVector<T> & v) override;
188 :
189 : /**
190 : * We override one NumericVector<T>::add_vector() method but don't
191 : * want to hide the other defaults.
192 : */
193 : using NumericVector<T>::add_vector;
194 :
195 0 : virtual void add_vector (const NumericVector<T> &,
196 : const SparseMatrix<T> &) override
197 0 : { libmesh_not_implemented(); }
198 :
199 0 : virtual void add_vector_transpose (const NumericVector<T> &,
200 : const SparseMatrix<T> &) override
201 0 : { libmesh_not_implemented(); }
202 :
203 : virtual void scale (const T factor) override;
204 :
205 : virtual void abs() override;
206 :
207 : virtual T dot(const NumericVector<T> & V) const override;
208 :
209 : virtual void localize (std::vector<T> & v_local) const override;
210 :
211 : virtual void localize (NumericVector<T> & v_local) const override;
212 :
213 : virtual void localize (NumericVector<T> & v_local,
214 : const std::vector<numeric_index_type> & send_list) const override;
215 :
216 : virtual void localize (std::vector<T> & v_local,
217 : const std::vector<numeric_index_type> & indices) const override;
218 :
219 : virtual void localize (const numeric_index_type first_local_idx,
220 : const numeric_index_type last_local_idx,
221 : const std::vector<numeric_index_type> & send_list) override;
222 :
223 : virtual void localize_to_one (std::vector<T> & v_local,
224 : const processor_id_type proc_id=0) const override;
225 :
226 : virtual void pointwise_mult (const NumericVector<T> & vec1,
227 : const NumericVector<T> & vec2) override;
228 :
229 : virtual void pointwise_divide (const NumericVector<T> & vec1,
230 : const NumericVector<T> & vec2) override;
231 :
232 : virtual void swap (NumericVector<T> & v) override;
233 :
234 : virtual std::size_t max_allowed_id() const override;
235 :
236 : private:
237 :
238 : /**
239 : * Actual vector datatype to hold vector entries.
240 : */
241 : std::vector<T> _values;
242 :
243 : /**
244 : * Entries to add or set on remote processors during the next
245 : * close()
246 : */
247 : std::vector<std::pair<numeric_index_type, T>> _remote_values;
248 :
249 : /**
250 : * Whether we are adding or setting remote values or neither - this
251 : * determines behavior at the next close();
252 : */
253 : enum UnclosedState { DO_NOTHING = 0,
254 : ADD_VALUES,
255 : SET_VALUES,
256 : INVALID_UNCLOSEDSTATE };
257 :
258 : /**
259 : * The current state of this vector, if it is unclosed.
260 : */
261 : UnclosedState _unclosed_state;
262 :
263 : /**
264 : * The global vector size.
265 : */
266 : numeric_index_type _global_size;
267 :
268 : /**
269 : * The local vector size.
270 : */
271 : numeric_index_type _local_size;
272 :
273 : /**
274 : * The first component stored locally.
275 : */
276 : numeric_index_type _first_local_index;
277 :
278 : /**
279 : * The last component (+1) stored locally.
280 : */
281 : numeric_index_type _last_local_index;
282 : };
283 :
284 :
285 : //--------------------------------------------------------------------------
286 : // DistributedVector inline methods
287 : template <typename T>
288 : inline
289 0 : DistributedVector<T>::DistributedVector (const Parallel::Communicator & comm_in,
290 : const ParallelType ptype) :
291 : NumericVector<T>(comm_in, ptype),
292 0 : _unclosed_state (DO_NOTHING),
293 0 : _global_size (0),
294 0 : _local_size (0),
295 0 : _first_local_index(0),
296 0 : _last_local_index (0)
297 : {
298 0 : this->_type = ptype;
299 0 : }
300 :
301 :
302 :
303 : template <typename T>
304 : inline
305 0 : DistributedVector<T>::DistributedVector (const Parallel::Communicator & comm_in,
306 : const numeric_index_type n,
307 : const ParallelType ptype)
308 0 : : NumericVector<T>(comm_in, ptype)
309 : {
310 0 : this->init(n, n, false, ptype);
311 0 : }
312 :
313 :
314 :
315 : template <typename T>
316 : inline
317 0 : DistributedVector<T>::DistributedVector (const Parallel::Communicator & comm_in,
318 : const numeric_index_type n,
319 : const numeric_index_type n_local,
320 : const ParallelType ptype)
321 0 : : NumericVector<T>(comm_in, ptype)
322 : {
323 0 : this->init(n, n_local, false, ptype);
324 0 : }
325 :
326 :
327 :
328 : template <typename T>
329 : inline
330 0 : DistributedVector<T>::DistributedVector (const Parallel::Communicator & comm_in,
331 : const numeric_index_type n,
332 : const numeric_index_type n_local,
333 : const std::vector<numeric_index_type> & ghost,
334 : const ParallelType ptype)
335 0 : : NumericVector<T>(comm_in, ptype)
336 : {
337 0 : this->init(n, n_local, ghost, false, ptype);
338 0 : }
339 :
340 :
341 :
342 : template <typename T>
343 : inline
344 0 : void DistributedVector<T>::init (const numeric_index_type n,
345 : const numeric_index_type n_local,
346 : const bool fast,
347 : const ParallelType ptype)
348 : {
349 : // This function must be run on all processors at once
350 0 : parallel_object_only();
351 :
352 0 : libmesh_assert_less_equal (n_local, n);
353 :
354 0 : if (ptype == AUTOMATIC)
355 : {
356 0 : if (n == n_local)
357 0 : this->_type = SERIAL;
358 : else
359 0 : this->_type = PARALLEL;
360 : }
361 0 : else if (ptype == GHOSTED &&
362 0 : n == n_local) // We can support GHOSTED with no ghosts...
363 0 : this->_type = SERIAL;
364 : else
365 0 : this->_type = ptype;
366 :
367 0 : libmesh_assert ((this->_type==SERIAL && n==n_local) ||
368 : this->_type==PARALLEL);
369 :
370 : // Clear the data structures if already initialized
371 0 : if (this->initialized())
372 0 : this->clear();
373 :
374 : // Initialize data structures
375 0 : _unclosed_state = DO_NOTHING;
376 0 : _values.resize(n_local);
377 0 : _remote_values.clear();
378 0 : _local_size = n_local;
379 0 : _global_size = n;
380 :
381 0 : _first_local_index = 0;
382 :
383 : #ifdef LIBMESH_HAVE_MPI
384 :
385 0 : std::vector<numeric_index_type> local_sizes (this->n_processors(), 0);
386 :
387 0 : local_sizes[this->processor_id()] = n_local;
388 :
389 0 : this->comm().sum(local_sizes);
390 :
391 : // _first_local_index is the sum of _local_size
392 : // for all processor ids less than ours
393 0 : for (auto p : make_range(this->processor_id()))
394 0 : _first_local_index += local_sizes[p];
395 :
396 :
397 : # ifdef DEBUG
398 : // Make sure all the local sizes sum up to the global
399 : // size, otherwise there is big trouble!
400 0 : numeric_index_type dbg_sum=0;
401 :
402 0 : for (auto p : make_range(this->n_processors()))
403 0 : dbg_sum += local_sizes[p];
404 :
405 0 : libmesh_assert_equal_to (dbg_sum, n);
406 :
407 : # endif
408 :
409 : #else
410 :
411 : // No other options without MPI!
412 : libmesh_error_msg_if(n != n_local, "ERROR: MPI is required for n != n_local!");
413 :
414 : #endif
415 :
416 0 : _last_local_index = _first_local_index + n_local;
417 :
418 : // Set the initialized flag
419 0 : this->_is_initialized = true;
420 0 : this->_is_closed = true;
421 :
422 : // Zero the components unless directed otherwise
423 0 : if (!fast)
424 0 : this->zero();
425 0 : }
426 :
427 :
428 : template <typename T>
429 : inline
430 0 : void DistributedVector<T>::init (const numeric_index_type n,
431 : const numeric_index_type n_local,
432 : const std::vector<numeric_index_type> & /*ghost*/,
433 : const bool fast,
434 : const ParallelType ptype)
435 : {
436 : // TODO: we shouldn't ignore the ghost sparsity pattern
437 0 : this->init(n, n_local, fast, ptype);
438 0 : }
439 :
440 :
441 :
442 : /* Default implementation for solver packages for which ghosted
443 : vectors are not yet implemented. */
444 : template <class T>
445 0 : void DistributedVector<T>::init (const NumericVector<T> & other,
446 : const bool fast)
447 : {
448 0 : this->init(other.size(),other.local_size(),fast,other.type());
449 0 : }
450 :
451 :
452 :
453 : template <typename T>
454 : inline
455 0 : void DistributedVector<T>::init (const numeric_index_type n,
456 : const bool fast,
457 : const ParallelType ptype)
458 : {
459 0 : this->init(n,n,fast,ptype);
460 0 : }
461 :
462 :
463 :
464 : template <typename T>
465 : inline
466 0 : void DistributedVector<T>::clear ()
467 : {
468 0 : _values.clear();
469 0 : _remote_values.clear();
470 :
471 0 : _unclosed_state = DO_NOTHING;
472 :
473 0 : _global_size =
474 0 : _local_size =
475 0 : _first_local_index =
476 0 : _last_local_index = 0;
477 :
478 0 : this->_is_closed = this->_is_initialized = false;
479 0 : }
480 :
481 :
482 :
483 : template <typename T>
484 : inline
485 0 : void DistributedVector<T>::zero ()
486 : {
487 0 : libmesh_assert (this->initialized());
488 0 : libmesh_assert_equal_to (_values.size(), _local_size);
489 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
490 :
491 : // We should be closed
492 0 : libmesh_assert (this->closed());
493 0 : libmesh_assert_equal_to (_unclosed_state, DO_NOTHING);
494 0 : libmesh_assert (_remote_values.empty());
495 :
496 0 : std::fill (_values.begin(),
497 : _values.end(),
498 0 : 0.);
499 0 : }
500 :
501 :
502 :
503 : template <typename T>
504 : inline
505 0 : std::unique_ptr<NumericVector<T>> DistributedVector<T>::zero_clone () const
506 : {
507 0 : std::unique_ptr<NumericVector<T>> cloned_vector =
508 : std::make_unique<DistributedVector<T>>(this->comm());
509 0 : cloned_vector->init(*this);
510 0 : return cloned_vector;
511 0 : }
512 :
513 :
514 :
515 : template <typename T>
516 : inline
517 0 : std::unique_ptr<NumericVector<T>> DistributedVector<T>::clone () const
518 : {
519 0 : std::unique_ptr<NumericVector<T>> cloned_vector =
520 : std::make_unique<DistributedVector<T>>(this->comm());
521 0 : cloned_vector->init(*this, true);
522 0 : *cloned_vector = *this;
523 0 : return cloned_vector;
524 0 : }
525 :
526 :
527 :
528 : template <typename T>
529 : inline
530 0 : numeric_index_type DistributedVector<T>::size () const
531 : {
532 0 : libmesh_assert (this->initialized());
533 0 : libmesh_assert_equal_to (_values.size(), _local_size);
534 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
535 :
536 0 : return _global_size;
537 : }
538 :
539 :
540 :
541 : template <typename T>
542 : inline
543 0 : numeric_index_type DistributedVector<T>::local_size () const
544 : {
545 0 : libmesh_assert (this->initialized());
546 0 : libmesh_assert_equal_to (_values.size(), _local_size);
547 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
548 :
549 414 : return _local_size;
550 : }
551 :
552 :
553 :
554 : template <typename T>
555 : inline
556 0 : numeric_index_type DistributedVector<T>::first_local_index () const
557 : {
558 0 : libmesh_assert (this->initialized());
559 0 : libmesh_assert_equal_to (_values.size(), _local_size);
560 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
561 :
562 0 : return _first_local_index;
563 : }
564 :
565 :
566 :
567 : template <typename T>
568 : inline
569 0 : numeric_index_type DistributedVector<T>::last_local_index () const
570 : {
571 0 : libmesh_assert (this->initialized());
572 0 : libmesh_assert_equal_to (_values.size(), _local_size);
573 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
574 :
575 0 : return _last_local_index;
576 : }
577 :
578 :
579 :
580 : template <typename T>
581 : inline
582 0 : T DistributedVector<T>::operator() (const numeric_index_type i) const
583 : {
584 0 : libmesh_assert (this->initialized());
585 0 : libmesh_assert_equal_to (_values.size(), _local_size);
586 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
587 0 : libmesh_assert ( ((i >= first_local_index()) &&
588 : (i < last_local_index())) );
589 :
590 0 : return _values[i - _first_local_index];
591 : }
592 :
593 :
594 :
595 : template <typename T>
596 : inline
597 0 : void DistributedVector<T>::set (const numeric_index_type i, const T value)
598 : {
599 0 : libmesh_assert (this->initialized());
600 0 : libmesh_assert_equal_to (_values.size(), _local_size);
601 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
602 0 : libmesh_assert_less (i, size());
603 0 : libmesh_assert (_unclosed_state == DO_NOTHING ||
604 : _unclosed_state == SET_VALUES);
605 :
606 0 : std::scoped_lock lock(this->_numeric_vector_mutex);
607 0 : if (i >= first_local_index() && i < last_local_index())
608 : {
609 0 : _values[i - _first_local_index] = value;
610 : }
611 : else
612 : {
613 0 : _remote_values.emplace_back(i, value);
614 0 : _unclosed_state = SET_VALUES;
615 : }
616 :
617 0 : this->_is_closed = false;
618 0 : }
619 :
620 :
621 :
622 :
623 : template <typename T>
624 : inline
625 0 : void DistributedVector<T>::add (const numeric_index_type i, const T value)
626 : {
627 0 : libmesh_assert (this->initialized());
628 0 : libmesh_assert_equal_to (_values.size(), _local_size);
629 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
630 0 : libmesh_assert_less (i, size());
631 :
632 0 : std::scoped_lock lock(this->_numeric_vector_mutex);
633 :
634 0 : if (i >= first_local_index() && i < last_local_index())
635 : {
636 0 : _values[i - _first_local_index] += value;
637 : }
638 : else
639 : {
640 0 : _remote_values.emplace_back(i, value);
641 0 : _unclosed_state = ADD_VALUES;
642 : }
643 :
644 0 : this->_is_closed = false;
645 0 : }
646 :
647 :
648 :
649 : template <typename T>
650 : inline
651 0 : Real DistributedVector<T>::min () const
652 : {
653 : // This function must be run on all processors at once
654 0 : parallel_object_only();
655 :
656 0 : libmesh_assert (this->initialized());
657 0 : libmesh_assert (this->closed());
658 0 : libmesh_assert_equal_to (_values.size(), _local_size);
659 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
660 :
661 0 : Real local_min = std::numeric_limits<Real>::max();
662 0 : for (auto v : _values)
663 0 : local_min = std::min(libmesh_real(v), local_min);
664 :
665 0 : this->comm().min(local_min);
666 :
667 0 : return local_min;
668 : }
669 :
670 :
671 :
672 : template <typename T>
673 : inline
674 0 : Real DistributedVector<T>::max() const
675 : {
676 : // This function must be run on all processors at once
677 0 : parallel_object_only();
678 :
679 0 : libmesh_assert (this->initialized());
680 0 : libmesh_assert (this->closed());
681 0 : libmesh_assert_equal_to (_values.size(), _local_size);
682 0 : libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
683 :
684 0 : Real local_max = -std::numeric_limits<Real>::max();
685 0 : for (auto v : _values)
686 0 : local_max = std::max(libmesh_real(v), local_max);
687 :
688 0 : this->comm().max(local_max);
689 :
690 0 : return local_max;
691 : }
692 :
693 :
694 : template <typename T>
695 : inline
696 0 : void DistributedVector<T>::swap (NumericVector<T> & other)
697 : {
698 0 : DistributedVector<T> & v = cast_ref<DistributedVector<T> &>(other);
699 :
700 0 : std::swap(_global_size, v._global_size);
701 0 : std::swap(_local_size, v._local_size);
702 0 : std::swap(_first_local_index, v._first_local_index);
703 0 : std::swap(_last_local_index, v._last_local_index);
704 :
705 0 : std::swap(this->_is_initialized, v._is_initialized);
706 0 : std::swap(this->_is_closed, v._is_closed);
707 0 : std::swap(this->_unclosed_state, v._unclosed_state);
708 :
709 : // This should be O(1) with any reasonable STL implementation
710 0 : std::swap(_values, v._values);
711 0 : std::swap(_remote_values, v._remote_values);
712 0 : }
713 :
714 : template <typename T>
715 : inline
716 0 : std::size_t DistributedVector<T>::max_allowed_id () const
717 : {
718 : // Uses a std:vector<T>, so our indexing matches that
719 0 : return std::numeric_limits<typename std::vector<T>::size_type>::max();
720 : }
721 :
722 : } // namespace libMesh
723 :
724 :
725 : #endif // LIBMESH_DISTRIBUTED_VECTOR_H
|
