distributed_vector.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2025 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



// Local includes
#include "libmesh/distributed_vector.h"

// libMesh includes
#include "libmesh/dense_vector.h"
#include "libmesh/dense_subvector.h"
#include "libmesh/int_range.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/tensor_tools.h"

// TIMPI includes
#include "timpi/parallel_implementation.h"
#include "timpi/parallel_sync.h"

// C++ includes
#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
#include <cmath> // for std::abs
#include <limits> // std::numeric_limits<T>::min()


namespace libMesh
{



//--------------------------------------------------------------------------
// DistributedVector methods
template <typename T>
T DistributedVector<T>::sum () const
{
  // This function must be run on all processors at once
  parallel_object_only();

  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  T local_sum = 0.;

  for (auto & val : _values)
    local_sum += val;

  this->comm().sum(local_sum);

  return local_sum;
}



template <typename T>
Real DistributedVector<T>::l1_norm () const
{
  // This function must be run on all processors at once
  parallel_object_only();

  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  Real local_l1 = 0.;

  for (auto & val : _values)
    local_l1 += std::abs(val);

  this->comm().sum(local_l1);

  return local_l1;
}



template <typename T>
Real DistributedVector<T>::l2_norm () const
{
  // This function must be run on all processors at once
  parallel_object_only();

  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  Real local_l2 = 0.;

  for (auto & val : _values)
    local_l2 += TensorTools::norm_sq(val);

  this->comm().sum(local_l2);

  return std::sqrt(local_l2);
}



template <typename T>
Real DistributedVector<T>::linfty_norm () const
{
  // This function must be run on all processors at once
  parallel_object_only();

  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  Real local_linfty = 0.;

  for (auto & val : _values)
    local_linfty  = std::max(local_linfty,
                             static_cast<Real>(std::abs(val))
                             ); // Note we static_cast so that both
                                // types are the same, as required
                                // by std::max

  this->comm().max(local_linfty);

  return local_linfty;
}



template <typename T>
NumericVector<T> & DistributedVector<T>::operator += (const NumericVector<T> & v)
{
  libmesh_assert (this->closed());
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  add(1., v);

  return *this;
}



template <typename T>
NumericVector<T> & DistributedVector<T>::operator -= (const NumericVector<T> & v)
{
  libmesh_assert (this->closed());
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  add(-1., v);

  return *this;
}



template <typename T>
NumericVector<T> & DistributedVector<T>::operator *= (const NumericVector<T> & v)
{
  libmesh_assert_equal_to(size(), v.size());

  const DistributedVector<T> & v_vec = cast_ref<const DistributedVector<T> &>(v);

  for (auto i : index_range(_values))
    _values[i] *= v_vec._values[i];

  return *this;
}



template <typename T>
NumericVector<T> & DistributedVector<T>::operator /= (const NumericVector<T> & v)
{
  libmesh_assert_equal_to(size(), v.size());

  const DistributedVector<T> & v_vec = cast_ref<const DistributedVector<T> &>(v);

  for (auto i : index_range(_values))
    _values[i] /= v_vec._values[i];

  return *this;
}




template <typename T>
void DistributedVector<T>::reciprocal()
{
  for (auto & val : _values)
    {
      // Don't divide by zero
      libmesh_assert_not_equal_to (val, T(0));

      val = 1. / val;
    }
}




template <typename T>
void DistributedVector<T>::conjugate()
{
  // Replace values by complex conjugate
  for (auto & val : _values)
    val = libmesh_conj(val);
}





template <typename T>
void DistributedVector<T>::add (const T v)
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  for (auto & val : _values)
    val += v;
}



template <typename T>
void DistributedVector<T>::add (const NumericVector<T> & v)
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  add (1., v);
}



