diff --git a/src/utils/Messenger.hpp b/src/utils/Messenger.hpp
index cfe312a0884ae0718ee76c41855ce8b063fecffd..99de952b4d9e084a7981e0e0f43fe5feb4f8b4ab 100644
--- a/src/utils/Messenger.hpp
+++ b/src/utils/Messenger.hpp
@@ -93,6 +93,94 @@ class Messenger
   size_t m_rank{0};
   size_t m_size{1};
 
+  template <typename DataType>
+  void
+  _gather(const DataType& data, Array<DataType> gather, size_t rank) const
+  {
+    static_assert(std::is_arithmetic_v<DataType>);
+    Assert(gather.size() == m_size * (rank == m_rank));   // LCOV_EXCL_LINE
+
+#ifdef PUGS_HAS_MPI
+    MPI_Datatype mpi_datatype = Messenger::helper::mpiType<DataType>();
+
+    auto gather_address = (gather.size() > 0) ? &(gather[0]) : nullptr;
+
+    MPI_Gather(&data, 1, mpi_datatype, gather_address, 1, mpi_datatype, rank, MPI_COMM_WORLD);
+#else    // PUGS_HAS_MPI
+    gather[0] = data;
+#endif   // PUGS_HAS_MPI
+  }
+
+  template <template <typename... SendT> typename SendArrayType,
+            template <typename... RecvT>
+            typename RecvArrayType,
+            typename... SendT,
+            typename... RecvT>
+  void
+  _gather(const SendArrayType<SendT...>& data_array, RecvArrayType<RecvT...> gather_array, size_t rank) const
+  {
+    Assert(gather_array.size() == data_array.size() * m_size * (rank == m_rank));   // LCOV_EXCL_LINE
+
+    using SendDataType = typename SendArrayType<SendT...>::data_type;
+    using RecvDataType = typename RecvArrayType<RecvT...>::data_type;
+
+    static_assert(std::is_same_v<std::remove_const_t<SendDataType>, RecvDataType>);
+    static_assert(std::is_arithmetic_v<SendDataType>);
+
+#ifdef PUGS_HAS_MPI
+    MPI_Datatype mpi_datatype = Messenger::helper::mpiType<RecvDataType>();
+
+    auto data_address   = (data_array.size() > 0) ? &(data_array[0]) : nullptr;
+    auto gather_address = (gather_array.size() > 0) ? &(gather_array[0]) : nullptr;
+
+    MPI_Gather(data_address, data_array.size(), mpi_datatype, gather_address, data_array.size(), mpi_datatype, rank,
+               MPI_COMM_WORLD);
+#else    // PUGS_HAS_MPI
+    value_copy(data_array, gather_array);
+#endif   // PUGS_HAS_MPI
+  }
+
+  template <template <typename... SendT> typename SendArrayType,
+            template <typename... RecvT>
+            typename RecvArrayType,
+            typename... SendT,
+            typename... RecvT>
+  void
+  _gatherVariable(const SendArrayType<SendT...>& data_array,
+                  RecvArrayType<RecvT...> gather_array,
+                  Array<int> sizes_array,
+                  size_t rank) const
+  {
+    Assert(gather_array.size() - sum(sizes_array) * (rank == m_rank) == 0);   // LCOV_EXCL_LINE
+
+    using SendDataType = typename SendArrayType<SendT...>::data_type;
+    using RecvDataType = typename RecvArrayType<RecvT...>::data_type;
+
+    static_assert(std::is_same_v<std::remove_const_t<SendDataType>, RecvDataType>);
+    static_assert(std::is_arithmetic_v<SendDataType>);
+
+#ifdef PUGS_HAS_MPI
+    MPI_Datatype mpi_datatype = Messenger::helper::mpiType<RecvDataType>();
+    Array<int> start_positions{sizes_array.size()};
+    if (start_positions.size() > 0) {
+      start_positions[0] = 0;
+      for (size_t i = 1; i < start_positions.size(); ++i) {
+        start_positions[i] = start_positions[i - 1] + sizes_array[i - 1];
+      }
+    }
+
+    auto data_address      = (data_array.size() > 0) ? &(data_array[0]) : nullptr;
+    auto sizes_address     = (sizes_array.size() > 0) ? &(sizes_array[0]) : nullptr;
+    auto positions_address = (start_positions.size() > 0) ? &(start_positions[0]) : nullptr;
+    auto gather_address    = (gather_array.size() > 0) ? &(gather_array[0]) : nullptr;
+
+    MPI_Gatherv(data_address, data_array.size(), mpi_datatype, gather_address, sizes_address, positions_address,
+                mpi_datatype, rank, MPI_COMM_WORLD);
+#else    // PUGS_HAS_MPI
+    value_copy(data_array, gather_array);
+#endif   // PUGS_HAS_MPI
+  }
+
   template <typename DataType>
   void
   _allGather(const DataType& data, Array<DataType> gather) const
