distributed_vector.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
 
 
#include "libmesh/libmesh_common.h"
 
 
 
#ifndef LIBMESH_DISTRIBUTED_VECTOR_H
#define LIBMESH_DISTRIBUTED_VECTOR_H
 
// Local includes
#include "libmesh/int_range.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/parallel.h"
 
// C++ includes
#include <vector>
#include <algorithm>
#include <limits>
#include <mutex>
 
namespace libMesh
{
 
template <typename T>
class DistributedVector final : public NumericVector<T>
{
public:
 
  explicit
  DistributedVector (const Parallel::Communicator & comm,
                     const ParallelType = AUTOMATIC);
 
  explicit
  DistributedVector (const Parallel::Communicator & comm,
                     const numeric_index_type n,
                     const ParallelType ptype = AUTOMATIC);
 
  DistributedVector (const Parallel::Communicator & comm,
                     const numeric_index_type n,
                     const numeric_index_type n_local,
                     const ParallelType ptype = AUTOMATIC);
 
  DistributedVector (const Parallel::Communicator & comm,
                     const numeric_index_type N,
                     const numeric_index_type n_local,
                     const std::vector<numeric_index_type> & ghost,
                     const ParallelType ptype = AUTOMATIC);
 
  DistributedVector & operator= (const DistributedVector &);
 
  DistributedVector (DistributedVector &&) = default;
  DistributedVector (const DistributedVector &) = default;
  DistributedVector & operator= (DistributedVector &&) = default;
  virtual ~DistributedVector () = default;
 
  virtual void close () override;
 
  virtual void clear () override;
 
  virtual void zero () override;
 
  virtual std::unique_ptr<NumericVector<T>> zero_clone () const override;
 
  virtual std::unique_ptr<NumericVector<T>> clone () const override;
 
  virtual void init (const numeric_index_type N,
                     const numeric_index_type n_local,
                     const bool fast=false,
                     const ParallelType ptype=AUTOMATIC) override;
 
  virtual void init (const numeric_index_type N,
                     const bool fast=false,
                     const ParallelType ptype=AUTOMATIC) override;
 
  virtual void init (const numeric_index_type N,
                     const numeric_index_type n_local,
                     const std::vector<numeric_index_type> & ghost,
                     const bool fast = false,
                     const ParallelType = AUTOMATIC) override;
 
  virtual void init (const NumericVector<T> & other,
                     const bool fast = false) override;
 
  virtual NumericVector<T> & operator= (const T s) override;
 
  virtual NumericVector<T> & operator= (const NumericVector<T> & v) override;
 
  virtual NumericVector<T> & operator= (const std::vector<T> & v) override;
 
  virtual Real min () const override;
 
  virtual Real max () const override;
 
  virtual T sum() const override;
 
  virtual Real l1_norm () const override;
 
  virtual Real l2_norm () const override;
 
  virtual Real linfty_norm () const override;
 
  virtual numeric_index_type size () const override;
 
  virtual numeric_index_type local_size() const override;
 
  virtual numeric_index_type first_local_index() const override;
 
  virtual numeric_index_type last_local_index() const override;
 
  virtual T operator() (const numeric_index_type i) const override;
 
  virtual NumericVector<T> & operator += (const NumericVector<T> & v) override;
 
  virtual NumericVector<T> & operator -= (const NumericVector<T> & v) override;
 
  virtual NumericVector<T> & operator *= (const NumericVector<T> & v) override;
 
  virtual NumericVector<T> & operator /= (const NumericVector<T> & v) override;
 
  virtual void reciprocal() override;
 
  virtual void conjugate() override;
 
  virtual void set (const numeric_index_type i, const T value) override;
 
  virtual void add (const numeric_index_type i, const T value) override;
 
  virtual void add (const T s) override;
 
  virtual void add (const NumericVector<T> & V) override;
 
  virtual void add (const T a, const NumericVector<T> & v) override;
 
  using NumericVector<T>::add_vector;
 