template <typename T>
void DistributedVector<T>::add (const T a, const NumericVector<T> & v_in)
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  // Make sure the NumericVector passed in is really a DistributedVector
  const DistributedVector<T> * v = cast_ptr<const DistributedVector<T> *>(&v_in);
  libmesh_error_msg_if(!v, "Cannot add different types of NumericVectors.");

  for (auto i : index_range(_values))
    _values[i] += a * v->_values[i];
}



template <typename T>
void DistributedVector<T>::scale (const T factor)
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  for (auto & val : _values)
    val *= factor;
}

template <typename T>
void DistributedVector<T>::abs()
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  for (auto & val : _values)
    val = std::abs(val);
}





template <typename T>
T DistributedVector<T>::dot (const NumericVector<T> & V) const
{
  // This function must be run on all processors at once
  parallel_object_only();

  // Make sure the NumericVector passed in is really a DistributedVector
  const DistributedVector<T> * v = cast_ptr<const DistributedVector<T> *>(&V);

  // Make sure that the two vectors are distributed in the same way.
  libmesh_assert_equal_to ( this->first_local_index(), v->first_local_index() );
  libmesh_assert_equal_to ( this->last_local_index(), v->last_local_index()  );

  // The result of dotting together the local parts of the vector.
  T local_dot = 0;

  for (auto i : index_range(_values))
    local_dot += this->_values[i] * v->_values[i];

  // The local dot products are now summed via MPI
  this->comm().sum(local_dot);

  return local_dot;
}



template <typename T>
NumericVector<T> &
DistributedVector<T>::operator = (const T s)
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  for (auto & val : _values)
    val = s;

  return *this;
}



template <typename T>
NumericVector<T> &
DistributedVector<T>::operator = (const NumericVector<T> & v_in)
{
  // Make sure the NumericVector passed in is really a DistributedVector
  const DistributedVector<T> * v = cast_ptr<const DistributedVector<T> *>(&v_in);

  *this = *v;

  return *this;
}



template <typename T>
DistributedVector<T> &
DistributedVector<T>::operator = (const DistributedVector<T> & v)
{
  this->_is_initialized    = v._is_initialized;
  this->_is_closed         = v._is_closed;

  _global_size       = v._global_size;
  _local_size        = v._local_size;
  _first_local_index = v._first_local_index;
  _last_local_index  = v._last_local_index;

  if (v.local_size() == this->local_size())
    _values = v._values;

  else
    libmesh_error_msg("v.local_size() = " << v.local_size() << " must be equal to this->local_size() = " << this->local_size());

  return *this;
}



template <typename T>
NumericVector<T> &
DistributedVector<T>::operator = (const std::vector<T> & v)
{
  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  if (v.size() == local_size())
    _values = v;

  else if (v.size() == size())
    for (auto i : index_range(*this))
      _values[i-first_local_index()] = v[i];

  else
    libmesh_error_msg("Incompatible sizes in DistributedVector::operator=");

  return *this;
}



template <typename T>
inline
void DistributedVector<T>::close ()
{
  libmesh_assert (this->initialized());

  parallel_object_only();

  bool have_remote_values = !_remote_values.empty();
  this->comm().max(have_remote_values);

  if (have_remote_values)
    {
      bool someone_is_setting = (_unclosed_state == SET_VALUES);
      this->comm().max(someone_is_setting);

#ifndef NDEBUG
      // If *I* have remote values, I must know what to do with them.
      if (!_remote_values.empty())
        libmesh_assert(_unclosed_state == SET_VALUES ||
                       _unclosed_state == ADD_VALUES);

      // If *anyone* had remote values, we must be doing something,
      // and all ranks must agree on what we're doing.
      bool someone_is_adding = (_unclosed_state == ADD_VALUES);
      this->comm().max(someone_is_adding);

      libmesh_assert_not_equal_to(someone_is_setting, someone_is_adding);
#endif

      // We want to traverse in id order later, but we can't compare
      // values with default < in the case where Number==complex.
      std::sort(_remote_values.begin(), _remote_values.end(),
                [](auto a, auto b)
                { return a.first < b.first; });

      std::vector<numeric_index_type> last_local_indices;
      this->comm().allgather(_last_local_index, last_local_indices);

      std::map<processor_id_type, std::vector<std::pair<numeric_index_type,T>>>
        updates_to_send;

      processor_id_type p = 0;
      for (auto [i, v] : _remote_values)
        {
          while (i >= last_local_indices[p])
            {
              libmesh_assert_less(p, this->n_processors());
              ++p;
            }
          updates_to_send[p].emplace_back(i, v);
        }
      _remote_values.clear();

      auto action_functor =
        [this, someone_is_setting]
        (processor_id_type,
         const std::vector<std::pair<numeric_index_type,T>> & incoming_values)
        {
          for (auto [i, v] : incoming_values)
            {
              libmesh_assert_greater_equal(i, _first_local_index);
              libmesh_assert_less(i, _last_local_index);
              if (someone_is_setting)
                _values[i-_first_local_index] = v;
              else
                _values[i-_first_local_index] += v;
            }
        };

      Parallel::push_parallel_vector_data
        (this->comm(), updates_to_send, action_functor);

      _unclosed_state = DO_NOTHING;
    }

  this->_is_closed = true;
}