@@ -128,8 +216,51 @@ class Messenger
 #ifdef PUGS_HAS_MPI
     MPI_Datatype mpi_datatype = Messenger::helper::mpiType<RecvDataType>();
 
-    MPI_Allgather(&(data_array[0]), data_array.size(), mpi_datatype, &(gather_array[0]), data_array.size(),
-                  mpi_datatype, MPI_COMM_WORLD);
+    auto data_address   = (data_array.size() > 0) ? &(data_array[0]) : nullptr;
+    auto gather_address = (gather_array.size() > 0) ? &(gather_array[0]) : nullptr;
+
+    MPI_Allgather(data_address, data_array.size(), mpi_datatype, gather_address, data_array.size(), mpi_datatype,
+                  MPI_COMM_WORLD);
+#else    // PUGS_HAS_MPI
+    value_copy(data_array, gather_array);
+#endif   // PUGS_HAS_MPI
+  }
+
+  template <template <typename... SendT> typename SendArrayType,
+            template <typename... RecvT>
+            typename RecvArrayType,
+            typename... SendT,
+            typename... RecvT>
+  void
+  _allGatherVariable(const SendArrayType<SendT...>& data_array,
+                     RecvArrayType<RecvT...> gather_array,
+                     Array<int> sizes_array) const
+  {
+    Assert(gather_array.size() == static_cast<size_t>(sum(sizes_array)));   // LCOV_EXCL_LINE
+
+    using SendDataType = typename SendArrayType<SendT...>::data_type;
+    using RecvDataType = typename RecvArrayType<RecvT...>::data_type;
+
+    static_assert(std::is_same_v<std::remove_const_t<SendDataType>, RecvDataType>);
+    static_assert(std::is_arithmetic_v<SendDataType>);
+
+#ifdef PUGS_HAS_MPI
+    MPI_Datatype mpi_datatype = Messenger::helper::mpiType<RecvDataType>();
+    Array<int> start_positions{sizes_array.size()};
+    if (start_positions.size() > 0) {
+      start_positions[0] = 0;
+      for (size_t i = 1; i < start_positions.size(); ++i) {
+        start_positions[i] = start_positions[i - 1] + sizes_array[i - 1];
+      }
+    }
+
+    auto data_address      = (data_array.size() > 0) ? &(data_array[0]) : nullptr;
+    auto sizes_address     = (sizes_array.size() > 0) ? &(sizes_array[0]) : nullptr;
+    auto positions_address = (start_positions.size() > 0) ? &(start_positions[0]) : nullptr;
+    auto gather_address    = (gather_array.size() > 0) ? &(gather_array[0]) : nullptr;
+
+    MPI_Allgatherv(data_address, data_array.size(), mpi_datatype, gather_address, sizes_address, positions_address,
+                   mpi_datatype, MPI_COMM_WORLD);
 #else    // PUGS_HAS_MPI
     value_copy(data_array, gather_array);
 #endif   // PUGS_HAS_MPI
@@ -158,8 +289,10 @@ class Messenger
     static_assert(std::is_arithmetic_v<DataType>);
 
 #ifdef PUGS_HAS_MPI
+    auto array_address = (array.size() > 0) ? &(array[0]) : nullptr;
+
     MPI_Datatype mpi_datatype = Messenger::helper::mpiType<DataType>();
-    MPI_Bcast(&(array[0]), array.size(), mpi_datatype, root_rank, MPI_COMM_WORLD);
+    MPI_Bcast(array_address, array.size(), mpi_datatype, root_rank, MPI_COMM_WORLD);
 #endif   // PUGS_HAS_MPI
   }
 
