From 5ca181e7b55a7386910cd062f4076b80281496d9 Mon Sep 17 00:00:00 2001
From: Magdalena Schreter <schreter.magdalena@gmail.com>
Date: Thu, 24 Nov 2022 11:39:53 +0100
Subject: [PATCH] add utility function MPI::scatter

---
 doc/news/changes/minor/20221124Schreter   |  3 +
 include/deal.II/base/mpi.h                | 88 +++++++++++++++++++++++
 tests/base/mpi_scatter_01.cc              | 71 ++++++++++++++++++
 tests/base/mpi_scatter_01.mpirun=2.output |  2 +
 tests/base/mpi_scatter_01.mpirun=6.output |  2 +
 tests/base/mpi_scatter_01.output          |  2 +
 6 files changed, 168 insertions(+)
 create mode 100644 doc/news/changes/minor/20221124Schreter
 create mode 100644 tests/base/mpi_scatter_01.cc
 create mode 100644 tests/base/mpi_scatter_01.mpirun=2.output
 create mode 100644 tests/base/mpi_scatter_01.mpirun=6.output
 create mode 100644 tests/base/mpi_scatter_01.output

diff --git a/doc/news/changes/minor/20221124Schreter b/doc/news/changes/minor/20221124Schreter
new file mode 100644
index 0000000000..35c184a46c
--- /dev/null
+++ b/doc/news/changes/minor/20221124Schreter
@@ -0,0 +1,3 @@
+New: Add utility function Utilities::MPI::scatter.
+<br>
+(Magdalena Schreter, Peter Munch, 2022/11/24)
diff --git a/include/deal.II/base/mpi.h b/include/deal.II/base/mpi.h
index d460f1cab2..45432ac0d2 100644
--- a/include/deal.II/base/mpi.h
+++ b/include/deal.II/base/mpi.h
@@ -1281,6 +1281,25 @@ namespace Utilities
            const T &          object_to_send,
            const unsigned int root_process = 0);
 