template <typename T>
void DistributedVector<T>::localize (NumericVector<T> & v_local_in) 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);

  DistributedVector<T> * v_local = cast_ptr<DistributedVector<T> *>(&v_local_in);

  v_local->_first_local_index = 0;

  v_local->_global_size =
    v_local->_local_size =
    v_local->_last_local_index = size();

  v_local->_is_initialized =
    v_local->_is_closed = true;

  // Call localize on the vector's values.  This will help
  // prevent code duplication
  localize (v_local->_values);

#ifndef LIBMESH_HAVE_MPI

  libmesh_assert_equal_to (local_size(), size());

#endif
}



template <typename T>
void DistributedVector<T>::localize (NumericVector<T> & v_local_in,
                                     const std::vector<numeric_index_type> &) 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);

  // TODO: We don't yet support the send list; this is inefficient:
  localize (v_local_in);
}



template <typename T>
void DistributedVector<T>::localize (std::vector<T> & v_local,
                                     const std::vector<numeric_index_type> & indices) const
{
  // Resize v_local so there is enough room to hold all the local values.
  v_local.resize(indices.size());

  // We need to know who has the values we want, so get everyone's _local_size
  std::vector<numeric_index_type> local_sizes;
  this->comm().allgather (_local_size, local_sizes);

  // Make a vector of partial sums of local sizes
  std::vector<numeric_index_type> local_size_sums(this->n_processors());
  local_size_sums[0] = local_sizes[0];
  for (auto i : IntRange<numeric_index_type>(1, local_sizes.size()))
    local_size_sums[i] = local_size_sums[i-1] + local_sizes[i];

  // We now fill in 'requested_ids' based on the indices.  Also keep
  // track of the local index (in the indices vector) for each of
  // these, since we need that when unpacking.
  std::map<processor_id_type, std::vector<numeric_index_type>>
    requested_ids, local_requested_ids;

  // We'll use this typedef a couple of times below.
  typedef typename std::vector<numeric_index_type>::iterator iter_t;

  // For each index in indices, determine which processor it is on.
  // This is an O(N*log(p)) algorithm that uses std::upper_bound().
  // Note: upper_bound() returns an iterator to the first entry which is
  // greater than the given value.
  for (auto i : index_range(indices))
    {
      iter_t ub = std::upper_bound(local_size_sums.begin(),
                                   local_size_sums.end(),
                                   indices[i]);

      processor_id_type on_proc = cast_int<processor_id_type>
        (std::distance(local_size_sums.begin(), ub));

      requested_ids[on_proc].push_back(indices[i]);
      local_requested_ids[on_proc].push_back(i);
    }

  auto gather_functor =
    [this]
    (processor_id_type, const std::vector<dof_id_type> & ids,
     std::vector<T> & values)
    {
      // The first send/receive we did was for indices, the second one will be
      // for corresponding floating point values, so create storage for that now...
      const std::size_t ids_size = ids.size();
      values.resize(ids_size);

      for (std::size_t i=0; i != ids_size; i++)
        {
          // The index of the requested value
          const numeric_index_type requested_index = ids[i];

          // Transform into local numbering, and get requested value.
          values[i] = _values[requested_index - _first_local_index];
        }
    };

  auto action_functor =
    [& v_local, & local_requested_ids]
    (processor_id_type pid,
     const std::vector<dof_id_type> &,
     const std::vector<T> & values)
    {
      // Now write the received values to the appropriate place(s) in v_local
      for (auto i : index_range(values))
        {
          libmesh_assert(local_requested_ids.count(pid));
          libmesh_assert_less(i, local_requested_ids[pid].size());

          // Get the index in v_local where this value needs to be inserted.
          const numeric_index_type local_requested_index =
            local_requested_ids[pid][i];

          // Actually set the value in v_local
          v_local[local_requested_index] = values[i];
        }
    };

  const T * ex = nullptr;
  Parallel::pull_parallel_vector_data
    (this->comm(), requested_ids, gather_functor, action_functor, ex);
}



