libMesh
distributed_vector.h
Go to the documentation of this file.
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 
53 template <typename T>
54 class DistributedVector final : public NumericVector<T>
55 {
56 public:
57 
61  explicit
63  const ParallelType = AUTOMATIC);
64 
68  explicit
70  const numeric_index_type n,
71  const ParallelType ptype = AUTOMATIC);
72 
78  const numeric_index_type n,
79  const numeric_index_type n_local,
80  const ParallelType ptype = AUTOMATIC);
81 
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 
101 
106  DistributedVector (DistributedVector &&) = default;
107  DistributedVector (const DistributedVector &) = default;
109  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 
194 
195  virtual void add_vector (const NumericVector<T> &,
196  const SparseMatrix<T> &) override
197  { libmesh_not_implemented(); }
198 
199  virtual void add_vector_transpose (const NumericVector<T> &,
200  const SparseMatrix<T> &) override
201  { 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 
241  std::vector<T> _values;
242 
247  std::vector<std::pair<numeric_index_type, T>> _remote_values;
248 
257 
262 
267 
272 
277 
282 };
283 
284 
285 //--------------------------------------------------------------------------
286 // DistributedVector inline methods
287 template <typename T>
288 inline
290  const ParallelType ptype) :
291  NumericVector<T>(comm_in, ptype),
292  _unclosed_state (DO_NOTHING),
293  _global_size (0),
294  _local_size (0),
295  _first_local_index(0),
296  _last_local_index (0)
297 {
298  this->_type = ptype;
299 }
300 
301 
302 
303 template <typename T>
304 inline
306  const numeric_index_type n,
307  const ParallelType ptype)
308  : NumericVector<T>(comm_in, ptype)
309 {
310  this->init(n, n, false, ptype);
311 }
312 
313 
314 
315 template <typename T>
316 inline
318  const numeric_index_type n,
319  const numeric_index_type n_local,
320  const ParallelType ptype)
321  : NumericVector<T>(comm_in, ptype)
322 {
323  this->init(n, n_local, false, ptype);
324 }
325 
326 
327 
328 template <typename T>
329 inline
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  : NumericVector<T>(comm_in, ptype)
336 {
337  this->init(n, n_local, ghost, false, ptype);
338 }
339 
340 
341 
342 template <typename T>
343 inline
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  parallel_object_only();
351 
352  libmesh_assert_less_equal (n_local, n);
353 
354  if (ptype == AUTOMATIC)
355  {
356  if (n == n_local)
357  this->_type = SERIAL;
358  else
359  this->_type = PARALLEL;
360  }
361  else if (ptype == GHOSTED &&
362  n == n_local) // We can support GHOSTED with no ghosts...
363  this->_type = SERIAL;
364  else
365  this->_type = ptype;
366 
367  libmesh_assert ((this->_type==SERIAL && n==n_local) ||
368  this->_type==PARALLEL);
369 
370  // Clear the data structures if already initialized
371  if (this->initialized())
372  this->clear();
373 
374  // Initialize data structures
375  _unclosed_state = DO_NOTHING;
376  _values.resize(n_local);
377  _remote_values.clear();
378  _local_size = n_local;
379  _global_size = n;
380 
381  _first_local_index = 0;
382 
383 #ifdef LIBMESH_HAVE_MPI
384 
385  std::vector<numeric_index_type> local_sizes (this->n_processors(), 0);
386 
387  local_sizes[this->processor_id()] = n_local;
388 
389  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  for (auto p : make_range(this->processor_id()))
394  _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  numeric_index_type dbg_sum=0;
401 
402  for (auto p : make_range(this->n_processors()))
403  dbg_sum += local_sizes[p];
404 
405  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  _last_local_index = _first_local_index + n_local;
417 
418  // Set the initialized flag
419  this->_is_initialized = true;
420  this->_is_closed = true;
421 
422  // Zero the components unless directed otherwise
423  if (!fast)
424  this->zero();
425 }
426 
427 
428 template <typename T>
429 inline
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  this->init(n, n_local, fast, ptype);
438 }
439 
440 
441 
442 /* Default implementation for solver packages for which ghosted
443  vectors are not yet implemented. */
444 template <class T>
446  const bool fast)
447 {
448  this->init(other.size(),other.local_size(),fast,other.type());
449 }
450 
451 
452 
453 template <typename T>
454 inline
456  const bool fast,
457  const ParallelType ptype)
458 {
459  this->init(n,n,fast,ptype);
460 }
461 
462 
463 
464 template <typename T>
465 inline
467 {
468  _values.clear();
469  _remote_values.clear();
470 
471  _unclosed_state = DO_NOTHING;
472 
473  _global_size =
474  _local_size =
475  _first_local_index =
476  _last_local_index = 0;
477 
478  this->_is_closed = this->_is_initialized = false;
479 }
480 
481 
482 
483 template <typename T>
484 inline
486 {
487  libmesh_assert (this->initialized());
488  libmesh_assert_equal_to (_values.size(), _local_size);
489  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
490 
491  // We should be closed
492  libmesh_assert (this->closed());
493  libmesh_assert_equal_to (_unclosed_state, DO_NOTHING);
494  libmesh_assert (_remote_values.empty());
495 
496  std::fill (_values.begin(),
497  _values.end(),
498  0.);
499 }
500 
501 
502 
503 template <typename T>
504 inline
505 std::unique_ptr<NumericVector<T>> DistributedVector<T>::zero_clone () const
506 {
507  std::unique_ptr<NumericVector<T>> cloned_vector =
508  std::make_unique<DistributedVector<T>>(this->comm());
509  cloned_vector->init(*this);
510  return cloned_vector;
511 }
512 
513 
514 
515 template <typename T>
516 inline
517 std::unique_ptr<NumericVector<T>> DistributedVector<T>::clone () const
518 {
519  std::unique_ptr<NumericVector<T>> cloned_vector =
520  std::make_unique<DistributedVector<T>>(this->comm());
521  cloned_vector->init(*this, true);
522  *cloned_vector = *this;
523  return cloned_vector;
524 }
525 
526 
527 
528 template <typename T>
529 inline
531 {
532  libmesh_assert (this->initialized());
533  libmesh_assert_equal_to (_values.size(), _local_size);
534  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
535 
536  return _global_size;
537 }
538 
539 
540 
541 template <typename T>
542 inline
544 {
545  libmesh_assert (this->initialized());
546  libmesh_assert_equal_to (_values.size(), _local_size);
547  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
548 
549  return _local_size;
550 }
551 
552 
553 
554 template <typename T>
555 inline
557 {
558  libmesh_assert (this->initialized());
559  libmesh_assert_equal_to (_values.size(), _local_size);
560  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
561 
562  return _first_local_index;
563 }
564 
565 
566 
567 template <typename T>
568 inline
570 {
571  libmesh_assert (this->initialized());
572  libmesh_assert_equal_to (_values.size(), _local_size);
573  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
574 
575  return _last_local_index;
576 }
577 
578 
579 
580 template <typename T>
581 inline
583 {
584  libmesh_assert (this->initialized());
585  libmesh_assert_equal_to (_values.size(), _local_size);
586  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
587  libmesh_assert ( ((i >= first_local_index()) &&
588  (i < last_local_index())) );
589 
590  return _values[i - _first_local_index];
591 }
592 
593 
594 
595 template <typename T>
596 inline
598 {
599  libmesh_assert (this->initialized());
600  libmesh_assert_equal_to (_values.size(), _local_size);
601  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
602  libmesh_assert_less (i, size());
603  libmesh_assert (_unclosed_state == DO_NOTHING ||
604  _unclosed_state == SET_VALUES);
605 
606  std::scoped_lock lock(this->_numeric_vector_mutex);
607  if (i >= first_local_index() && i < last_local_index())
608  {
609  _values[i - _first_local_index] = value;
610  }
611  else
612  {
613  _remote_values.emplace_back(i, value);
614  _unclosed_state = SET_VALUES;
615  }
616 
617  this->_is_closed = false;
618 }
619 
620 
621 
622 
623 template <typename T>
624 inline
626 {
627  libmesh_assert (this->initialized());
628  libmesh_assert_equal_to (_values.size(), _local_size);
629  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
630  libmesh_assert_less (i, size());
631 
632  std::scoped_lock lock(this->_numeric_vector_mutex);
633 
634  if (i >= first_local_index() && i < last_local_index())
635  {
636  _values[i - _first_local_index] += value;
637  }
638  else
639  {
640  _remote_values.emplace_back(i, value);
641  _unclosed_state = ADD_VALUES;
642  }
643 
644  this->_is_closed = false;
645 }
646 
647 
648 
649 template <typename T>
650 inline
652 {
653  // This function must be run on all processors at once
654  parallel_object_only();
655 
656  libmesh_assert (this->initialized());
657  libmesh_assert (this->closed());
658  libmesh_assert_equal_to (_values.size(), _local_size);
659  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
660 
661  Real local_min = std::numeric_limits<Real>::max();
662  for (auto v : _values)
663  local_min = std::min(libmesh_real(v), local_min);
664 
665  this->comm().min(local_min);
666 
667  return local_min;
668 }
669 
670 
671 
672 template <typename T>
673 inline
675 {
676  // This function must be run on all processors at once
677  parallel_object_only();
678 
679  libmesh_assert (this->initialized());
680  libmesh_assert (this->closed());
681  libmesh_assert_equal_to (_values.size(), _local_size);
682  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
683 
684  Real local_max = -std::numeric_limits<Real>::max();
685  for (auto v : _values)
686  local_max = std::max(libmesh_real(v), local_max);
687 
688  this->comm().max(local_max);
689 
690  return local_max;
691 }
692 
693 
694 template <typename T>
695 inline
697 {
698  DistributedVector<T> & v = cast_ref<DistributedVector<T> &>(other);
699 
700  std::swap(_global_size, v._global_size);
701  std::swap(_local_size, v._local_size);
702  std::swap(_first_local_index, v._first_local_index);
703  std::swap(_last_local_index, v._last_local_index);
704 
705  std::swap(this->_is_initialized, v._is_initialized);
706  std::swap(this->_is_closed, v._is_closed);
707  std::swap(this->_unclosed_state, v._unclosed_state);
708 
709  // This should be O(1) with any reasonable STL implementation
710  std::swap(_values, v._values);
711  std::swap(_remote_values, v._remote_values);
712 }
713 
714 template <typename T>
715 inline
717 {
718  // Uses a std:vector<T>, so our indexing matches that
719  return std::numeric_limits<typename std::vector<T>::size_type>::max();
720 }
721 
722 } // namespace libMesh
723 
724 
725 #endif // LIBMESH_DISTRIBUTED_VECTOR_H
T libmesh_real(T a)
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:283
UnclosedState
Whether we are adding or setting remote values or neither - this determines behavior at the next clos...
virtual Real max() const override
virtual T dot(const NumericVector< T > &V) const override
numeric_index_type _last_local_index
The last component (+1) stored locally.
virtual void scale(const T factor) override
Scale each element of the vector by the given factor.
virtual void add_vector(const NumericVector< T > &, const SparseMatrix< T > &) override
Computes , i.e.
virtual numeric_index_type last_local_index() const override
virtual void close() override
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
virtual numeric_index_type size() const =0
virtual void set(const numeric_index_type i, const T value) override
Sets v(i) = value.
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) override
Subtracts v from *this, .
DistributedVector & operator=(const DistributedVector &)
Copy assignment operator.
std::vector< std::pair< numeric_index_type, T > > _remote_values
Entries to add or set on remote processors during the next close()
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: vector_fe_ex5.C:44
const Parallel::Communicator & comm() const
virtual NumericVector< T > & operator/=(const NumericVector< T > &v) override
Computes the component-wise division of this vector&#39;s entries by another&#39;s, .
virtual void add(const numeric_index_type i, const T value) override
Adds value to the vector entry specified by i.
bool _is_initialized
true once init() has been called.
virtual Real linfty_norm() const override
The libMesh namespace provides an interface to certain functionality in the library.
numeric_index_type _global_size
The global vector size.
virtual std::unique_ptr< NumericVector< T > > zero_clone() const override
const Number zero
.
Definition: libmesh.h:304
virtual numeric_index_type local_size() const override
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) override
Adds v to *this, .
uint8_t processor_id_type
Definition: id_types.h:104
Generic sparse matrix.
Definition: vector_fe_ex5.C:46
virtual void reciprocal() override
Computes the component-wise reciprocal, .
virtual void zero() override
Set all entries to zero.
virtual Real l1_norm() const override
This class provides a simple parallel, distributed vector datatype which is specific to libmesh...
dof_id_type numeric_index_type
Definition: id_types.h:99
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:257
virtual ~DistributedVector()=default
virtual numeric_index_type size() const override
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
libmesh_assert(ctx)
virtual std::size_t max_allowed_id() const override
virtual void pointwise_divide(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
ParallelType _type
Type of vector.
virtual std::unique_ptr< NumericVector< T > > clone() const override
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) override
Computes (summation not implied) i.e.
virtual void conjugate() override
Negates the imaginary component of each entry in the vector.
ParallelType type() const
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const override
Creates a local copy of the global vector in v_local only on processor proc_id.
numeric_index_type _local_size
The local vector size.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual numeric_index_type first_local_index() const override
virtual NumericVector< T > & operator*=(const NumericVector< T > &v) override
Computes the component-wise multiplication of this vector&#39;s entries by another&#39;s, ...
virtual void localize(std::vector< T > &v_local) const override
Creates a copy of the global vector in the local vector v_local.
virtual Real l2_norm() const override
virtual numeric_index_type local_size() const =0
virtual void abs() override
Sets for each entry in the vector.
static const bool value
Definition: xdr_io.C:54
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:276
virtual void init(const numeric_index_type N, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC) override
Change the dimension of the vector to n.
virtual void clear() override
Restores the NumericVector<T> to a pristine state.
std::vector< T > _values
Actual vector datatype to hold vector entries.
virtual Real min() const override
UnclosedState _unclosed_state
The current state of this vector, if it is unclosed.
virtual T sum() const override
virtual void swap(NumericVector< T > &v) override
Swaps the contents of this with v.
virtual void add_vector_transpose(const NumericVector< T > &, const SparseMatrix< T > &) override
Computes , i.e.
DistributedVector(const Parallel::Communicator &comm, const ParallelType=AUTOMATIC)
Dummy-Constructor.
numeric_index_type _first_local_index
The first component stored locally.
bool _is_closed
Flag which tracks whether the vector&#39;s values are consistent on all processors after insertion or add...
ParallelType
Defines an enum for parallel data structure types.
virtual T operator()(const numeric_index_type i) const override