scatter.hpp

Go to the documentation of this file.

// Copyright (c) 2020-2023 Simons Foundation
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0.txt
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
//
// Authors: Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "../basic_functions.hpp"
#include "../concepts.hpp"
#include "../exceptions.hpp"
#include "../traits.hpp"
 
#include <mpi/mpi.hpp>
 
#include <type_traits>
#include <utility>
#include <vector>

template <nda::Array A>
struct mpi::lazy<mpi::tag::scatter, A> {
  using value_type = typename std::decay_t<A>::value_type;

  using const_view_type = decltype(std::declval<const A>()());

  const_view_type rhs;

  mpi::communicator comm;

  const int root{0}; // NOLINT (const is fine here)

  const bool all{false}; // NOLINT (const is fine here)

  [[nodiscard]] auto shape() const {
    auto dims = rhs.shape();
    long dim0 = dims[0];
    mpi::broadcast(dim0, comm, root);
    dims[0] = mpi::chunk_length(dim0, comm.size(), comm.rank());
    return dims;
  }

  template <nda::Array T>
  void invoke(T &&target) const { // NOLINT (temporary views are allowed here)
    if (not target.is_contiguous()) NDA_RUNTIME_ERROR << "Error in MPI scatter for nda::Array: Target array needs to be contiguous";
    static_assert(std::decay_t<A>::layout_t::stride_order_encoded == std::decay_t<T>::layout_t::stride_order_encoded,
                  "Error in MPI scatter for nda::Array: Incompatible stride orders");
 
    // special case for non-mpi runs
    if (not mpi::has_env) {
      target = rhs;
      return;
    }
 
    // get target shape and resize or check the target array
    auto dims = shape();
    resize_or_check_if_view(target, dims);
 
    // compute send counts, receive counts and memory displacements
    auto dim0       = rhs.extent(0);
    auto stride0    = rhs.indexmap().strides()[0];
    auto sendcounts = std::vector<int>(comm.size());
    auto displs     = std::vector<int>(comm.size() + 1, 0);
    int recvcount   = mpi::chunk_length(dim0, comm.size(), comm.rank()) * stride0;
    for (int r = 0; r < comm.size(); ++r) {
      sendcounts[r] = mpi::chunk_length(dim0, comm.size(), r) * stride0;
      displs[r + 1] = sendcounts[r] + displs[r];
    }
 
    // scatter the data
    auto mpi_value_type = mpi::mpi_type<value_type>::get();
    MPI_Scatterv((void *)rhs.data(), &sendcounts[0], &displs[0], mpi_value_type, (void *)target.data(), recvcount, mpi_value_type, root, comm.get());
  }
};
 
namespace nda {

  template <typename A>
  ArrayInitializer<std::remove_reference_t<A>> auto mpi_scatter(A &&a, mpi::communicator comm = {}, int root = 0, bool all = false)
    requires(is_regular_or_view_v<A>)
  {
    if (not a.is_contiguous()) NDA_RUNTIME_ERROR << "Error in MPI scatter for nda::Array: Array needs to be contiguous";
    return mpi::lazy<mpi::tag::scatter, A>{std::forward<A>(a), comm, root, all};
  }
 
} // namespace nda