template <typename T>
void DistributedVector<T>::localize (const numeric_index_type first_local_idx,
                                     const numeric_index_type last_local_idx,
                                     const std::vector<numeric_index_type> & send_list)
{
  // Only good for serial vectors
  libmesh_assert_equal_to (this->size(), this->local_size());
  libmesh_assert_greater (last_local_idx, first_local_idx);
  libmesh_assert_less_equal (send_list.size(), this->size());
  libmesh_assert_less (last_local_idx, this->size());

  const numeric_index_type my_size       = this->size();
  const numeric_index_type my_local_size = (last_local_idx - first_local_idx + 1);

  // Don't bother for serial cases
  if ((first_local_idx == 0) &&
      (my_local_size == my_size))
    return;


  // Build a parallel vector, initialize it with the local
  // parts of (*this)
  DistributedVector<T> parallel_vec(this->comm());

  parallel_vec.init (my_size, my_local_size, true, PARALLEL);

  // Copy part of *this into the parallel_vec
  for (numeric_index_type i=first_local_idx; i<=last_local_idx; i++)
    parallel_vec._values[i-first_local_idx] = _values[i];

  // localize like normal
  parallel_vec.localize (*this, send_list);
}



template <typename T>
void DistributedVector<T>::localize (std::vector<T> & v_local) const
{
  // This function must be run on all processors at once
  parallel_object_only();

  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  v_local = this->_values;

  this->comm().allgather (v_local);

#ifndef LIBMESH_HAVE_MPI
  libmesh_assert_equal_to (local_size(), size());
#endif
}



template <typename T>
void DistributedVector<T>::localize_to_one (std::vector<T> & v_local,
                                            const processor_id_type pid) const
{
  // This function must be run on all processors at once
  parallel_object_only();

  libmesh_assert (this->initialized());
  libmesh_assert_equal_to (_values.size(), _local_size);
  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);

  v_local = this->_values;

  this->comm().gather (pid, v_local);

#ifndef LIBMESH_HAVE_MPI
  libmesh_assert_equal_to (local_size(), size());
#endif
}



template <typename T>
void DistributedVector<T>::pointwise_mult (const NumericVector<T> &,
                                           const NumericVector<T> &)
//void DistributedVector<T>::pointwise_mult (const NumericVector<T> & vec1,
//   const NumericVector<T> & vec2)
{
  libmesh_not_implemented();
}

template <typename T>
void DistributedVector<T>::pointwise_divide (const NumericVector<T> &,
                                             const NumericVector<T> &)
{
  libmesh_not_implemented();
}

//--------------------------------------------------------------
// Explicit instantiations
template class LIBMESH_EXPORT DistributedVector<Number>;

} // namespace libMesh