range.h

Go to the documentation of this file.

/*
 *  This file is a part of TiledArray.
 *  Copyright (C) 2013  Virginia Tech
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef TILEDARRAY_RANGE_H__INCLUDED
#define TILEDARRAY_RANGE_H__INCLUDED
 
#include <TiledArray/permutation.h>
#include <TiledArray/range_iterator.h>
#include <TiledArray/size_array.h>
#include <TiledArray/util/vector.h>
 
namespace TiledArray {
 
 
class Range {
 public:
  typedef Range Range_;                
  typedef TA_1INDEX_TYPE index1_type;  
  typedef container::svector<index1_type>
      index_type;            
  typedef index_type index;  
  typedef detail::SizeArray<const index1_type>
      index_view_type;  
  typedef index_type
      extent_type;  
  typedef std::size_t ordinal_type;  
  typedef std::make_signed_t<ordinal_type> distance_type;  
  typedef ordinal_type size_type;  
  typedef detail::RangeIterator<index1_type, Range_>
      const_iterator;  
  friend class detail::RangeIterator<index1_type, Range_>;
 
  static_assert(detail::is_range_v<index_type>);  // index is a Range
 
 protected:
  container::svector<index1_type, 4 * TA_MAX_SOO_RANK_METADATA> datavec_;
  distance_type offset_ = 0l;  
  ordinal_type volume_ = 0ul;  
  unsigned int rank_ = 0u;  
 
  void init_datavec(unsigned int rank) { datavec_.resize(rank << 2); }
  const index1_type* data() const { return datavec_.data(); }
  index1_type* data() { return datavec_.data(); }
 
  index1_type* lobound_data_nc() { return data(); }
  index1_type* upbound_data_nc() { return data() + rank_; }
  index1_type* extent_data_nc() { return data() + (rank_ << 1); }
  index1_type* stride_data_nc() { return extent_data_nc() + rank_; }
 
 private:
 
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<detail::is_integral_range_v<Index1> &&
                                        detail::is_integral_range_v<Index2>>>
  void init_range_data(const Index1& lower_bound, const Index2& upper_bound) {
    // Construct temp pointers
    auto* MADNESS_RESTRICT const lower = data();
    auto* MADNESS_RESTRICT const upper = lower + rank_;
    auto* MADNESS_RESTRICT const extent = upper + rank_;
    auto* MADNESS_RESTRICT const stride = extent + rank_;
 
    // Set the volume seed
    volume_ = 1ul;
    offset_ = 0l;
 
    // initialize bounds and extents
    auto lower_it = std::begin(lower_bound);
    auto upper_it = std::begin(upper_bound);
    auto lower_end = std::end(lower_bound);
    auto upper_end = std::end(upper_bound);
    for (int d = 0; lower_it != lower_end && upper_it != upper_end;
         ++lower_it, ++upper_it, ++d) {
      // Compute data for element d of lower, upper, and extent
      const auto lower_bound_d = *lower_it;
      const auto upper_bound_d = *upper_it;
      lower[d] = lower_bound_d;
      upper[d] = upper_bound_d;
      // Check input dimensions
      TA_ASSERT(lower[d] <= upper[d]);
      extent[d] = upper[d] - lower[d];
      TA_ASSERT(extent[d] ==
                static_cast<index1_type>(upper_bound_d - lower_bound_d));
    }
 
    // Set the volume seed
    volume_ = 1ul;
    offset_ = 0l;
 
    // Compute strides, volume, and offset, starting with last (least
    // significant) dimension
    for (int d = int(rank_) - 1; d >= 0; --d) {
      stride[d] = volume_;
      offset_ += lower[d] * stride[d];
      volume_ *= extent[d];
    }
  }
 
  // clang-format off
 
  // clang-format on
  template <typename PairRange,
            typename = std::enable_if_t<detail::is_gpair_range_v<PairRange>>>
  void init_range_data(const PairRange& bounds) {
    // Construct temp pointers
    auto* MADNESS_RESTRICT const lower = data();
    auto* MADNESS_RESTRICT const upper = lower + rank_;
    auto* MADNESS_RESTRICT const extent = upper + rank_;
    auto* MADNESS_RESTRICT const stride = extent + rank_;
 
    // Compute range data
    int d = 0;
    for (auto&& bound_d : bounds) {
      // Compute data for element i of lower, upper, and extent
      const auto lower_bound_d = detail::at(bound_d, 0);
      const auto upper_bound_d = detail::at(bound_d, 1);
      lower[d] = lower_bound_d;
      upper[d] = upper_bound_d;
      // Check input dimensions
      TA_ASSERT(lower[d] <= upper[d]);
      extent[d] = upper[d] - lower[d];
      ++d;
    }
    // Compute strides, volume, and offset, starting with last (least
    // significant) dimension
    volume_ = 1ul;
    offset_ = 0l;
    for (int d = int(rank_) - 1; d >= 0; --d) {
      stride[d] = volume_;
      offset_ += lower[d] * stride[d];
      volume_ *= extent[d];
    }
  }
 
 
  template <typename Index, typename std::enable_if_t<
                                detail::is_integral_range_v<Index>>* = nullptr>
  void init_range_data(const Index& extents) {
    // Construct temp pointers
    auto* MADNESS_RESTRICT const lower = data();
    auto* MADNESS_RESTRICT const upper = lower + rank_;
    auto* MADNESS_RESTRICT const extent = upper + rank_;
    auto* MADNESS_RESTRICT const stride = extent + rank_;
 
    // Initialize extents and bounds
    auto it = std::begin(extents);
    auto end = std::end(extents);
    for (int d = 0; it != end; ++it, ++d) {
      const auto extent_d = *it;
      lower[d] = 0ul;
      upper[d] = extent_d;
      extent[d] = extent_d;
      // Check bounds of the input extent
      TA_ASSERT(extent[d] >= 0);
    }
 
    // Compute strides and volume, starting with last (least significant)
    // dimension
    volume_ = 1ul;
    offset_ = 0ul;
    for (int d = int(rank_) - 1; d >= 0; --d) {
      stride[d] = volume_;
      volume_ *= extent[d];
    }
  }
 
 
  template <typename... Indices,
            typename std::enable_if<
                detail::is_integral_list<Indices...>::value>::type* = nullptr>
  void init_range_data(const std::tuple<Indices...>& extents) {
    const constexpr std::size_t rank =
        std::tuple_size<std::tuple<Indices...>>::value;
    TA_ASSERT(rank_ == rank);
 
    // Set the offset and volume initial values
    volume_ = 1ul;
    offset_ = 0ul;
 
    // initialize by recursion
    init_range_data_helper(extents, std::make_index_sequence<rank>{});
  }
 
  template <typename... Indices, std::size_t... Is>
  void init_range_data_helper(const std::tuple<Indices...>& extents,
                              std::index_sequence<Is...>) {
    int workers[] = {0, (init_range_data_helper_iter<Is>(extents), 0)...};
    ++workers[0];
  }
 
  template <std::size_t I, typename... Indices>
  void init_range_data_helper_iter(const std::tuple<Indices...>& extents) {
    // Get extent i
    const auto extent_i = std::get<I>(extents);
 
    auto* MADNESS_RESTRICT const lower = data();
    auto* MADNESS_RESTRICT const upper = lower + rank_;
    auto* MADNESS_RESTRICT const extent = upper + rank_;
    auto* MADNESS_RESTRICT const stride = extent + rank_;
 
    lower[I] = 0ul;
    upper[I] = extent_i;
    extent[I] = extent_i;
    // Check bounds of the input extent
    TA_ASSERT(extent[I] >= 0ul);
    stride[I] = volume_;
    volume_ *= extent[I];
  }
 
 
  void init_range_data(
      const Permutation& perm,
      const index1_type* MADNESS_RESTRICT const other_lower_bound,
      const index1_type* MADNESS_RESTRICT const other_upper_bound) {
    // Create temporary pointers to this range data
    auto* MADNESS_RESTRICT const lower = data();
    auto* MADNESS_RESTRICT const upper = lower + rank_;
    auto* MADNESS_RESTRICT const extent = upper + rank_;
    auto* MADNESS_RESTRICT const stride = extent + rank_;
 
    // Copy the permuted lower, upper, and extent into this range.
    for (unsigned int i = 0u; i < rank_; ++i) {
      const auto perm_i = perm[i];
 
      // Get the lower bound, upper bound, and extent from other for rank i.
      const auto other_lower_bound_i = other_lower_bound[i];
      const auto other_upper_bound_i = other_upper_bound[i];
      const auto other_extent_i = other_upper_bound_i - other_lower_bound_i;
 
      // Store the permuted lower bound, upper bound, and extent
      lower[perm_i] = other_lower_bound_i;
      upper[perm_i] = other_upper_bound_i;
      extent[perm_i] = other_extent_i;
    }
 
    // Recompute stride, offset, and volume
    volume_ = 1ul;
    offset_ = 0ul;
    for (int i = int(rank_) - 1; i >= 0; --i) {
      const auto lower_i = lower[i];
      const auto extent_i = extent[i];
      stride[i] = volume_;
      offset_ += lower_i * volume_;
      volume_ *= extent_i;
    }
  }
 
 public:
 
  Range() {}
 
 
  template <typename Index1, typename Index2,
            typename std::enable_if_t<
                detail::is_integral_sized_range_v<Index1> &&
                detail::is_integral_sized_range_v<Index2>>* = nullptr>
  Range(const Index1& lower_bound, const Index2& upper_bound) {
    using std::size;
    const auto n = size(lower_bound);
    TA_ASSERT(n == size(upper_bound));
    if (n) {
      // Initialize array memory
      init_datavec(n);
      rank_ = n;
      init_range_data(lower_bound, upper_bound);
    }  // rank-0 Range is null
  }
 
  // clang-format off
 
  // clang-format on
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<std::is_integral_v<Index1> &&
                                        std::is_integral_v<Index2>>>
  Range(const std::initializer_list<Index1>& lower_bound,
        const std::initializer_list<Index2>& upper_bound) {
    using std::size;
    const auto n = size(lower_bound);
    TA_ASSERT(n == size(upper_bound));
    init_datavec(n);
    rank_ = n;
    if (n) {
      init_range_data(lower_bound, upper_bound);
    }  // rank-0 Range is null
  }
 
 
  template <typename Index,
            typename std::enable_if_t<
                detail::is_integral_sized_range_v<Index>>* = nullptr>
  explicit Range(const Index& extent) {
    using std::size;
    const auto n = size(extent);
    if (n) {
      // Initialize array memory
      init_datavec(n);
      rank_ = n;
      init_range_data(extent);
    }  // rank-0 Range is null
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  explicit Range(const std::initializer_list<Index>& extent) {
    using std::size;
    const auto n = size(extent);
    if (n) {
      // Initialize array memory
      init_datavec(n);
      rank_ = n;
      init_range_data(extent);
    }  // rank-0 Range is null
  }
 
  // clang-format off
 
  // clang-format on
  template <typename PairRange,
            typename = std::enable_if_t<detail::is_sized_range_v<PairRange> &&
                                        detail::is_gpair_range_v<PairRange>>>
  explicit Range(const PairRange& bounds) {
    const auto n = std::size(bounds);
    if (n) {
      // Initialize array memory
      init_datavec(n);
      rank_ = n;
      init_range_data(bounds);
    }  // rank-0 Range is null
  }
 
  // clang-format off
 
  // clang-format on
  template <typename GPair>
  explicit Range(const std::initializer_list<GPair>& bounds,
                 std::enable_if_t<detail::is_gpair_v<GPair>>* = nullptr) {
    using std::size;
#ifndef NDEBUG
    if constexpr (detail::is_contiguous_range_v<GPair>) {
      for (auto&& bound_d : bounds) {
        TA_ASSERT(size(bound_d) == 2);
      }
    }
#endif
    const auto n = size(bounds);
    if (n) {
      // Initialize array memory
      init_datavec(n);
      rank_ = n;
      init_range_data(bounds);
    }  // rank-0 Range is null
  }
 
  // clang-format off
 
  // clang-format on
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  explicit Range(
      const std::initializer_list<std::initializer_list<Index>>& bounds) {
    using std::size;
    const auto n = size(bounds);
    if (n) {
#ifndef NDEBUG
      for (auto&& bound_d : bounds) {
        TA_ASSERT(size(bound_d) == 2);
      }
#endif
      // Initialize array memory
      init_datavec(n);
      rank_ = n;
      init_range_data(bounds);
    }  // rank-0 Range is null
  }
 
 
  template <typename... Index, typename std::enable_if<detail::is_integral_list<
                                   Index...>::value>::type* = nullptr>
  explicit Range(const Index... extents)
      : Range(std::array<size_t, sizeof...(Index)>{
            {static_cast<std::size_t>(extents)...}}) {}
 
 
  template <typename... IndexPairs,
            std::enable_if_t<detail::is_integral_pair_list_v<IndexPairs...>>* =
                nullptr>
  explicit Range(const IndexPairs... bounds)
      : Range(std::array<std::pair<std::size_t, std::size_t>,
                         sizeof...(IndexPairs)>{
            {static_cast<std::pair<std::size_t, std::size_t>>(bounds)...}}) {}
 
 
  Range(const Range_& other) = default;
 
 
  Range(Range_&& other)
      : datavec_(std::move(other.datavec_)),
        offset_(other.offset_),
        volume_(other.volume_),
        rank_(other.rank_) {
    // put other into null state
    other.datavec_.clear();
    other.datavec_.shrink_to_fit();
    other.offset_ = 0ul;
    other.volume_ = 0ul;
    other.rank_ = 0u;
  }
 
 
  Range(const Permutation& perm, const Range_& other) {
    TA_ASSERT(perm.size() == other.rank_);
 
    if (other.rank_ > 0ul) {
      rank_ = other.rank_;
 
      if (perm) {
        init_datavec(other.rank_);
        init_range_data(perm, other.lobound_data(), other.upbound_data());
      } else {
        // Simple copy will do
        datavec_ = other.datavec_;
        offset_ = other.offset_;
        volume_ = other.volume_;
      }
    } else  // handle null and rank-0 case
      volume_ = other.volume_;
  }
 
  ~Range() = default;
 
 
  Range_& operator=(const Range_& other) = default;
 
 
  Range_& operator=(Range_&& other) {
    datavec_ = std::move(other.datavec_);
    offset_ = other.offset_;
    volume_ = other.volume_;
    rank_ = other.rank_;
 
    // put other into null state
    other.datavec_.clear();
    other.datavec_.shrink_to_fit();
    other.offset_ = 0l;
    other.volume_ = 0ul;
    other.rank_ = 0u;
 
    return *this;
  }
 
 
  explicit operator bool() const { return rank() != 0; }
 
 
  unsigned int rank() const { return rank_; }
 
 
  std::pair<index1_type, index1_type> dim(std::size_t d) const {
    TA_ASSERT(d < rank());
    return std::make_pair(lobound(d), upbound(d));
  }
 
 
  const index1_type* lobound_data() const { return data(); }
 
 
  index_view_type lobound() const {
    return index_view_type(lobound_data(), rank_);
  }
 
 
  index1_type lobound(size_t dim) const {
    TA_ASSERT(dim < rank_);
    return *(lobound_data() + dim);
  }
 
 
  const index1_type* upbound_data() const { return data() + rank_; }
 
 
  index_view_type upbound() const {
    return index_view_type(upbound_data(), rank_);
  }
 
 
  index1_type upbound(size_t dim) const {
    TA_ASSERT(dim < rank_);
    return *(upbound_data() + dim);
  }
 
 
  const index1_type* extent_data() const { return data() + (rank_ << 1); }
 
 
  index_view_type extent() const {
    return index_view_type(extent_data(), rank_);
  }
 
 
  index1_type extent(size_t dim) const {
    TA_ASSERT(dim < rank_);
    return *(extent_data() + dim);
  }
 
 
  const index1_type* stride_data() const { return extent_data() + rank_; }
 
 
  index_view_type stride() const {
    return index_view_type(stride_data(), rank_);
  }
 
 
  index1_type stride(size_t dim) const {
    TA_ASSERT(dim < rank_);
    return *(stride_data() + dim);
  }
 
 
  ordinal_type volume() const { return volume_; }
 
  ordinal_type area() const { return volume_; }
 
 
  distance_type offset() const { return offset_; }
 
 
  const_iterator begin() const {
    return (volume_ > 0) ? const_iterator(lobound_data(), this) : end();
  }
 
 
  const_iterator end() const { return const_iterator(upbound_data(), this); }
 
 
  template <typename Index,
            typename std::enable_if<detail::is_integral_range_v<Index>,
                                    bool>::type* = nullptr>
  bool includes(const Index& index) const {
    TA_ASSERT(*this);
    const auto* MADNESS_RESTRICT const lower = lobound_data();
    const auto* MADNESS_RESTRICT const upper = upbound_data();
 
    bool result = (rank_ > 0u);
    unsigned int d = 0;
    for (auto&& index_d : index) {
      TA_ASSERT(d < rank_);
      const auto lower_d = lower[d];
      const auto upper_d = upper[d];
      result = result && (index_d >= lower_d) && (index_d < upper_d);
#ifdef NDEBUG
      if (!result) {
        d = rank_;
        break;
      }
#endif
      ++d;
    }
    TA_ASSERT(d == rank_);
 
    return result;
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  bool includes(const std::initializer_list<Index>& index) const {
    return includes<std::initializer_list<Index>>(index);
  }
 
 
  template <typename Ordinal>
  typename std::enable_if<std::is_integral_v<Ordinal>, bool>::type includes(
      Ordinal i) const {
    TA_ASSERT(*this);
    return include_ordinal_(i);
  }
 
  template <typename... Index>
  std::enable_if_t<
      (sizeof...(Index) > 1ul) && (std::is_integral_v<Index> && ...), bool>
  includes(const Index&... index) const {
    const index1_type i[sizeof...(Index)] = {
        static_cast<index1_type>(index)...};
    return includes(i);
  }
 
 
  Range_& operator*=(const Permutation& perm);
 
 
  template <
      typename Index1, typename Index2,
      typename = std::enable_if_t<detail::is_integral_sized_range_v<Index1> &&
                                  detail::is_integral_sized_range_v<Index2>>>
  Range_& resize(const Index1& lower_bound, const Index2& upper_bound) {
    using std::size;
    const auto n = size(lower_bound);
    TA_ASSERT(n == size(upper_bound));
 
    // Reallocate memory for range arrays
    if (rank_ != n) {
      init_datavec(n);
      rank_ = n;
    }
    if (n > 0ul)
      init_range_data(lower_bound, upper_bound);
    else
      volume_ = 0ul;
 
    return *this;
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<detail::is_integral_range_v<Index>>>
  Range_& inplace_shift(const Index& bound_shift) {
    auto* MADNESS_RESTRICT const lower = lobound_data_nc();
    auto* MADNESS_RESTRICT const upper = upbound_data_nc();
    const auto* MADNESS_RESTRICT const stride = upper + rank_ + rank_;
 
    // update the data
    offset_ = 0ul;
    unsigned int d = 0;
    for (auto&& bound_shift_d : bound_shift) {
      TA_ASSERT(d < rank_);
      // Load range data
      auto lower_d = lower[d];
      auto upper_d = upper[d];
      const auto stride_d = stride[d];
 
      // Compute new range bounds
      lower_d += bound_shift_d;
      upper_d += bound_shift_d;
 
      // Update range data
      offset_ += lower_d * stride_d;
      lower[d] = lower_d;
      upper[d] = upper_d;
 
      ++d;
    }
    TA_ASSERT(d == rank_);
 
    return *this;
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  Range_& inplace_shift(const std::initializer_list<Index>& bound_shift) {
    return inplace_shift<std::initializer_list<Index>>(bound_shift);
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<detail::is_integral_range_v<Index>>>
  Range_ shift(const Index& bound_shift) {
    Range_ result(*this);
    result.inplace_shift(bound_shift);
    return result;
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  Range_ shift(const std::initializer_list<Index>& bound_shift) {
    Range_ result(*this);
    result.inplace_shift(bound_shift);
    return result;
  }
 
 
  ordinal_type ordinal(const ordinal_type index) const {
    TA_ASSERT(includes(index));
    return index;
  }
 
 
  template <typename Index, typename std::enable_if_t<
                                detail::is_integral_range_v<Index>>* = nullptr>
  ordinal_type ordinal(const Index& index) const {
    TA_ASSERT(includes(index));
 
    auto* MADNESS_RESTRICT const stride = stride_data();
 
    ordinal_type result = 0ul;
    unsigned int d = 0;
    for (auto&& index_d : index) {
      TA_ASSERT(d < rank_);
      const auto stride_d = stride[d];
      result += index_d * stride_d;
      ++d;
    }
    TA_ASSERT(d == rank_);
 
    return result - offset_;
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  ordinal_type ordinal(const std::initializer_list<Index>& index) const {
    return this->ordinal<std::initializer_list<Index>>(index);
  }
 
 
  template <
      typename... Index,
      typename std::enable_if_t<(sizeof...(Index) > 1ul) &&
                                (std::is_integral_v<Index> && ...)>* = nullptr>
  ordinal_type ordinal(const Index&... index) const {
    const index1_type temp_index[sizeof...(Index)] = {
        static_cast<index1_type>(index)...};
    return ordinal(temp_index);
  }
 
 
  index_type idx(ordinal_type index) const {
    // Check that index is contained by range.
    // N.B. this will fail if any extent is zero
    TA_ASSERT(includes(index));
 
    // Construct result coordinate index object and allocate its memory.
    Range_::index result(rank_, 0);
 
    // Get pointers to the data
    auto* MADNESS_RESTRICT const result_data = result.data();
    const auto* MADNESS_RESTRICT const lower = lobound_data();
    const auto* MADNESS_RESTRICT const size = extent_data();
 
    // Compute the coordinate index of index in range.
    for (int i = int(rank_) - 1; i >= 0; --i) {
      const auto lower_i = lower[i];
      const auto size_i = size[i];
 
      // Compute result index element i
      const auto result_i = (index % size_i) + lower_i;
      index /= size_i;
 
      // Store result
      result_data[i] = result_i;
    }
 
    return result;
  }
 
 
  template <typename Index, typename std::enable_if_t<
                                detail::is_integral_range_v<Index>>* = nullptr>
  const Index& idx(const Index& i) const {
    TA_ASSERT(includes(i));
    return i;
  }
 
  template <typename Archive>
  void serialize(Archive& ar) {
    ar& rank_& datavec_& offset_& volume_;
  }
 
  void swap(Range_& other) {
    // Get temp data
    std::swap(datavec_, other.datavec_);
    std::swap(offset_, other.offset_);
    std::swap(volume_, other.volume_);
    std::swap(rank_, other.rank_);
  }
 
 private:
 
  template <typename Index>
  typename std::enable_if<std::is_integral_v<Index> && std::is_signed_v<Index>,
                          bool>::type
  include_ordinal_(Index i) const {
    return (i >= Index(0)) && (i < Index(volume_));
  }
 
 
  template <typename Index>
  typename std::enable_if<
      std::is_integral_v<Index> && !std::is_signed<Index>::value, bool>::type
  include_ordinal_(Index i) const {
    return i < volume_;
  }
 
 
  void increment(index_type& i) const {
    TA_ASSERT(includes(i));
 
    const auto* MADNESS_RESTRICT const lower = lobound_data();
    const auto* MADNESS_RESTRICT const upper = upbound_data();
 
    for (int d = int(rank_) - 1; d >= 0; --d) {
      // increment coordinate
      ++i[d];
 
      // break if done
      if (i[d] < upper[d]) return;
 
      // Reset current index to lower bound.
      i[d] = lower[d];
    }
 
    // if the current location was set to lower then it was at the end and
    // needs to be reset to equal upper.
    std::copy(upper, upper + rank_, i.begin());
  }
 
 
  void advance(index_type& i, std::ptrdiff_t n) const {
    TA_ASSERT(includes(i));
    const auto o = ordinal(i) + n;
    TA_ASSERT(includes(o));
    i = idx(o);
  }
 
 
  std::ptrdiff_t distance_to(const index_type& first,
                             const index_type& last) const {
    TA_ASSERT(includes(first));
    TA_ASSERT(includes(last));
    return ordinal(last) - ordinal(first);
  }
 
};  // class Range
 
inline Range& Range::operator*=(const Permutation& perm) {
  TA_ASSERT(perm.size() == rank_);
  if (rank_ > 1ul) {
    // Copy the lower and upper bound data into a temporary array
    container::svector<index1_type, 2 * TA_MAX_SOO_RANK_METADATA> temp_lower(
        rank_ << 1);
    const auto* MADNESS_RESTRICT const temp_upper = temp_lower.data() + rank_;
    std::copy(lobound_data(), lobound_data() + (rank_ << 1), temp_lower.data());
 
    init_range_data(perm, temp_lower.data(), temp_upper);
  }
  return *this;
}
 
inline void swap(Range& r0, Range& r1) {  // no throw
  r0.swap(r1);
}
 
 
inline Range operator*(const Permutation& perm, const Range& r) {
  return Range(perm, r);
}
 
 
template <typename I, typename = std::enable_if_t<std::is_integral_v<I>>>
Range permute(const Range& r, std::initializer_list<I> perm) {
  return Permutation(perm) * r;
}
 
 
inline bool operator==(const Range& r1, const Range& r2) {
  return (r1.rank() == r2.rank()) &&
         !std::memcmp(r1.lobound_data(), r2.lobound_data(),
                      r1.rank() * (2u * sizeof(Range::index1_type)));
}
 
 
inline bool operator!=(const Range& r1, const Range& r2) {
  return (r1.rank() != r2.rank()) ||
         std::memcmp(r1.lobound_data(), r2.lobound_data(),
                     r1.rank() * (2u * sizeof(Range::index1_type)));
}
 
 
inline std::ostream& operator<<(std::ostream& os, const Range& r) {
  os << "[ ";
  detail::print_array(os, r.lobound_data(), r.rank());
  os << ", ";
  detail::print_array(os, r.upbound_data(), r.rank());
  os << " )";
  return os;
}
 
 
inline bool is_congruent(const Range& r1, const Range& r2) {
  return (r1.rank() == r2.rank()) &&
         std::equal(r1.extent_data(), r1.extent_data() + r1.rank(),
                    r2.extent_data());
}
 
 
inline bool is_contiguous(const Range& range) { return true; }
 
}  // namespace TiledArray
#endif  // TILEDARRAY_RANGE_H__INCLUDED