paper demo

INSTALL.md

@@ -42,7 +42,7 @@ Both methods are supported. However, for most users we _strongly_ recommend to b
   - Boost.Range: header-only, *only used for unit testing*
 - [BTAS](http://github.com/ValeevGroup/BTAS), tag 561fe1bff7f3374814111a15e28c7a141ab9b67a . 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 91fff76deba20c751d0646c54f2f1c1e07bd6156 .
+- [MADNESS](https://github.com/m-a-d-n-e-s-s/madness), tag 58b3e2c623d772f6e4a2e9cf5758073de32ecc50 .
   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*.

examples/demo/demo2.cpp

@@ -56,6 +56,10 @@ int main(int argc, char* argv[]) {
   assert((ρ.upbound() == Index{11, 9}));
   // extent of $\mathbb{Z}_{1,11} \otimes \mathbb{Z}_{-1,9}$
   assert((ρ.extent() == Index{10, 10}));
+  // 1st dimension of ρ is $\mathbb{Z}_{1,11}$
+  assert((ρ.dim(0) == Range1{1, 11}));
+  // 2nd dimension of ρ is $\mathbb{Z}_{-1,9}$
+  assert((ρ.dim(1) == Range1{-1, 9}));
   // the number of elements in $\mathbb{Z}_{1,11} \otimes \mathbb{Z}_{-1,9}$
   assert(ρ.volume() == 100);
   // row-major order
@@ -100,12 +104,51 @@ int main(int argc, char* argv[]) {
   auto tile_0_0 = τ.tile({0, 0});
   assert((tile_0_0 == Range{{1, 5}, {-1, 5}}));
 
-  // default instance of $\code{DistArray}$ is a dense array of $\code{double}$s
-  DistArray array0(world, τ);
-  array0.fill(1.0);  // fill $\code{array0}$ with 1's
-
-  // grab a tile NB this returns a ${\bf future}$ to a tile; see Section 3.2.
-  auto t00 = array0.find({0, 0});
+  // clang-format off
+
+  // 2-d array of $\code{double}$ 0s, indexed by ρ
+  Tensor<double> t0(ρ, 0.);
+  // same as $\code{t0}$ but filled with ordinals
+  TensorD t1(ρ, [&ρ](auto&& idx) {
+    return ρ.ordinal(idx);
+  });
+  // print out "0 1 .. 99 "
+  for (auto&& v : t1) cout << v << " ";
+  // same as $\code{t0}$, using existing buffer
+  shared_ptr<double[]> v(new double[ρ.volume()]);
+  TensorD t2(ρ, v); // t2 and v co-manage buffer lifetime
+  v[0] = 1.;
+  assert(t2(1, -1) == 1.);
+  // same as $\code{t0}$, using existing (unmanaged) buffer
+  auto t3 = make_map(v.get(), ρ);
+  v[0] = 2.;
+  assert(t3(1, -1) == 2.);
+  // Tensor has shallow-copy semantics
+  auto t4 = t0;
+  t0(1, -1) = 3.;
+  assert(t4(1, -1) == 3.);
+
+  // clang-format on
+
+  // default instance of $\code{DistArray}$ is
+  // a {\em dense} array of $\code{double}$s
+  // NB can use TArrayD instead of DistArray<>
+  DistArray<> a0(τ);
+  a0.fill(1.);  // fill $\code{da}$ with 1s
+  // every tile exists in a dense array
+  assert(!a0.is_zero({0, 0}));
+  // grab a ${\em future}$ to the {0,0} tile
+  auto t00 = a0.find({0, 0});
+
+  // shape of a {\em sparse} array over τ
+  // tiles with even ordinals ({0,0}, {0,2}, {1,1}) are zero
+  SparseShape s(TensorF(τ.tiles_range(), {0, 1, 0, 1, 0, 1}), τ);
+  // a sparse array of $\code{double}$s
+  // TSpArrayX $\equiv$ DistArray<TensorX, SparsePolicy>
+  TSpArrayD a1(τ, s);
+  // only some tiles are nonzero in sparse array
+  assert(a1.is_zero({0, 0}));
+  assert(!a1.is_zero({0, 1}));
 
 #endif  // defined(TILEDARRAY_CXX_COMPILER_SUPPORT_UNICODE_VARIABLES)
 

external/versions.cmake

@@ -19,8 +19,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 91fff76deba20c751d0646c54f2f1c1e07bd6156)
-set(TA_TRACKED_MADNESS_PREVIOUS_TAG 0b44ef319643cb9721fbe17d294987c146e6460e)
+set(TA_TRACKED_MADNESS_TAG 58b3e2c623d772f6e4a2e9cf5758073de32ecc50)
+set(TA_TRACKED_MADNESS_PREVIOUS_TAG dc3294160209cbd683bfb57cb2b933bd5f86e07e)
 set(TA_TRACKED_MADNESS_VERSION 0.10.1)
 set(TA_TRACKED_MADNESS_PREVIOUS_VERSION 0.10.1)
 

src/TiledArray/dist_array.h

@@ -158,7 +158,7 @@ class DistArray : public madness::archive::ParallelSerializableObject {
   using rebind_numeric_t = typename rebind_numeric<Numeric>::type;
 
  private:
-  pimpl_type pimpl_;  ///< managed ptr to Array implementation
+  pimpl_type pimpl_ = {};  ///< managed ptr to Array implementation
   bool defer_deleter_to_next_fence_ =
       false;  ///< if true, the impl object is scheduled to be destroyed in the
               ///< next fence
@@ -277,7 +277,7 @@ class DistArray : public madness::archive::ParallelSerializableObject {
   /// array is uninitialized, but these arrays may be assign via a tensor
   /// expression assignment or the copy construction.
 
-  DistArray() : pimpl_() {}
+  DistArray() = default;
 
   /// Copy constructor
 
@@ -298,6 +298,19 @@ class DistArray : public madness::archive::ParallelSerializableObject {
             const std::shared_ptr<const pmap_interface>& pmap = {})
       : pimpl_(init(world, trange, shape_type(1, trange), pmap)) {}
 
+  /// Dense array constructor
+
+  /// Constructs an array with the given meta data in default World.
+  /// This constructor only initializes the array meta data;
+  /// the array tiles are empty and must be assigned by the user.
+  /// \param trange The tiled range object that will be used to set the array
+  /// tiling.
+  /// \param pmap The tile index -> process map
+  explicit DistArray(const trange_type& trange,
+                     const std::shared_ptr<const pmap_interface>& pmap = {})
+      : pimpl_(init(get_default_world(), trange, shape_type(1, trange), pmap)) {
+  }
+
   /// Sparse array constructor
 
   /// Constructs an array with the given meta data. This constructor only
@@ -312,6 +325,19 @@ class DistArray : public madness::archive::ParallelSerializableObject {
                 std::shared_ptr<const pmap_interface>())
       : pimpl_(init(world, trange, shape, pmap)) {}
 
+  /// Sparse array constructor
+
+  /// Constructs an array with the given meta data in default World.
+  /// This constructor only initializes the array meta data; the array tiles
+  /// are empty and must be assigned by the user.
+  /// \param trange The tiled range object that will be used to set the array
+  /// tiling. \param shape The array shape that defines zero and non-zero tiles
+  /// \param pmap The tile index -> process map
+  DistArray(const trange_type& trange, const shape_type& shape,
+            const std::shared_ptr<const pmap_interface>& pmap =
+                std::shared_ptr<const pmap_interface>())
+      : pimpl_(init(get_default_world(), trange, shape, pmap)) {}
+
   /// \name Initializer list constructors
   /// \brief Creates a new tensor containing the elements in the provided
   ///         `std::initializer_list`.
@@ -374,6 +400,41 @@ class DistArray : public madness::archive::ParallelSerializableObject {
               std::initializer_list<std::initializer_list<T>>>>>>
           il)
       : DistArray(array_from_il<DistArray>(world, il)) {}
+
+  template <typename T>
+  explicit DistArray(std::initializer_list<T> il)  // N.B. clang does not like
+                                                   // detail::vector_il<T> here
+      : DistArray(array_from_il<DistArray>(get_default_world(), il)) {}
+
+  template <typename T>
+  explicit DistArray(std::initializer_list<std::initializer_list<T>> il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), il)) {}
+
+  template <typename T>
+  explicit DistArray(
+      std::initializer_list<std::initializer_list<std::initializer_list<T>>> il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), il)) {}
+
+  template <typename T>
+  explicit DistArray(std::initializer_list<std::initializer_list<
+                         std::initializer_list<std::initializer_list<T>>>>
+                         il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), il)) {}
+
+  template <typename T>
+  explicit DistArray(
+      std::initializer_list<std::initializer_list<std::initializer_list<
+          std::initializer_list<std::initializer_list<T>>>>>
+          il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), il)) {}
+
+  template <typename T>
+  explicit DistArray(
+      std::initializer_list<
+          std::initializer_list<std::initializer_list<std::initializer_list<
+              std::initializer_list<std::initializer_list<T>>>>>>
+          il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), il)) {}
   ///@}
 
   /// \name Tiling initializer list constructors