@@ -185,7 +318,10 @@ class Messenger
 
     MPI_Datatype mpi_datatype = Messenger::helper::mpiType<SendDataType>();
 
-    MPI_Alltoall(&(sent_array[0]), count, mpi_datatype, &(recv_array[0]), count, mpi_datatype, MPI_COMM_WORLD);
+    auto sent_address = (sent_array.size() > 0) ? &(sent_array[0]) : nullptr;
+    auto recv_address = (recv_array.size() > 0) ? &(recv_array[0]) : nullptr;
+
+    MPI_Alltoall(sent_address, count, mpi_datatype, recv_address, count, mpi_datatype, MPI_COMM_WORLD);
 #else    // PUGS_HAS_MPI
     value_copy(sent_array, recv_array);
 #endif   // PUGS_HAS_MPI
@@ -400,6 +536,83 @@ class Messenger
 #endif   // PUGS_HAS_MPI
   }
 
+  template <typename DataType>
+  PUGS_INLINE Array<DataType>
+  gather(const DataType& data, size_t rank) const
+  {
+    static_assert(not std::is_const_v<DataType>);
+
+    Array<DataType> gather_array((m_rank == rank) ? m_size : 0);
+
+    if constexpr (std::is_arithmetic_v<DataType>) {
+      _gather(data, gather_array, rank);
+    } else if constexpr (is_trivially_castable<DataType>) {
+      using CastType = helper::split_cast_t<DataType>;
+
+      CastArray cast_value_array  = cast_value_to<const CastType>::from(data);
+      CastArray cast_gather_array = cast_array_to<CastType>::from(gather_array);
+
+      _gather(cast_value_array, cast_gather_array, rank);
+    } else {
+      static_assert(is_false_v<DataType>, "unexpected type of data");
+    }
+    return gather_array;
+  }
+
+  template <typename DataType>
+  PUGS_INLINE Array<std::remove_const_t<DataType>>
+  gather(const Array<DataType>& array, size_t rank) const
+  {
+    using MutableDataType = std::remove_const_t<DataType>;
+    Array<MutableDataType> gather_array((m_rank == rank) ? (m_size * array.size()) : 0);
+
+    if constexpr (std::is_arithmetic_v<DataType>) {
+      _gather(array, gather_array, rank);
+    } else if constexpr (is_trivially_castable<DataType>) {
+      using CastType        = helper::split_cast_t<DataType>;
+      using MutableCastType = helper::split_cast_t<MutableDataType>;
+
+      CastArray cast_array        = cast_array_to<CastType>::from(array);
+      CastArray cast_gather_array = cast_array_to<MutableCastType>::from(gather_array);
+
+      _gather(cast_array, cast_gather_array, rank);
+    } else {
+      static_assert(is_false_v<DataType>, "unexpected type of data");
+    }
+    return gather_array;
+  }
+
+  template <typename DataType>
+  PUGS_INLINE Array<std::remove_const_t<DataType>>
+  gatherVariable(const Array<DataType>& array, size_t rank) const
+  {
+    int send_size          = array.size();
+    Array<int> sizes_array = gather(send_size, rank);
+
+    using MutableDataType = std::remove_const_t<DataType>;
+    Array<MutableDataType> gather_array(sum(sizes_array));
+
+    if constexpr (std::is_arithmetic_v<DataType>) {
+      _gatherVariable(array, gather_array, sizes_array, rank);
+    } else if constexpr (is_trivially_castable<DataType>) {
+      using CastType        = helper::split_cast_t<DataType>;
+      using MutableCastType = helper::split_cast_t<MutableDataType>;
+
+      int size_ratio = sizeof(DataType) / sizeof(CastType);
+      for (size_t i = 0; i < sizes_array.size(); ++i) {
+        sizes_array[i] *= size_ratio;
+      }
+
+      CastArray cast_array        = cast_array_to<CastType>::from(array);
+      CastArray cast_gather_array = cast_array_to<MutableCastType>::from(gather_array);
+
+      _gatherVariable(cast_array, cast_gather_array, sizes_array, rank);
+    } else {
+      static_assert(is_false_v<DataType>, "unexpected type of data");
+    }
+    return gather_array;
+  }
+
   template <typename DataType>
   PUGS_INLINE Array<DataType>
   allGather(const DataType& data) const