  virtual void add_vector (const NumericVector<T> &,
                           const SparseMatrix<T> &) override
  { libmesh_not_implemented(); }
 
  virtual void add_vector_transpose (const NumericVector<T> &,
                                     const SparseMatrix<T> &) override
  { libmesh_not_implemented(); }
 
  virtual void scale (const T factor) override;
 
  virtual void abs() override;
 
  virtual T dot(const NumericVector<T> & V) const override;
 
  virtual void localize (std::vector<T> & v_local) const override;
 
  virtual void localize (NumericVector<T> & v_local) const override;
 
  virtual void localize (NumericVector<T> & v_local,
                         const std::vector<numeric_index_type> & send_list) const override;
 
  virtual void localize (std::vector<T> & v_local,
                         const std::vector<numeric_index_type> & indices) const override;
 
  virtual void localize (const numeric_index_type first_local_idx,
                         const numeric_index_type last_local_idx,
                         const std::vector<numeric_index_type> & send_list) override;
 
  virtual void localize_to_one (std::vector<T> & v_local,
                                const processor_id_type proc_id=0) const override;
 
  virtual void pointwise_mult (const NumericVector<T> & vec1,
                               const NumericVector<T> & vec2) override;
 
  virtual void pointwise_divide (const NumericVector<T> & vec1,
                                 const NumericVector<T> & vec2) override;
 
  virtual void swap (NumericVector<T> & v) override;
 
  virtual std::size_t max_allowed_id() const override;
 
private:
 
  std::vector<T> _values;
 
  std::vector<std::pair<numeric_index_type, T>> _remote_values;
 
  enum UnclosedState { DO_NOTHING = 0,
                       ADD_VALUES,
                       SET_VALUES,
                       INVALID_UNCLOSEDSTATE };
 
  UnclosedState _unclosed_state;
 
  numeric_index_type _global_size;
 
  numeric_index_type _local_size;
 
  numeric_index_type _first_local_index;
 
