diff --git a/src/TiledArray/conversions/foreach.h b/src/TiledArray/conversions/foreach.h
index d62772e553..1eb672fcc8 100644
--- a/src/TiledArray/conversions/foreach.h
+++ b/src/TiledArray/conversions/foreach.h
@@ -348,8 +348,7 @@ inline std::
       shape_type(world, tile_norms, arg.trange(), op_returns_void),
       arg.pmap());  // if Op returns void tile_norms contains scaled norms, so
                     // do not scale again
-  for (typename std::vector<datum_type>::const_iterator it = tiles.begin();
-       it != tiles.end(); ++it) {
+  for (auto it = tiles.begin(); it != tiles.end(); ++it) {
     const auto index = it->first;
     if (!result.is_zero(index)) result.set(it->first, it->second);
   }
diff --git a/src/TiledArray/dist_eval/contraction_eval.h b/src/TiledArray/dist_eval/contraction_eval.h
index 37e9e50516..af363caecb 100644
--- a/src/TiledArray/dist_eval/contraction_eval.h
+++ b/src/TiledArray/dist_eval/contraction_eval.h
@@ -895,7 +895,12 @@ class Summa
       // Initialize the reduction task
       ReducePairTask<op_type>* MADNESS_RESTRICT const reduce_task =
           reduce_tasks_ + t;
-      new (reduce_task) ReducePairTask<op_type>(TensorImpl_::world(), op_);
+      new (reduce_task) ReducePairTask<op_type>(TensorImpl_::world(), op_
+#ifdef TILEDARRAY_ENABLE_SUMMA_TRACE_INITIALIZE
+                                                ,
+                                                nullptr, t
+#endif  // TILEDARRAY_ENABLE_SUMMA_TRACE_INITIALIZE
+      );
     }
 
     return proc_grid_.local_size();
@@ -939,7 +944,12 @@ class Summa
           ss << index << " ";
 #endif  // TILEDARRAY_ENABLE_SUMMA_TRACE_INITIALIZE
 