@@ -446,6 +659,37 @@ class Messenger
     return gather_array;
   }
 
+  template <typename DataType>
+  PUGS_INLINE Array<std::remove_const_t<DataType>>
+  allGatherVariable(const Array<DataType>& array) const
+  {
+    int send_size          = array.size();
+    Array<int> sizes_array = allGather(send_size);
+
+    using MutableDataType = std::remove_const_t<DataType>;
+    Array<MutableDataType> gather_array(sum(sizes_array));
+
+    if constexpr (std::is_arithmetic_v<DataType>) {
+      _allGatherVariable(array, gather_array, sizes_array);
+    } else if constexpr (is_trivially_castable<DataType>) {
+      using CastType        = helper::split_cast_t<DataType>;
+      using MutableCastType = helper::split_cast_t<MutableDataType>;
+
+      int size_ratio = sizeof(DataType) / sizeof(CastType);
+      for (size_t i = 0; i < sizes_array.size(); ++i) {
+        sizes_array[i] *= size_ratio;
+      }
+
+      CastArray cast_array        = cast_array_to<CastType>::from(array);
+      CastArray cast_gather_array = cast_array_to<MutableCastType>::from(gather_array);
+
+      _allGatherVariable(cast_array, cast_gather_array, sizes_array);
+    } else {
+      static_assert(is_false_v<DataType>, "unexpected type of data");
+    }
+    return gather_array;
+  }
+
   template <typename SendDataType>
   PUGS_INLINE Array<std::remove_const_t<SendDataType>>
   allToAll(const Array<SendDataType>& sent_array) const
@@ -624,6 +868,27 @@ allReduceSum(const DataType& data)
   return messenger().allReduceSum(data);
 }
 
+template <typename DataType>
+PUGS_INLINE Array<DataType>
+gather(const DataType& data, size_t rank)
+{
+  return messenger().gather(data, rank);
+}
+
+template <typename DataType>
+PUGS_INLINE Array<std::remove_const_t<DataType>>
+gather(const Array<DataType>& array, size_t rank)
+{
+  return messenger().gather(array, rank);
+}
+
+template <typename DataType>
+PUGS_INLINE Array<std::remove_const_t<DataType>>
+gatherVariable(const Array<DataType>& array, size_t rank)
+{
+  return messenger().gatherVariable(array, rank);
+}
+
 template <typename DataType>
 PUGS_INLINE Array<DataType>
 allGather(const DataType& data)
@@ -638,6 +903,13 @@ allGather(const Array<DataType>& array)
   return messenger().allGather(array);
 }
 
+template <typename DataType>
+PUGS_INLINE Array<std::remove_const_t<DataType>>
+allGatherVariable(const Array<DataType>& array)
+{
+  return messenger().allGatherVariable(array);
+}
+
 template <typename DataType>
 PUGS_INLINE Array<std::remove_const_t<DataType>>
 allToAll(const Array<DataType>& array)
diff --git a/tests/test_Messenger.cpp b/tests/test_Messenger.cpp
index 22945c15e2e8ad244cdc25f27a40fdee7cba56f5..eebd5dc939c23c7875c3b92ff697cb1f00afca63 100644
--- a/tests/test_Messenger.cpp
+++ b/tests/test_Messenger.cpp
@@ -253,6 +253,372 @@ TEST_CASE("Messenger", "[mpi]")
     }
   }
 