+    /**
+     * A generalization of the classic MPI_Scatter function, that accepts
+     * arbitrary data types T, as long as boost::serialize accepts T as an
+     * argument.
+     *
+     * @param[in] comm MPI communicator.
+     * @param[in] object_to_send a vector of objects, with size equal to the
+     * number of processes.
+     * @param[in] root_process The process, which sends the objects to all
+     * processes. By default the process with rank 0 is the root process.
+     *
+     * @return Every process receives an object from the root_process.
+     */
+    template <typename T>
+    T
+    scatter(const MPI_Comm &      comm,
+            const std::vector<T> &object_to_send,
+            const unsigned int    root_process = 0);
+
     /**
      * This function sends an object @p object_to_send from the process @p root_process
      * to all other processes.
@@ -2093,6 +2112,75 @@ namespace Utilities
 
 
 
+    template <typename T>
+    T
+    scatter(const MPI_Comm &      comm,
+            const std::vector<T> &objects_to_send,
+            const unsigned int    root_process)
+    {
+#  ifndef DEAL_II_WITH_MPI
+      (void)comm;
+      (void)root_process;
+      return objects_to_send[0];
+#  else
+      const auto n_procs = dealii::Utilities::MPI::n_mpi_processes(comm);
+      const auto my_rank = dealii::Utilities::MPI::this_mpi_process(comm);
+
+      AssertIndexRange(root_process, n_procs);
+      AssertThrow(
+        my_rank != root_process || objects_to_send.size() == n_procs,
+        ExcMessage(
+          "The number of objects to be scattered must correspond to the number processes."));
+
+      std::vector<char> send_buffer;
+      std::vector<int>  send_counts;
+      std::vector<int>  send_displacements;
+
+      if (my_rank == root_process)
+        {
+          send_counts.resize(n_procs, 0);
+          send_displacements.resize(n_procs + 1, 0);
+
+          for (unsigned int i = 0; i < n_procs; ++i)
+            {
+              const auto packed_data = Utilities::pack(objects_to_send[i]);
+              send_buffer.insert(send_buffer.end(),
+                                 packed_data.begin(),
+                                 packed_data.end());
+              send_counts[i] = packed_data.size();
+            }
+
+          for (unsigned int i = 0; i < n_procs; ++i)
+            send_displacements[i + 1] = send_displacements[i] + send_counts[i];
+        }
+
+      int n_local_data;
+      MPI_Scatter(send_counts.data(),
+                  1,
+                  MPI_INT,
+                  &n_local_data,
+                  1,
+                  MPI_INT,
+                  root_process,
+                  comm);
+
+      std::vector<char> recv_buffer(n_local_data);
+
+      MPI_Scatterv(send_buffer.data(),
+                   send_counts.data(),
+                   send_displacements.data(),
+                   MPI_CHAR,
+                   recv_buffer.data(),
+                   n_local_data,
+                   MPI_CHAR,
+                   root_process,
+                   comm);
+
+      return Utilities::unpack<T>(recv_buffer);
+#  endif
+    }
+
+
     template <typename T>
     void
     broadcast(T *                buffer,
diff --git a/tests/base/mpi_scatter_01.cc b/tests/base/mpi_scatter_01.cc
new file mode 100644
index 0000000000..9c1dd7fe6f
--- /dev/null
+++ b/tests/base/mpi_scatter_01.cc
@@ -0,0 +1,71 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, 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.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test for Utilities::MPI::scatter
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/point.h>
+
+#include <vector>
+
+#include "../tests.h"
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  initlog();
+
+  const unsigned int root_process = 0;
+  const auto         n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  const auto         my_proc = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  // Creating the local array of points
+  std::vector<Point<3>> local_points(my_proc + 1);
+  for (unsigned int i = 0; i < my_proc + 1; ++i)
+    local_points[i] = Point<3>(my_proc, -my_proc, i);
+
+  // send to process 0
+  const auto gathered_points =
+    Utilities::MPI::gather(MPI_COMM_WORLD, local_points, root_process);
+
+  // scatter from process 0
+  const auto scattered_points =
+    Utilities::MPI::scatter(MPI_COMM_WORLD, gathered_points, root_process);
+
+  int test_passed = 1;
+
+  if (scattered_points.size() != my_proc + 1)
+    {
+      test_passed = 0;
+      deallog << "Error: Points on rank " << my_proc << " have wrong size. "
+              << std::endl;
+    }
+  for (unsigned int p = 0; p < scattered_points.size(); ++p)
+    if (scattered_points[p][0] != (double)my_proc ||
+        scattered_points[p][1] != (double)-my_proc ||
+        scattered_points[p][2] != (double)p)
+      {
+        test_passed = 0;
+        deallog << "Error with point " << p << " on rank " << my_proc
+                << std::endl;
+      }
+
+  if (Utilities::MPI::min(test_passed, MPI_COMM_WORLD))
+    deallog << "Test: ok" << std::endl;
+  else
+    deallog << "Test: FAILED" << std::endl;
+}
diff --git a/tests/base/mpi_scatter_01.mpirun=2.output b/tests/base/mpi_scatter_01.mpirun=2.output
new file mode 100644
index 0000000000..186b780267
--- /dev/null
+++ b/tests/base/mpi_scatter_01.mpirun=2.output
@@ -0,0 +1,2 @@
+
+DEAL::Test: ok
diff --git a/tests/base/mpi_scatter_01.mpirun=6.output b/tests/base/mpi_scatter_01.mpirun=6.output
new file mode 100644
index 0000000000..186b780267
--- /dev/null
+++ b/tests/base/mpi_scatter_01.mpirun=6.output
@@ -0,0 +1,2 @@
+
+DEAL::Test: ok
diff --git a/tests/base/mpi_scatter_01.output b/tests/base/mpi_scatter_01.output
new file mode 100644
index 0000000000..186b780267
--- /dev/null
+++ b/tests/base/mpi_scatter_01.output
@@ -0,0 +1,2 @@
+
+DEAL::Test: ok
-- 
2.39.5