  numeric_index_type _last_local_index;
};
 
 
//--------------------------------------------------------------------------
// DistributedVector inline methods
template <typename T>
inline
DistributedVector<T>::DistributedVector (const Parallel::Communicator & comm_in,
                                         const ParallelType ptype) :
  NumericVector<T>(comm_in, ptype),
  _unclosed_state   (DO_NOTHING),
  _global_size      (0),
  _local_size       (0),
  _first_local_index(0),
  _last_local_index (0)
{
  this->_type = ptype;
}
 
 
 
template <typename T>
inline
DistributedVector<T>::DistributedVector (const Parallel::Communicator & comm_in,
                                         const numeric_index_type n,
                                         const ParallelType ptype)
  : NumericVector<T>(comm_in, ptype)
{
  this->init(n, n, false, ptype);
}
 
 
 
template <typename T>
inline
DistributedVector<T>::DistributedVector (const Parallel::Communicator & comm_in,
                                         const numeric_index_type n,
                                         const numeric_index_type n_local,
                                         const ParallelType ptype)
  : NumericVector<T>(comm_in, ptype)
{
  this->init(n, n_local, false, ptype);
}
 
 
 
template <typename T>
inline
DistributedVector<T>::DistributedVector (const Parallel::Communicator & comm_in,
                                         const numeric_index_type n,
                                         const numeric_index_type n_local,
                                         const std::vector<numeric_index_type> & ghost,
                                         const ParallelType ptype)
  : NumericVector<T>(comm_in, ptype)
{
  this->init(n, n_local, ghost, false, ptype);
}
 
 
 
template <typename T>
inline
void DistributedVector<T>::init (const numeric_index_type n,
                                 const numeric_index_type n_local,
                                 const bool fast,
                                 const ParallelType ptype)
{
  // This function must be run on all processors at once
  parallel_object_only();
 
  libmesh_assert_less_equal (n_local, n);
 
  if (ptype == AUTOMATIC)
    {
      if (n == n_local)
        this->_type = SERIAL;
      else
        this->_type = PARALLEL;
    }
  else if (ptype == GHOSTED &&
           n == n_local) // We can support GHOSTED with no ghosts...
    this->_type = SERIAL;
  else
    this->_type = ptype;
 
  libmesh_assert ((this->_type==SERIAL && n==n_local) ||
                  this->_type==PARALLEL);
 
  // Clear the data structures if already initialized
  if (this->initialized())
    this->clear();
 
  // Initialize data structures
  _unclosed_state = DO_NOTHING;
  _values.resize(n_local);
  _remote_values.clear();
  _local_size  = n_local;
  _global_size = n;
 
  _first_local_index = 0;
 
#ifdef LIBMESH_HAVE_MPI
 
  std::vector<numeric_index_type> local_sizes (this->n_processors(), 0);
 
  local_sizes[this->processor_id()] = n_local;
 
  this->comm().sum(local_sizes);
 
  // _first_local_index is the sum of _local_size
  // for all processor ids less than ours
  for (auto p : make_range(this->processor_id()))
    _first_local_index += local_sizes[p];
 
 
#  ifdef DEBUG
  // Make sure all the local sizes sum up to the global
  // size, otherwise there is big trouble!
  numeric_index_type dbg_sum=0;
 
  for (auto p : make_range(this->n_processors()))
    dbg_sum += local_sizes[p];
 
  libmesh_assert_equal_to (dbg_sum, n);
 
#  endif
 
#else
 
  // No other options without MPI!
  libmesh_error_msg_if(n != n_local, "ERROR:  MPI is required for n != n_local!");
 
#endif
 
  _last_local_index = _first_local_index + n_local;
 
  // Set the initialized flag
  this->_is_initialized = true;
  this->_is_closed = true;
 
  // Zero the components unless directed otherwise
  if (!fast)
    this->zero();
}
 
 
template <typename T>
inline
void DistributedVector<T>::init (const numeric_index_type n,
                                 const numeric_index_type n_local,
                                 const std::vector<numeric_index_type> & /*ghost*/,
                                 const bool fast,
                                 const ParallelType ptype)
{
  // TODO: we shouldn't ignore the ghost sparsity pattern
  this->init(n, n_local, fast, ptype);
}
 
 
 
/* Default implementation for solver packages for which ghosted
   vectors are not yet implemented.  */
template <class T>
void DistributedVector<T>::init (const NumericVector<T> & other,
                                 const bool fast)
{
  this->init(other.size(),other.local_size(),fast,other.type());
}
 
 
 
template <typename T>
inline
void DistributedVector<T>::init (const numeric_index_type n,
                                 const bool fast,
                                 const ParallelType ptype)
{
  this->init(n,n,fast,ptype);
}
 
 
 
template <typename T>
inline
void DistributedVector<T>::clear ()
{
  _values.clear();
  _remote_values.clear();
 
  _unclosed_state = DO_NOTHING;
 
  _global_size =
    _local_size =
    _first_local_index =
    _last_local_index = 0;
 
  this->_is_closed = this->_is_initialized = false;
}
 
 
 
template <typename T>
inline
void DistributedVector<T>::zero ()
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
 
  // We should be closed
  libmesh_assert (this->closed());
  libmesh_assert_equal_to (_unclosed_state, DO_NOTHING);
  libmesh_assert (_remote_values.empty());
 
  std::fill (_values.begin(),
             _values.end(),
             0.);
}
 
 
 
template <typename T>
inline
std::unique_ptr<NumericVector<T>> DistributedVector<T>::zero_clone () const
{
  std::unique_ptr<NumericVector<T>> cloned_vector =
    std::make_unique<DistributedVector<T>>(this->comm());
  cloned_vector->init(*this);
  return cloned_vector;
}
 
 
 
template <typename T>
inline
std::unique_ptr<NumericVector<T>> DistributedVector<T>::clone () const
{
  std::unique_ptr<NumericVector<T>> cloned_vector =
    std::make_unique<DistributedVector<T>>(this->comm());
  cloned_vector->init(*this, true);
  *cloned_vector = *this;
  return cloned_vector;
}
 
 
 