@@ -440,6 +501,44 @@ class DistArray : public madness::archive::ParallelSerializableObject {
               std::initializer_list<std::initializer_list<T>>>>>>
           il)
       : DistArray(array_from_il<DistArray>(world, trange, il)) {}
+
+  template <typename T>
+  DistArray(const trange_type& trange, std::initializer_list<T> il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), trange, il)) {}
+
+  template <typename T>
+  DistArray(const trange_type& trange,
+            std::initializer_list<std::initializer_list<T>> il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), trange, il)) {}
+
+  template <typename T>
+  DistArray(
+      const trange_type& trange,
+      std::initializer_list<std::initializer_list<std::initializer_list<T>>> il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), trange, il)) {}
+
+  template <typename T>
+  DistArray(const trange_type& trange,
+            std::initializer_list<std::initializer_list<
+                std::initializer_list<std::initializer_list<T>>>>
+                il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), trange, il)) {}
+
+  template <typename T>
+  DistArray(const trange_type& trange,
+            std::initializer_list<std::initializer_list<std::initializer_list<
+                std::initializer_list<std::initializer_list<T>>>>>
+                il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), trange, il)) {}
+
+  template <typename T>
+  DistArray(
+      const trange_type& trange,
+      std::initializer_list<
+          std::initializer_list<std::initializer_list<std::initializer_list<
+              std::initializer_list<std::initializer_list<T>>>>>>
+          il)
+      : DistArray(array_from_il<DistArray>(get_default_world(), trange, il)) {}
   /// @}
 
   /// converting copy constructor