+  SECTION("gather value to")
+  {
+    {
+      const size_t target_rank = 10 % parallel::size();
+      // simple type
+      size_t value{(3 + parallel::rank()) * 2};
+      Array<size_t> gather_array = parallel::gather(value, target_rank);
+      if (parallel::rank() == target_rank) {
+        REQUIRE(gather_array.size() == parallel::size());
+
+        for (size_t i = 0; i < gather_array.size(); ++i) {
+          REQUIRE((gather_array[i] == (3 + i) * 2));
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+
+    {
+      const size_t target_rank = 7 % parallel::size();
+      // trivial simple type
+      mpi_check::integer value{static_cast<int>((3 + parallel::rank()) * 2)};
+      Array<mpi_check::integer> gather_array = parallel::gather(value, target_rank);
+      if (parallel::rank() == target_rank) {
+        REQUIRE(gather_array.size() == parallel::size());
+
+        for (size_t i = 0; i < gather_array.size(); ++i) {
+          const int expected_value = (3 + i) * 2;
+          REQUIRE(gather_array[i] == expected_value);
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+
+    {
+      const size_t target_rank = 5 % parallel::size();
+      // compound trivial type
+      mpi_check::tri_int value{static_cast<int>((3 + parallel::rank()) * 2), static_cast<int>(2 + parallel::rank()),
+                               static_cast<int>(4 - parallel::rank())};
+      Array<mpi_check::tri_int> gather_array = parallel::gather(value, target_rank);
+
+      if (parallel::rank() == target_rank) {
+        REQUIRE(gather_array.size() == parallel::size());
+        for (size_t i = 0; i < gather_array.size(); ++i) {
+          mpi_check::tri_int expected_value{static_cast<int>((3 + i) * 2), static_cast<int>(2 + i),
+                                            static_cast<int>(4 - i)};
+          REQUIRE((gather_array[i] == expected_value));
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+  }
+
+  SECTION("gather array to")
+  {
+    {
+      const size_t target_rank = 17 % parallel::size();
+      // simple type
+      Array<int> array(3);
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (3 + parallel::rank()) * 2 + i;
+      }
+      Array<int> gather_array = parallel::gather(array, target_rank);
+      if (parallel::rank() == target_rank) {
+        REQUIRE(gather_array.size() == array.size() * parallel::size());
+
+        for (size_t i = 0; i < gather_array.size(); ++i) {
+          const int expected_value = (3 + i / array.size()) * 2 + (i % array.size());
+          REQUIRE((gather_array[i] == expected_value));
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+
+    {
+      const size_t target_rank = 13 % parallel::size();
+      // trivial simple type
+      Array<mpi_check::integer> array(3);
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (3 + parallel::rank()) * 2 + i;
+      }
+      Array<mpi_check::integer> gather_array = parallel::gather(array, target_rank);
+      if (parallel::rank() == target_rank) {
+        REQUIRE(gather_array.size() == array.size() * parallel::size());
+
+        for (size_t i = 0; i < gather_array.size(); ++i) {
+          const int expected_value = (3 + i / array.size()) * 2 + (i % array.size());
+          REQUIRE((gather_array[i] == expected_value));
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+
+    {
+      const size_t target_rank = 11 % parallel::size();
+      // compound trivial type
+      Array<mpi_check::tri_int> array(3);
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] =
+          mpi_check::tri_int{static_cast<int>((3 + parallel::rank()) * 2), static_cast<int>(2 + parallel::rank() + i),
+                             static_cast<int>(4 - parallel::rank() - i)};
+      }
+      Array<mpi_check::tri_int> gather_array = parallel::gather(array, target_rank);
+
+      if (parallel::rank() == target_rank) {
+        REQUIRE(gather_array.size() == array.size() * parallel::size());
+        for (size_t i = 0; i < gather_array.size(); ++i) {
+          mpi_check::tri_int expected_value{static_cast<int>((3 + i / array.size()) * 2),
+                                            static_cast<int>(2 + i / array.size() + (i % array.size())),
+                                            static_cast<int>(4 - i / array.size() - (i % array.size()))};
+          REQUIRE((gather_array[i] == expected_value));
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+  }
+
+  SECTION("gather variable array to")
+  {
+    {
+      const size_t target_rank = 17 % parallel::size();
+      // simple type
+      Array<size_t> array(3 + parallel::rank());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (parallel::rank() + i) * 2 + i;
+      }
+      Array<size_t> gather_array = parallel::gatherVariable(array, target_rank);
+
+      if (parallel::rank() == target_rank) {
+        const size_t gather_array_size = [&]() {
+          size_t size = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            size += (3 + i_rank);
+          }
+          return size;
+        }();
+
+        REQUIRE(gather_array.size() == gather_array_size);
+
+        {
+          size_t i = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            for (size_t j = 0; j < i_rank + 3; ++j) {
+              REQUIRE(gather_array[i++] == (i_rank + j) * 2 + j);
+            }
+          }
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+
+    {
+      const size_t target_rank = 0;
+      // trivial simple type
+      Array<mpi_check::integer> array(2 + 2 * parallel::rank());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (3 + parallel::rank()) * 2 + i;
+      }
+      Array<mpi_check::integer> gather_array = parallel::gatherVariable(array, target_rank);
+
+      if (parallel::rank() == target_rank) {
+        const size_t gather_array_size = [&]() {
+          size_t size = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            size += (2 + 2 * i_rank);
+          }
+          return size;
+        }();
+
+        REQUIRE(gather_array.size() == gather_array_size);
+
+        {
+          size_t i = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            for (size_t j = 0; j < 2 + 2 * i_rank; ++j) {
+              REQUIRE(gather_array[i++] == (3 + i_rank) * 2 + j);
+            }
+          }
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+
+    {
+      const size_t target_rank = 13 % parallel::size();
+      // compound trivial type
+      Array<mpi_check::tri_int> array(3 * parallel::rank() + 1);
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] =
+          mpi_check::tri_int{static_cast<int>((3 + parallel::rank()) * 2), static_cast<int>(2 + parallel::rank() + i),
+                             static_cast<int>(4 - parallel::rank() - i)};
+      }
+      Array<mpi_check::tri_int> gather_array = parallel::gatherVariable(array, target_rank);
+
+      if (parallel::rank() == target_rank) {
+        const size_t gather_array_size = [&]() {
+          size_t size = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            size += (3 * i_rank + 1);
+          }
+          return size;
+        }();
+
+        REQUIRE(gather_array.size() == gather_array_size);
+
+        {
+          size_t i = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            for (size_t j = 0; j < 3 * i_rank + 1; ++j) {
+              auto value = mpi_check::tri_int{static_cast<int>((3 + i_rank) * 2), static_cast<int>(2 + i_rank + j),
+                                              static_cast<int>(4 - i_rank - j)};
+
+              REQUIRE(gather_array[i++] == value);
+            }
+          }
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+  }
+
+  SECTION("gather variable array to (with some empty)")
+  {
+    {
+      const size_t target_rank = 17 % parallel::size();
+      // simple type
+      Array<size_t> array(3 * parallel::rank());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (parallel::rank() + i) * 2 + i;
+      }
+      Array<size_t> gather_array = parallel::gatherVariable(array, target_rank);
+
+      if (parallel::rank() == target_rank) {
+        const size_t gather_array_size = [&]() {
+          size_t size = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            size += (3 * i_rank);
+          }
+          return size;
+        }();
+
+        REQUIRE(gather_array.size() == gather_array_size);
+
+        {
+          size_t i = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            for (size_t j = 0; j < i_rank * 3; ++j) {
+              REQUIRE(gather_array[i++] == (i_rank + j) * 2 + j);
+            }
+          }
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+
+    {
+      const size_t target_rank = 0;
+      // trivial simple type
+      Array<mpi_check::integer> array((parallel::rank() < parallel::size() - 1) ? (2 + 2 * parallel::rank()) : 0);
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (3 + parallel::rank()) * 2 + i;
+      }
+      Array<mpi_check::integer> gather_array = parallel::gatherVariable(array, target_rank);
+
+      if (parallel::rank() == target_rank) {
+        const size_t gather_array_size = [&]() {
+          size_t size = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            size += (i_rank < parallel::size() - 1) ? (2 + 2 * i_rank) : 0;
+          }
+          return size;
+        }();
+
+        REQUIRE(gather_array.size() == gather_array_size);
+
+        {
+          size_t i = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            for (size_t j = 0; j < ((i_rank < parallel::size() - 1) ? (2 + 2 * i_rank) : 0); ++j) {
+              REQUIRE(gather_array[i++] == (3 + i_rank) * 2 + j);
+            }
+          }
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+
+    {
+      const size_t target_rank = 13 % parallel::size();
+      // compound trivial type
+      Array<mpi_check::tri_int> array(3 * parallel::rank());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] =
+          mpi_check::tri_int{static_cast<int>((3 + parallel::rank()) * 2), static_cast<int>(2 + parallel::rank() + i),
+                             static_cast<int>(4 - parallel::rank() - i)};
+      }
+      Array<mpi_check::tri_int> gather_array = parallel::gatherVariable(array, target_rank);
+
+      if (parallel::rank() == target_rank) {
+        const size_t gather_array_size = [&]() {
+          size_t size = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            size += (3 * i_rank);
+          }
+          return size;
+        }();
+
+        REQUIRE(gather_array.size() == gather_array_size);
+
+        {
+          size_t i = 0;
+          for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+            for (size_t j = 0; j < 3 * i_rank; ++j) {
+              auto value = mpi_check::tri_int{static_cast<int>((3 + i_rank) * 2), static_cast<int>(2 + i_rank + j),
+                                              static_cast<int>(4 - i_rank - j)};
+
+              REQUIRE(gather_array[i++] == value);
+            }
+          }
+        }
+      } else {
+        REQUIRE(gather_array.size() == 0);
+      }
+    }
+  }
+
+  SECTION("gather variable array to (with all empty arrays)")
+  {
+    {
+      const size_t target_rank = 7 % parallel::size();
+      // simple type [empty arrays]
+      Array<size_t> array;
+      Array<size_t> gather_array = parallel::gatherVariable(array, target_rank);
+
+      REQUIRE(gather_array.size() == 0);
+    }
+
+    {
+      const size_t target_rank = 0;
+      // trivial simple type
+      Array<mpi_check::integer> array;
+      Array<mpi_check::integer> gather_array = parallel::gatherVariable(array, target_rank);
+
+      REQUIRE(gather_array.size() == 0);
+    }
+
+    {
+      const size_t target_rank = 13 % parallel::size();
+      // compound trivial type
+      Array<mpi_check::tri_int> array;
+      Array<mpi_check::tri_int> gather_array = parallel::gatherVariable(array, target_rank);
+
+      REQUIRE(gather_array.size() == 0);
+    }
+  }
+
   SECTION("all gather value")
   {
     {
@@ -345,6 +711,242 @@ TEST_CASE("Messenger", "[mpi]")
     }
   }
 
+  SECTION("all gather empty array")
+  {
+    {
+      // simple type
+      Array<int> array;
+      Array<int> gather_array = parallel::allGather(array);
+      REQUIRE(gather_array.size() == 0);
+    }
+
+    {
+      // trivial simple type
+      Array<mpi_check::integer> array;
+      Array<mpi_check::integer> gather_array = parallel::allGather(array);
+      REQUIRE(gather_array.size() == 0);
+    }
+
+    {
+      // compound trivial type
+      Array<mpi_check::tri_int> array;
+      Array<mpi_check::tri_int> gather_array = parallel::allGather(array);
+
+      REQUIRE(gather_array.size() == 0);
+    }
+  }
+
+  SECTION("all gather variable array")
+  {
+    {
+      // simple type
+      Array<size_t> array(3 + parallel::rank());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (parallel::rank() + i) * 2 + i;
+      }
+      Array<size_t> gather_array = parallel::allGatherVariable(array);
+
+      const size_t gather_array_size = [&]() {
+        size_t size = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          size += (3 + i_rank);
+        }
+        return size;
+      }();
+
+      REQUIRE(gather_array.size() == gather_array_size);
+
+      {
+        size_t i = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          for (size_t j = 0; j < i_rank + 3; ++j) {
+            REQUIRE(gather_array[i++] == (i_rank + j) * 2 + j);
+          }
+        }
+      }
+    }
+
+    {
+      // trivial simple type
+      Array<mpi_check::integer> array(2 + 2 * parallel::rank());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (3 + parallel::rank()) * 2 + i;
+      }
+      Array<mpi_check::integer> gather_array = parallel::allGatherVariable(array);
+
+      const size_t gather_array_size = [&]() {
+        size_t size = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          size += (2 + 2 * i_rank);
+        }
+        return size;
+      }();
+
+      REQUIRE(gather_array.size() == gather_array_size);
+
+      {
+        size_t i = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          for (size_t j = 0; j < 2 + 2 * i_rank; ++j) {
+            REQUIRE(gather_array[i++] == (3 + i_rank) * 2 + j);
+          }
+        }
+      }
+    }
+
+    {
+      // compound trivial type
+      Array<mpi_check::tri_int> array(3 * parallel::rank() + 1);
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] =
+          mpi_check::tri_int{static_cast<int>((3 + parallel::rank()) * 2), static_cast<int>(2 + parallel::rank() + i),
+                             static_cast<int>(4 - parallel::rank() - i)};
+      }
+      Array<mpi_check::tri_int> gather_array = parallel::allGatherVariable(array);
+
+      const size_t gather_array_size = [&]() {
+        size_t size = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          size += (3 * i_rank + 1);
+        }
+        return size;
+      }();
+
+      REQUIRE(gather_array.size() == gather_array_size);
+
+      {
+        size_t i = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          for (size_t j = 0; j < 3 * i_rank + 1; ++j) {
+            auto value = mpi_check::tri_int{static_cast<int>((3 + i_rank) * 2), static_cast<int>(2 + i_rank + j),
+                                            static_cast<int>(4 - i_rank - j)};
+
+            REQUIRE(gather_array[i++] == value);
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("all gather variable array to (with some empty)")
+  {
+    {
+      // simple type
+      Array<size_t> array(3 * parallel::rank());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (parallel::rank() + i) * 2 + i;
+      }
+      Array<size_t> gather_array = parallel::allGatherVariable(array);
+
+      const size_t gather_array_size = [&]() {
+        size_t size = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          size += (3 * i_rank);
+        }
+        return size;
+      }();
+
+      REQUIRE(gather_array.size() == gather_array_size);
+
+      {
+        size_t i = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          for (size_t j = 0; j < i_rank * 3; ++j) {
+            REQUIRE(gather_array[i++] == (i_rank + j) * 2 + j);
+          }
+        }
+      }
+    }
+
+    {
+      // trivial simple type
+      Array<mpi_check::integer> array((parallel::rank() < parallel::size() - 1) ? (2 + 2 * parallel::rank()) : 0);
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] = (3 + parallel::rank()) * 2 + i;
+      }
+      Array<mpi_check::integer> gather_array = parallel::allGatherVariable(array);
+
+      const size_t gather_array_size = [&]() {
+        size_t size = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          size += (i_rank < parallel::size() - 1) ? (2 + 2 * i_rank) : 0;
+        }
+        return size;
+      }();
+
+      REQUIRE(gather_array.size() == gather_array_size);
+
+      {
+        size_t i = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          for (size_t j = 0; j < ((i_rank < parallel::size() - 1) ? (2 + 2 * i_rank) : 0); ++j) {
+            REQUIRE(gather_array[i++] == (3 + i_rank) * 2 + j);
+          }
+        }
+      }
+    }
+
+    {
+      // compound trivial type
+      Array<mpi_check::tri_int> array(3 * parallel::rank());
+      for (size_t i = 0; i < array.size(); ++i) {
+        array[i] =
+          mpi_check::tri_int{static_cast<int>((3 + parallel::rank()) * 2), static_cast<int>(2 + parallel::rank() + i),
+                             static_cast<int>(4 - parallel::rank() - i)};
+      }
+      Array<mpi_check::tri_int> gather_array = parallel::allGatherVariable(array);
+
+      const size_t gather_array_size = [&]() {
+        size_t size = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          size += (3 * i_rank);
+        }
+        return size;
+      }();
+
+      REQUIRE(gather_array.size() == gather_array_size);
+
+      {
+        size_t i = 0;
+        for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+          for (size_t j = 0; j < 3 * i_rank; ++j) {
+            auto value = mpi_check::tri_int{static_cast<int>((3 + i_rank) * 2), static_cast<int>(2 + i_rank + j),
+                                            static_cast<int>(4 - i_rank - j)};
+
+            REQUIRE(gather_array[i++] == value);
+          }
+        }
+      }
+    }
+  }
+
+  SECTION("gather variable array to (with all empty arrays)")
+  {
+    {
+      // simple type [empty arrays]
+      Array<size_t> array;
+      Array<size_t> gather_array = parallel::allGatherVariable(array);
+
+      REQUIRE(gather_array.size() == 0);
+    }
+
+    {
+      // trivial simple type
+      Array<mpi_check::integer> array;
+      Array<mpi_check::integer> gather_array = parallel::allGatherVariable(array);
+
+      REQUIRE(gather_array.size() == 0);
+    }
+
+    {
+      // compound trivial type
+      Array<mpi_check::tri_int> array;
+      Array<mpi_check::tri_int> gather_array = parallel::allGatherVariable(array);
+
+      REQUIRE(gather_array.size() == 0);
+    }
+  }
+
   SECTION("all array exchanges")
   {
     {   // simple type