Merge branch 'master' into kmp5/feature/CP

Author: kmp5VT

INSTALL.md

@@ -43,7 +43,7 @@ Both methods are supported. However, for most users we _strongly_ recommend to b
 - [Range-V3](https://github.com/ericniebler/range-v3.git) -- a Ranges library that served as the basis for Ranges component of C++20 and later.
 - [BTAS](http://github.com/ValeevGroup/BTAS), tag 1cfcb12647c768ccd83b098c64cda723e1275e49 . If usable BTAS installation is not found, TiledArray will download and compile
   BTAS from source. *This is the recommended way to compile BTAS for all users*.
-- [MADNESS](https://github.com/m-a-d-n-e-s-s/madness), tag 93a9a5cec2a8fa87fba3afe8056607e6062a9058 .
+- [MADNESS](https://github.com/m-a-d-n-e-s-s/madness), tag ef97ad1f0080da04f9592f03185c1a331cd5e001 .
   Only the MADworld runtime and BLAS/LAPACK C API component of MADNESS is used by TiledArray.
   If usable MADNESS installation is not found, TiledArray will download and compile
   MADNESS from source. *This is the recommended way to compile MADNESS for all users*.

cmake/modules/FindOrFetchMADWorld.cmake

@@ -41,13 +41,15 @@ if (NOT TARGET MADworld)
 
   # look for C and MPI here to make troubleshooting easier and be able to override defaults for MADNESS
   enable_language(C)
-  find_package(MPI REQUIRED COMPONENTS C CXX)
+  find_package(MPI REQUIRED COMPONENTS C)
 
+  set(FETCHCONTENT_QUIET FALSE)
   include(FetchContent)
   FetchContent_Declare(
           MADNESS
-	      GIT_REPOSITORY https://github.com/m-a-d-n-e-s-s/madness.git
+	  GIT_REPOSITORY https://github.com/m-a-d-n-e-s-s/madness.git
           GIT_TAG ${TA_TRACKED_MADNESS_TAG}
+	  GIT_PROGRESS TRUE
   )
   FetchContent_MakeAvailable(MADNESS)
   FetchContent_GetProperties(MADNESS

external/versions.cmake

@@ -11,8 +11,8 @@ set(TA_INSTALL_EIGEN_PREVIOUS_VERSION 3.3.7)
 set(TA_INSTALL_EIGEN_URL_HASH SHA256=b4c198460eba6f28d34894e3a5710998818515104d6e74e5cc331ce31e46e626)
 set(TA_INSTALL_EIGEN_PREVIOUS_URL_HASH MD5=b9e98a200d2455f06db9c661c5610496)
 
-set(TA_TRACKED_MADNESS_TAG 8da56b1fc0b3d6eabe155923fb844f3430fc7d05)
-set(TA_TRACKED_MADNESS_PREVIOUS_TAG 93a9a5cec2a8fa87fba3afe8056607e6062a9058)
+set(TA_TRACKED_MADNESS_TAG bd84a52766ab497dedc2f15f2162fb0eb7ec4653)
+set(TA_TRACKED_MADNESS_PREVIOUS_TAG ef97ad1f0080da04f9592f03185c1a331cd5e001)
 set(TA_TRACKED_MADNESS_VERSION 0.10.1)
 set(TA_TRACKED_MADNESS_PREVIOUS_VERSION 0.10.1)
 

python/src/TiledArray/python/module.h

@@ -40,11 +40,8 @@ static World &initialize() {
   // this loads MPI before TA tries to do it
   int initialized = 0;
   MPI_Initialized(&initialized);
-  MPI_Comm ta_comm;
 
-  if (!initialized) {
-    ta_comm = MPI_COMM_WORLD;
-  } else {
+  if (initialized) {
     int thread_level;
     MPI_Query_thread(&thread_level);
     if (thread_level != MPI_THREAD_MULTIPLE)
@@ -58,10 +55,10 @@ static World &initialize() {
   char *_argv[0];
   char **argv = _argv;
   if (!madness::initialized()) {
-    madness::initialize(argc, argv, ta_comm);
+    madness::initialize(argc, argv);
     initialized_madness = true;
   }
-  TiledArray::World &world = TiledArray::initialize(argc, argv, ta_comm);
+  TiledArray::World &world = TiledArray::initialize(argc, argv, MPI_COMM_WORLD);
   TiledArray::set_default_world(world);
   if (world.rank() == 0) {
     std::cout << "initialized TA in a world with " << world.size() << " ranks"

src/TiledArray/conversions/make_array.h

@@ -28,6 +28,7 @@
 
 #include "TiledArray/array_impl.h"
 #include "TiledArray/external/madness.h"
+#include "TiledArray/pmap/replicated_pmap.h"
 #include "TiledArray/shape.h"
 #include "TiledArray/type_traits.h"
 
@@ -73,7 +74,7 @@ template <typename Array, typename Op,
           typename std::enable_if<is_dense<Array>::value>::type* = nullptr>
 inline Array make_array(
     World& world, const detail::trange_t<Array>& trange,
-    const std::shared_ptr<const detail::pmap_t<Array> >& pmap, Op&& op) {
+    const std::shared_ptr<const detail::pmap_t<Array>>& pmap, Op&& op) {
   typedef typename Array::value_type value_type;
   typedef typename value_type::range_type range_type;
 
@@ -150,10 +151,10 @@ template <typename Array, typename Op,
           typename std::enable_if<!is_dense<Array>::value>::type* = nullptr>
 inline Array make_array(
     World& world, const detail::trange_t<Array>& trange,
-    const std::shared_ptr<const detail::pmap_t<Array> >& pmap, Op&& op) {
+    const std::shared_ptr<const detail::pmap_t<Array>>& pmap, Op&& op) {
   typedef typename Array::value_type value_type;
   typedef typename Array::ordinal_type ordinal_type;
-  typedef std::pair<ordinal_type, Future<value_type> > datum_type;
+  typedef std::pair<ordinal_type, Future<value_type>> datum_type;
 
   // Create a vector to hold local tiles
   std::vector<datum_type> tiles;
@@ -241,6 +242,41 @@ inline Array make_array(World& world, const detail::trange_t<Array>& trange,
                            op);
 }
 
+/// a make_array variant that uses a sequence of tiles
+/// to construct a DistArray with default pmap
+template <typename Array, typename Tiles>
+Array make_array(World& world, const detail::trange_t<Array>& tiled_range,
+                 Tiles begin, Tiles end, bool replicated) {
+  Array array;
+  using Tuple = std::remove_reference_t<decltype(*begin)>;
+  using Index = std::tuple_element_t<0, Tuple>;
+  using shape_type = typename Array::shape_type;
+
+  std::shared_ptr<typename Array::pmap_interface> pmap;
+  if (replicated) {
+    size_t ntiles = tiled_range.tiles_range().volume();
+    pmap = std::make_shared<detail::ReplicatedPmap>(world, ntiles);
+  }
+
+  if constexpr (shape_type::is_dense()) {
+    array = Array(world, tiled_range, pmap);
+  } else {
+    std::vector<std::pair<Index, float>> tile_norms;
+    for (Tiles it = begin; it != end; ++it) {
+      auto [index, tile] = *it;
+      tile_norms.push_back({index, tile.norm()});
+    }
+    shape_type shape(world, tile_norms, tiled_range);
+    array = Array(world, tiled_range, shape, pmap);
+  }
+  for (Tiles it = begin; it != end; ++it) {
+    auto [index, tile] = *it;
+    if (array.is_zero(index)) continue;
+    array.set(index, tile);
+  }
+  return array;
+}
+
 }  // namespace TiledArray
 
 #endif  // TILEDARRAY_CONVERSIONS_MAKE_ARRAY_H__INCLUDED

src/TiledArray/dense_shape.h

@@ -89,6 +89,17 @@ class DenseShape {
     return false;
   }
 
+  /// Check that a tile is zero
+
+  /// \tparam Integer an integer type
+  /// \param i the index
+  /// \return false
+  template <typename Integer>
+  std::enable_if_t<std::is_integral_v<Integer>, bool> is_zero(
+      const std::initializer_list<Integer>& i) const {
+    return false;
+  }
+
   /// Check density
 
   /// \return true

src/TiledArray/dist_array.h

@@ -487,7 +487,9 @@ class DistArray : public madness::archive::ParallelSerializableObject {
   /// initialized using the `op` function/functor, which transforms
   /// each tile in `other` using `op`
   /// \param other The array to be copied
-  template <typename OtherTile, typename Op>
+  template <typename OtherTile, typename Op,
+            typename = std::enable_if_t<
+                !std::is_same_v<detail::remove_cvr_t<Op>, TiledRange>>>
   DistArray(const DistArray<OtherTile, Policy>& other, Op&& op) : pimpl_() {
     *this = foreach<Tile>(other, std::forward<Op>(op));
   }
@@ -1878,39 +1880,6 @@ auto norm2(const DistArray<Tile, Policy>& a) {
   return std::sqrt(squared_norm(a));
 }
 
-template <typename Array, typename Tiles>
-Array make_array(World& world, const detail::trange_t<Array>& tiled_range,
-                 Tiles begin, Tiles end, bool replicated) {
-  Array array;
-  using Tuple = std::remove_reference_t<decltype(*begin)>;
-  using Index = std::tuple_element_t<0, Tuple>;
-  using shape_type = typename Array::shape_type;
-
-  std::shared_ptr<typename Array::pmap_interface> pmap;
-  if (replicated) {
-    size_t ntiles = tiled_range.tiles_range().volume();
-    pmap = std::make_shared<detail::ReplicatedPmap>(world, ntiles);
-  }
-
-  if constexpr (shape_type::is_dense()) {
-    array = Array(world, tiled_range, pmap);
-  } else {
-    std::vector<std::pair<Index, float>> tile_norms;
-    for (Tiles it = begin; it != end; ++it) {
-      auto [index, tile] = *it;
-      tile_norms.push_back({index, tile.norm()});
-    }
-    shape_type shape(world, tile_norms, tiled_range);
-    array = Array(world, tiled_range, shape, pmap);
-  }
-  for (Tiles it = begin; it != end; ++it) {
-    auto [index, tile] = *it;
-    if (array.is_zero(index)) continue;
-    array.set(index, tile);
-  }
-  return array;
-}
-
 template <typename T, typename P>
 DistArray<T, P> replicated(const DistArray<T, P>& a) {
   auto& world = a.world();

src/TiledArray/einsum/tiledarray.h

@@ -1,6 +1,7 @@
 #ifndef TILEDARRAY_EINSUM_TILEDARRAY_H__INCLUDED
 #define TILEDARRAY_EINSUM_TILEDARRAY_H__INCLUDED
 
+#include "TiledArray/conversions/make_array.h"
 #include "TiledArray/dist_array.h"
 #include "TiledArray/einsum/index.h"
 #include "TiledArray/einsum/range.h"
@@ -410,6 +411,12 @@ template <DeNest DeNestFlag = DeNest::False, typename ArrayA_, typename ArrayB_,
 auto einsum(expressions::TsrExpr<ArrayA_> A, expressions::TsrExpr<ArrayB_> B,
             std::tuple<Einsum::Index<std::string>, Indices...> cs,
             World &world) {
+  // hotfix: process all preceding tasks before entering this code with many
+  // blocking calls
+  // TODO figure out why having free threads left after blocking MPI split
+  //  still not enough to ensure progress
+  world.gop.fence();
+
   using ArrayA = std::remove_cv_t<ArrayA_>;
   using ArrayB = std::remove_cv_t<ArrayB_>;
 
@@ -739,7 +746,8 @@ auto einsum(expressions::TsrExpr<ArrayA_> A, expressions::TsrExpr<ArrayB_> B,
     for (Index h : H.tiles) {
       auto &[A, B] = AB;
       auto own = A.own(h) || B.own(h);
-      auto comm = world.mpi.comm().Split(own, world.rank());
+      auto comm = madness::blocking_invoke(&SafeMPI::Intracomm::Split,
+                                           world.mpi.comm(), own, world.rank());
       worlds.push_back(std::make_unique<World>(comm));
       auto &owners = worlds.back();
       if (!own) continue;

src/TiledArray/sparse_shape.h

@@ -509,6 +509,17 @@ class SparseShape {
     return tile_norms_[i] < my_threshold_;
   }
 
+  /// Check that a tile is zero
+
+  /// \tparam Integer an integer type
+  /// \param i the index
+  /// \return true if tile at position \p i is zero
+  template <typename Integer>
+  std::enable_if_t<std::is_integral_v<Integer>, bool> is_zero(
+      const std::initializer_list<Integer>& i) const {
+    return this->is_zero<std::initializer_list<Integer>>(i);
+  }
+
   /// Check density
 
   /// \return true

src/TiledArray/tensor/tensor.h

@@ -431,8 +431,7 @@ class Tensor {
         auto volume = total_size();
         for (decltype(volume) i = 0; i < volume; ++i) {
           auto& el = *(data() + i);
-          if (!el.empty())
-            el = p(el, inner_perm);
+          if (!el.empty()) el = p(el, inner_perm);
         }
       }
     }
@@ -1633,10 +1632,19 @@ class Tensor {
     if (right.empty()) return *this;
     return binary(
         right,
-        [](const value_type& l, const value_t<Right>& r) -> decltype(auto) {
+        [](const value_type& l, const value_t<Right>& r) -> decltype(l + r) {
           if constexpr (detail::is_tensor_v<value_type>) {
-            if (l.empty() && r.empty())
-              return value_type{};
+            if (l.empty()) {
+              if (r.empty())
+                return {};
+              else
+                return r;
+            } else {
+              if (r.empty())
+                return l;
+              else
+                return l + r;
+            }
           }
           return l + r;
         });
@@ -1799,8 +1807,23 @@ class Tensor {
                 detail::tensors_have_equal_nested_rank_v<Tensor, Right>>>
   Tensor subt(const Right& right) const {
     return binary(
-        right, [](const value_type& l, const value_type& r) -> decltype(auto) {
-          return l - r;
+        right,
+        [](const value_type& l, const value_t<Right>& r) -> decltype(l - r) {
+          if constexpr (detail::is_tensor_v<value_type>) {
+            if (l.empty()) {
+              if (r.empty())
+                return {};
+              else
+                return -r;
+            } else {
+              if (r.empty())
+                return l;
+              else
+                return l - r;
+            }
+          } else {
+            return l - r;
+          }
         });
   }
 
@@ -1936,13 +1959,14 @@ class Tensor {
             typename std::enable_if<detail::is_nested_tensor_v<Right>>::type* =
                 nullptr>
   decltype(auto) mult(const Right& right) const {
-
-    auto mult_op =[](const value_type& l, const value_t<Right>& r) -> decltype(auto) {
+    auto mult_op = [](const value_type& l,
+                      const value_t<Right>& r) -> decltype(auto) {
       return l * r;
     };
 
     if (empty() || right.empty()) {
-      using res_t = decltype(std::declval<Tensor>().binary(std::declval<Right>(), mult_op));
+      using res_t = decltype(std::declval<Tensor>().binary(
+          std::declval<Right>(), mult_op));
       return res_t{};
     }