src/TiledArray/external/btas.h

@@ -109,6 +109,34 @@ inline bool is_congruent(const btas::RangeNd<Order, Args...>& r1,
                     r2.extent_data());
 }
 
+/// Test if a BTAS range and a TA range are congruent
+
+/// This function tests that the rank and extent of
+/// \c r1 are equal to those of \c r2.
+/// \param r1 The first Range to compare
+/// \param r2 The second Range to compare
+template <blas::Layout Order, typename... Args>
+inline bool is_congruent(const btas::RangeNd<Order, Args...>& r1,
+                         const TiledArray::Range& r2) {
+  return (r1.rank() == r2.rank()) &&
+         std::equal(r1.extent_data(), r1.extent_data() + r1.rank(),
+                    r2.extent_data());
+}
+
+/// Test if a TA range and a BTAS range are congruent
+
+/// This function tests that the rank and extent of
+/// \c r1 are equal to those of \c r2.
+/// \param r1 The first Range to compare
+/// \param r2 The second Range to compare
+template <blas::Layout Order, typename... Args>
+inline bool is_congruent(const TiledArray::Range& r1,
+                         const btas::RangeNd<Order, Args...>& r2) {
+  return (r1.rank() == r2.rank()) &&
+         std::equal(r1.extent_data(), r1.extent_data() + r1.rank(),
+                    r2.extent_data());
+}
+
 template <typename T, typename Range, typename Storage>
 decltype(auto) make_ti(const btas::Tensor<T, Range, Storage>& arg) {
   return TiledArray::detail::TensorInterface<const T, TiledArray::Range,

src/TiledArray/external/madness.h

@@ -20,11 +20,6 @@
 #ifndef TILEDARRAY_EXTERNAL_MADNESS_H__INCLUDED
 #define TILEDARRAY_EXTERNAL_MADNESS_H__INCLUDED
 
-// This needs to be defined before world/worldreduce.h and world/worlddc.h
-#ifndef WORLD_INSTANTIATE_STATIC_TEMPLATES
-#define WORLD_INSTANTIATE_STATIC_TEMPLATES
-#endif  // WORLD_INSTANTIATE_STATIC_TEMPLATES
-
 #include <memory>
 
 #include <TiledArray/config.h>

src/TiledArray/external/umpire.h

@@ -91,14 +91,16 @@ class umpire_allocator_impl {
       : umpalloc_(umpalloc) {}
 
   template <class U>
-  umpire_allocator_impl(const umpire_allocator_impl<U>& rhs) noexcept
+  umpire_allocator_impl(
+      const umpire_allocator_impl<U, StaticLock>& rhs) noexcept
       : umpalloc_(rhs.umpalloc_) {}
 
   /// allocates memory using umpire dynamic pool
   pointer allocate(size_t n) {
     TA_ASSERT(umpalloc_);
 
-    size_t nbytes = n * sizeof(T);
+    // QuickPool::allocate_internal does not handle zero-size allocations
+    size_t nbytes = n == 0 ? 1 : n * sizeof(T);
     pointer result = nullptr;
     auto* allocation_strategy = umpalloc_->getAllocationStrategy();
 
@@ -117,7 +119,8 @@ class umpire_allocator_impl {
   void deallocate(pointer ptr, size_t n) {
     TA_ASSERT(umpalloc_);
 
-    const auto nbytes = n * sizeof(T);
+    // QuickPool::allocate_internal does not handle zero-size allocations
+    const auto nbytes = n == 0 ? 1 : n * sizeof(T);
     auto* allocation_strategy = umpalloc_->getAllocationStrategy();
 
     // N.B. with multiple threads would have to do this test in
@@ -137,15 +140,15 @@ class umpire_allocator_impl {
   umpire::Allocator* umpalloc_;
 };  // class umpire_allocator
 
-template <class T1, class T2>
-bool operator==(const umpire_allocator_impl<T1>& lhs,
-                const umpire_allocator_impl<T2>& rhs) noexcept {
+template <class T1, class T2, class StaticLock>
+bool operator==(const umpire_allocator_impl<T1, StaticLock>& lhs,
+                const umpire_allocator_impl<T2, StaticLock>& rhs) noexcept {
   return lhs.umpire_allocator() == rhs.umpire_allocator();
 }
 
-template <class T1, class T2>
-bool operator!=(const umpire_allocator_impl<T1>& lhs,
-                const umpire_allocator_impl<T2>& rhs) noexcept {
+template <class T1, class T2, class StaticLock>
+bool operator!=(const umpire_allocator_impl<T1, StaticLock>& lhs,
+                const umpire_allocator_impl<T2, StaticLock>& rhs) noexcept {
   return !(lhs == rhs);
 }
 

src/TiledArray/fwd.h

@@ -117,6 +117,11 @@ namespace symmetry {
 class Permutation;
 }
 
+// shapes
+class DenseShape;
+template <typename T = float>
+class SparseShape;
+
 // TiledArray Arrays
 template <typename, typename>
 class DistArray;

src/TiledArray/host/allocator.h

@@ -53,7 +53,9 @@ class host_allocator_impl
 
   template <class U>
   host_allocator_impl(const host_allocator_impl<U>& rhs) noexcept
-      : base_type(static_cast<const umpire_allocator_impl<U>&>(rhs)) {}
+      : base_type(static_cast<
+                  const umpire_allocator_impl<U, detail::MutexLock<hostEnv>>&>(
+            rhs)) {}
 
   template <typename T1, typename T2>
   friend bool operator==(const host_allocator_impl<T1>& lhs,