template <typename T>
inline
numeric_index_type DistributedVector<T>::size () const
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
 
  return _global_size;
}
 
 
 
template <typename T>
inline
numeric_index_type DistributedVector<T>::local_size () const
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
 
  return _local_size;
}
 
 
 
template <typename T>
inline
numeric_index_type DistributedVector<T>::first_local_index () const
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
 
  return _first_local_index;
}
 
 
 
template <typename T>
inline
numeric_index_type DistributedVector<T>::last_local_index () const
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
 
  return _last_local_index;
}
 
 
 
template <typename T>
inline
T DistributedVector<T>::operator() (const numeric_index_type i) const
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
  libmesh_assert ( ((i >= first_local_index()) &&
                    (i <  last_local_index())) );
 
  return _values[i - _first_local_index];
}
 
 
 
template <typename T>
inline
void DistributedVector<T>::set (const numeric_index_type i, const T value)
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
  libmesh_assert_less (i, size());
  libmesh_assert (_unclosed_state == DO_NOTHING ||
                  _unclosed_state == SET_VALUES);
 
  std::scoped_lock lock(this->_numeric_vector_mutex);
  if (i >= first_local_index() && i < last_local_index())
    {
      _values[i - _first_local_index] = value;
    }
  else
    {
      _remote_values.emplace_back(i, value);
      _unclosed_state = SET_VALUES;
    }
 
  this->_is_closed = false;
}
 
 
 
 
template <typename T>
inline
void DistributedVector<T>::add (const numeric_index_type i, const T value)
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
  libmesh_assert_less (i, size());
 
  std::scoped_lock lock(this->_numeric_vector_mutex);
 
  if (i >= first_local_index() && i < last_local_index())
    {
      _values[i - _first_local_index] += value;
    }
  else
    {
      _remote_values.emplace_back(i, value);
      _unclosed_state = ADD_VALUES;
    }
 
  this->_is_closed = false;
}
 
 
 
template <typename T>
inline
Real DistributedVector<T>::min () const
{
  // This function must be run on all processors at once
  parallel_object_only();
 
  libmesh_assert (this->initialized());
  libmesh_assert (this->closed());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
 
  Real local_min = std::numeric_limits<Real>::max();
  for (auto v : _values)
    local_min = std::min(libmesh_real(v), local_min);
 
  this->comm().min(local_min);
 
  return local_min;
}
 
 
 
template <typename T>
inline
Real DistributedVector<T>::max() const
{
  // This function must be run on all processors at once
  parallel_object_only();
 
  libmesh_assert (this->initialized());
  libmesh_assert (this->closed());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
 
  Real local_max = -std::numeric_limits<Real>::max();
  for (auto v : _values)
    local_max = std::max(libmesh_real(v), local_max);
 
  this->comm().max(local_max);
 
  return local_max;
}
 
 
template <typename T>
inline
void DistributedVector<T>::swap (NumericVector<T> & other)
{
  DistributedVector<T> & v = cast_ref<DistributedVector<T> &>(other);
 
  std::swap(_global_size, v._global_size);
  std::swap(_local_size, v._local_size);
  std::swap(_first_local_index, v._first_local_index);
  std::swap(_last_local_index, v._last_local_index);
 
  std::swap(this->_is_initialized, v._is_initialized);
  std::swap(this->_is_closed, v._is_closed);
  std::swap(this->_unclosed_state, v._unclosed_state);
 
  // This should be O(1) with any reasonable STL implementation
  std::swap(_values, v._values);
  std::swap(_remote_values, v._remote_values);
}
 
template <typename T>
inline
std::size_t DistributedVector<T>::max_allowed_id () const
{
  // Uses a std:vector<T>, so our indexing matches that
  return std::numeric_limits<typename std::vector<T>::size_type>::max();
}
 
} // namespace libMesh
 
 
#endif  // LIBMESH_DISTRIBUTED_VECTOR_H