-          new (reduce_task) ReducePairTask<op_type>(TensorImpl_::world(), op_);
+          new (reduce_task) ReducePairTask<op_type>(TensorImpl_::world(), op_
+#ifdef TILEDARRAY_ENABLE_SUMMA_TRACE_INITIALIZE
+                                                    ,
+                                                    nullptr, index
+#endif  // TILEDARRAY_ENABLE_SUMMA_TRACE_INITIALIZE
+          );
           ++tile_count;
         } else {
           // Construct an empty task to represent zero tiles.
diff --git a/src/TiledArray/expressions/blk_tsr_engine.h b/src/TiledArray/expressions/blk_tsr_engine.h
index 5157ccac09..c0dfe5d658 100644
--- a/src/TiledArray/expressions/blk_tsr_engine.h
+++ b/src/TiledArray/expressions/blk_tsr_engine.h
@@ -180,9 +180,9 @@ class BlkTsrEngineBase : public LeafEngine<Derived> {
   trange_type make_trange() const {
     const unsigned int rank = array_.trange().tiles_range().rank();
 
-    std::vector<TiledRange1> trange_data;
+    container::svector<TiledRange1> trange_data;
     trange_data.reserve(rank);
-    std::vector<std::size_t> trange1_data;
+    container::svector<std::size_t> trange1_data;
 
     // Get temporary data pointers
     const auto* MADNESS_RESTRICT const trange = array_.trange().data().data();
@@ -215,7 +215,7 @@ class BlkTsrEngineBase : public LeafEngine<Derived> {
   trange_type make_trange(const Permutation& perm) const {
     const unsigned int rank = array_.trange().tiles_range().rank();
 
-    std::vector<TiledRange1> trange_data;
+    container::svector<TiledRange1> trange_data;
     trange_data.reserve(rank);
     std::vector<std::size_t> trange1_data;
 
@@ -370,7 +370,7 @@ class BlkTsrEngine
     const unsigned int rank = trange_.tiles_range().rank();
 
     // Construct and allocate memory for the shift range
-    std::vector<long> range_shift;
+    container::svector<long> range_shift;
     range_shift.reserve(rank);
 
     // Get temporary data pointers
@@ -397,7 +397,7 @@ class BlkTsrEngine
     const unsigned int rank = trange_.tiles_range().rank();
 
     // Construct and allocate memory for the shift range
-    std::vector<long> range_shift(rank, 0l);
+    container::svector<long> range_shift(rank, 0l);
 
     // Get temporary data pointers
     const auto* MADNESS_RESTRICT const trange = array_.trange().data().data();
@@ -543,7 +543,7 @@ class ScalBlkTsrEngine
     const unsigned int rank = trange_.tiles_range().rank();
 
     // Construct and allocate memory for the shift range
-    std::vector<long> range_shift(rank, 0l);
+    container::svector<long> range_shift(rank, 0l);
 
     // Get temporary data pointers
     const auto* MADNESS_RESTRICT const trange = array_.trange().data().data();
diff --git a/src/TiledArray/expressions/contraction_helpers.h b/src/TiledArray/expressions/contraction_helpers.h
index b4096a5a6d..92696a365f 100644
--- a/src/TiledArray/expressions/contraction_helpers.h
+++ b/src/TiledArray/expressions/contraction_helpers.h
@@ -59,7 +59,7 @@ auto range_from_annotation(const IndexList_& target_idxs,
   using size_type = typename range_type::size_type;
   using extent_type = std::pair<size_type, size_type>;
 
-  std::vector<extent_type> ranges;  // Will be the ranges for each extent
+  container::svector<extent_type> ranges;  // Will be the ranges for each extent
   const auto& lrange = lhs.range();
   const auto& rrange = rhs.range();
 
@@ -82,7 +82,7 @@ auto trange_from_annotation(const IndexList_& target_idxs,
                             const IndexList_& lhs_idxs,
                             const IndexList_& rhs_idxs, LHSType&& lhs,
                             RHSType&& rhs) {
-  std::vector<TiledRange1> ranges;  // Will be the ranges for each extent
+  container::svector<TiledRange1> ranges;  // Will be the ranges for each extent
   const auto& lrange = lhs.trange();
   const auto& rrange = rhs.trange();
 
@@ -148,8 +148,8 @@ inline auto make_bound_annotation(const BipartiteIndexList& free_vars,
                                   const BipartiteIndexList& rhs_vars) {
   const auto bound_temp = bound_annotations(free_vars, lhs_vars, rhs_vars);
   BipartiteIndexList bound_vars(
-      std::vector<std::string>(bound_temp.begin(), bound_temp.end()),
-      std::vector<std::string>{});
+      container::svector<std::string>(bound_temp.begin(), bound_temp.end()),
+      container::svector<std::string>{});
   return bound_vars;
 }
 
diff --git a/src/TiledArray/expressions/expr.h b/src/TiledArray/expressions/expr.h
index 9cc004aaa6..d78aa26004 100644
--- a/src/TiledArray/expressions/expr.h
+++ b/src/TiledArray/expressions/expr.h
@@ -498,7 +498,8 @@ class Expr {
     // This step may involve communication when the tiles are moved from the
     // sub-block distribution to the array distribution.
     {
-      const std::vector<long> shift =
+      // N.B. must deep copy
+      const container::svector<long> shift =
           tsr.array().trange().make_tile_range(tsr.lower_bound()).lobound();
 
       std::shared_ptr<op_type> shift_op =
diff --git a/src/TiledArray/math/gemm_helper.h b/src/TiledArray/math/gemm_helper.h
index e3243e37c8..d33ec2029b 100644
--- a/src/TiledArray/math/gemm_helper.h
+++ b/src/TiledArray/math/gemm_helper.h
@@ -28,6 +28,7 @@
 
 #include <TiledArray/error.h>
 #include <TiledArray/math/blas.h>
+#include <TiledArray/util/vector.h>
 
 #include <vector>
 
@@ -170,7 +171,7 @@ class GemmHelper {
     const auto* MADNESS_RESTRICT const right_upper = right.upbound_data();
 
     // Create the start and finish indices
-    std::vector<std::size_t> lower, upper;
+    TiledArray::container::svector<std::size_t> lower, upper;
     lower.reserve(result_rank_);
     upper.reserve(result_rank_);
 
diff --git a/src/TiledArray/reduce_task.h b/src/TiledArray/reduce_task.h
index 03abd9c818..753ac5df58 100644
--- a/src/TiledArray/reduce_task.h
+++ b/src/TiledArray/reduce_task.h
@@ -668,6 +668,7 @@ class ReduceTask {
     Future<result_type> result_;  ///< The result of the reduction task
     madness::Spinlock lock_;      ///< Task lock
     madness::CallbackInterface* callback_;  ///< The completion callback
+    int task_id_;                           ///< Task id
 
    public:
     /// Implementation constructor
@@ -675,8 +676,10 @@ class ReduceTask {
     /// \param world The world that owns this task
     /// \param op The reduction operation
     /// \param callback The callback that will be invoked when this task
-    /// has completed
-    ReduceTaskImpl(World& world, opT op, madness::CallbackInterface* callback)
+    ///        has completed
+    /// \param task_id the task id (for debugging)
+    ReduceTaskImpl(World& world, opT op, madness::CallbackInterface* callback,
+                   int task_id = -1)
         : madness::TaskInterface(1, TaskAttributes::hipri()),
           world_(world),
           op_(op),
@@ -684,7 +687,8 @@ class ReduceTask {
           ready_object_(nullptr),
           result_(),
           lock_(),
-          callback_(callback) {}
+          callback_(callback),
+          task_id_(task_id) {}
 
     virtual ~ReduceTaskImpl() {}
 
@@ -745,11 +749,12 @@ class ReduceTask {
 
   /// \param world The world that owns this task
   /// \param op The reduction operation [ default = opT() ]
-  /// \param callback The callback that will be invoked when this task is
-  /// complete
+  /// \param callback The callback that will be invoked when
+  ///        this task is complete
+  /// \param task_id the task id (for debugging)
   ReduceTask(World& world, const opT& op = opT(),
-             madness::CallbackInterface* callback = nullptr)
-      : pimpl_(new ReduceTaskImpl(world, op, callback)), count_(0ul) {}
+             madness::CallbackInterface* callback = nullptr, int task_id = -1)
+      : pimpl_(new ReduceTaskImpl(world, op, callback, task_id)), count_(0ul) {}
 
   /// Move constructor
 
@@ -935,10 +940,12 @@ class ReducePairTask : public ReduceTask<ReducePairOpWrapper<opT> > {
   /// \param world The world that owns this task
   /// \param op The pair reduction operation [ default = opT() ]
   /// \param callback The callback that will be invoked when this task is
-  /// complete
+  ///        complete
+  /// \param task_id the task id (for debugging)
   ReducePairTask(World& world, const opT& op = opT(),
-                 madness::CallbackInterface* callback = nullptr)
-      : ReduceTask_(world, op_type(op), callback) {}
+                 madness::CallbackInterface* callback = nullptr,
+                 int task_id = -1)
+      : ReduceTask_(world, op_type(op), callback, task_id) {}
 
   /// Move constructor
 
diff --git a/src/TiledArray/special/diagonal_array.h b/src/TiledArray/special/diagonal_array.h
index 87e1c3f23f..825d66fd98 100644
--- a/src/TiledArray/special/diagonal_array.h
+++ b/src/TiledArray/special/diagonal_array.h
@@ -40,8 +40,9 @@ namespace detail {
 
 /// computes [min,max) describing the diagonal elements that rank-d Range
 /// contains; if the input Range contains no diagonal elements this return an
-/// empty Range \param[in] rng an input (rank-d) Range \return the range of
-/// diagonal elements, as a rank-1 Range
+/// empty Range
+/// \param[in] rng an input (rank-d) Range
+/// \return the range of diagonal elements, as a rank-1 Range
 inline Range diagonal_range(Range const &rng) {
   const auto rank = rng.rank();
   TA_ASSERT(rng.rank() > 0);
@@ -62,10 +63,13 @@ inline Range diagonal_range(Range const &rng) {
 }
 
 /// \brief computes shape data (i.e. Frobenius norms of the tiles) for a
-/// constant diagonal tensor \tparam T a numeric type \param trange a TiledRange
-/// of the result \param val value of the diagonal elements \return a
-/// Tensor<float> containing the Frobenius norms of the tiles of a DistArray
-/// with \p val on the diagonal and zeroes elsewhere
+/// constant diagonal tensor
+/// \tparam T a numeric type
+/// \param trange a TiledRange of the result
+/// \param val value of the diagonal elements
+/// \return a Tensor<float> containing the Frobenius norms of
+///         the tiles of a DistArray with \p val on the diagonal and
+///         zeroes elsewhere
 template <typename T>
 Tensor<float> diagonal_shape(TiledRange const &trange, T val) {
   Tensor<float> shape(trange.tiles_range(), 0.0);
@@ -79,7 +83,8 @@ Tensor<float> diagonal_shape(TiledRange const &trange, T val) {
   // shortest dimension
   while (diag_elem < diag_extent) {
     // Get the tile index corresponding to the current diagonal_elem
-    auto tile_idx = trange.element_to_tile(std::vector<int>(ndim, diag_elem));
+    auto tile_idx =
+        trange.element_to_tile(container::svector<long>(ndim, diag_elem));
     auto tile_range = trange.make_tile_range(tile_idx);
 
     // Compute the range of diagonal elements in the tile
@@ -98,13 +103,15 @@ Tensor<float> diagonal_shape(TiledRange const &trange, T val) {
 }
 
 /// \brief computes shape data (i.e. Frobenius norms of the tiles) for a
-/// non-constant diagonal tensor \tparam RandomAccessIterator an iterator over
-/// the range of diagonal elements \param[in] trange a TiledRange of the result
+/// non-constant diagonal tensor
+/// \tparam RandomAccessIterator an iterator over
+/// the range of diagonal elements
+/// \param[in] trange a TiledRange of the result
 /// \param[in] diagonals_begin the begin iterator of the range of the diagonals
 /// \param[in] diagonals_end the end iterator of the range of the diagonals; if
-/// not given, default initialized and thus will not be checked \return a
-/// Tensor<float> containing the Frobenius norms of the tiles of a DistArray
-/// with \p val on the diagonal and zeroes elsewhere
+/// not given, default initialized and thus will not be checked
+/// \return a Tensor<float> containing the Frobenius norms of the tiles of
+/// a DistArray with \p val on the diagonal and zeroes elsewhere
 template <typename RandomAccessIterator>
 std::enable_if_t<is_iterator<RandomAccessIterator>::value, Tensor<float>>
 diagonal_shape(TiledRange const &trange, RandomAccessIterator diagonals_begin,
@@ -221,9 +228,11 @@ write_diag_tiles_to_array_rng(Array &A, RandomAccessIterator diagonals_begin) {
 /// \brief Creates a constant diagonal DistArray
 
 /// Creates an array whose only nonzero values are the (hyper)diagonal elements
-/// (i.e. (n,n,n, ..., n) ), and they are all have the same value \tparam Policy
-/// the policy type of the resulting DistArray \tparam T a numeric type \param
-/// world The world for the array \param[in] trange The trange for the array
+/// (i.e. (n,n,n, ..., n) ), and they are all have the same value
+/// \tparam Policy the policy type of the resulting DistArray
+/// \tparam T a numeric type
+/// \param world The world for the array
+/// \param[in] trange The trange for the array
 /// \param[in] val The value of the diagonal elements
 /// \return a constant diagonal DistArray
 template <typename Array, typename T = double>
@@ -250,13 +259,15 @@ Array diagonal_array(World &world, TiledRange const &trange, T val = 1) {
 
 /// Creates an array whose only nonzero values are the (hyper)diagonal elements
 /// (i.e. (n,n,n, ..., n) ); the values of the diagonal elements are given by an
-/// input range \tparam Array a DistArray type \tparam RandomAccessIterator an
-/// iterator over the range of diagonal elements \param world The world for the
-/// array \param[in] trange The trange for the array \param[in] diagonals_begin
-/// the begin iterator of the range of the diagonals \param[in] diagonals_end
-/// the end iterator of the range of the diagonals; if not given, default
-/// initialized and thus will not be checked \return a constant diagonal
-/// DistArray
+/// input range
+/// \tparam Array a DistArray type
+/// \tparam RandomAccessIterator an iterator over the range of diagonal elements
+/// \param world The world for the array
+/// \param[in] trange The trange for the array
+/// \param[in] diagonals_begin the begin iterator of the range of the diagonals
+/// \param[in] diagonals_end the end iterator of the range of the diagonals;
+///            if not given, default initialized and thus will not be checked
+/// \return a constant diagonal DistArray
 template <typename Array, typename RandomAccessIterator>
 std::enable_if_t<detail::is_iterator<RandomAccessIterator>::value, Array>
 diagonal_array(World &world, TiledRange const &trange,
diff --git a/src/TiledArray/tile_interface/shift.h b/src/TiledArray/tile_interface/shift.h
index 2d19c8db46..6fad3b9428 100644
--- a/src/TiledArray/tile_interface/shift.h
+++ b/src/TiledArray/tile_interface/shift.h
@@ -91,12 +91,12 @@ using TiledArray::shift;
 using TiledArray::shift_to;
 
 template <typename T>
-using result_of_shift_t = typename std::decay<decltype(
-    shift(std::declval<T>(), std::declval<std::vector<long>>()))>::type;
+using result_of_shift_t = typename std::decay<decltype(shift(
+    std::declval<T>(), std::declval<container::svector<long>>()))>::type;
 
 template <typename T>
-using result_of_shift_to_t = typename std::decay<decltype(
-    shift_to(std::declval<T>(), std::declval<std::vector<long>>()))>::type;
+using result_of_shift_to_t = typename std::decay<decltype(shift_to(
+    std::declval<T>(), std::declval<container::svector<long>>()))>::type;
 
 template <typename Result, typename Arg, typename Enabler = void>
 class Shift {
diff --git a/src/TiledArray/tile_op/shift.h b/src/TiledArray/tile_op/shift.h
index 03f4c2cc90..316c8a323e 100644
--- a/src/TiledArray/tile_op/shift.h
+++ b/src/TiledArray/tile_op/shift.h
@@ -54,7 +54,7 @@ class Shift {
       Consumable && std::is_same<result_type, argument_type>::value;
 
  private:
-  std::vector<long> range_shift_;
+  container::svector<long> range_shift_;
 
   // Permuting tile evaluation function
   // These operations cannot consume the argument tile since this operation
@@ -105,7 +105,12 @@ class Shift {
   /// Default constructor
 
   /// Construct a no operation that does not permute the result tile
-  Shift(const std::vector<long>& range_shift) : range_shift_(range_shift) {}
+  template <typename IntegralRange,
+            typename = std::enable_if_t<
+                TiledArray::detail::is_integral_range_v<IntegralRange>>>
+  Shift(IntegralRange&& range_shift)
+      : range_shift_(std::forward<IntegralRange>(range_shift).begin(),
+                     std::forward<IntegralRange>(range_shift).end()) {}
 
   /// Shift and permute operator
 
@@ -164,8 +169,8 @@ class ScalShift {
       Consumable && std::is_same<result_type, argument_type>::value;
 
  private:
-  std::vector<long> range_shift_;  ///< Range shift array
-  scalar_type factor_;             ///< Scaling factor
+  container::svector<long> range_shift_;  ///< Range shift array
+  scalar_type factor_;                    ///< Scaling factor
 
  public:
   // Permuting tile evaluation function
@@ -214,8 +219,13 @@ class ScalShift {
   /// Default constructor
 
   /// Construct a no operation that does not permute the result tile
-  ScalShift(const std::vector<long>& range_shift, const scalar_type factor)
-      : range_shift_(range_shift), factor_(factor) {}
+  template <typename IntegralRange,
+            typename = std::enable_if_t<
+                TiledArray::detail::is_integral_range_v<IntegralRange>>>
+  ScalShift(IntegralRange&& range_shift, const scalar_type factor)
+      : range_shift_(std::forward<IntegralRange>(range_shift).begin(),
+                     std::forward<IntegralRange>(range_shift).end()),
+        factor_(factor) {}
 
   /// Shift and permute operator