block_range.h

Go to the documentation of this file.

/*
 *  This file is a part of TiledArray.
 *  Copyright (C) 2015  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/>.
 *
 *  Justus Calvin
 *  Department of Chemistry, Virginia Tech
 *
 *  block_range.h
 *  May 29, 2015
 *
 */
 
#ifndef TILEDARRAY_BLOCK_RANGE_H__INCLUDED
#define TILEDARRAY_BLOCK_RANGE_H__INCLUDED
 
#include <TiledArray/range.h>
 
namespace TiledArray {
 
class BlockRange : public Range {
 private:
  using Range::datavec_;
  using Range::offset_;
  using Range::rank_;
  using Range::volume_;
 
  Range::ordinal_type block_offset_ = 0ul;
 
  friend inline bool operator==(const BlockRange& r1, const BlockRange& r2);
 
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<detail::is_integral_range_v<Index1> &&
                                        detail::is_integral_range_v<Index2>>>
  void init(const Range& range, const Index1& lower_bound,
            const Index2& upper_bound) {
    TA_ASSERT(range.rank());
 
    // Initialize the block range data members
    init_datavec(range.rank());
    offset_ = range.offset();
    volume_ = 1ul;
    rank_ = range.rank();
    block_offset_ = 0ul;
 
    // Construct temp pointers
    const auto* MADNESS_RESTRICT const range_stride = range.stride_data();
    auto* MADNESS_RESTRICT const lower = lobound_data_nc();
    auto* MADNESS_RESTRICT const upper = upbound_data_nc();
    auto* MADNESS_RESTRICT const extent = extent_data_nc();
    auto* MADNESS_RESTRICT const stride = stride_data_nc();
 
    // 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 i of lower, upper, and extent
      const auto lower_bound_d = *lower_it;
      const auto upper_bound_d = *upper_it;
 
      // Set the block range data
      lower[d] = lower_bound_d;
      upper[d] = upper_bound_d;
      // Check input dimensions
      TA_ASSERT(lower[d] >= range.lobound(d));
      TA_ASSERT(lower[d] < upper[d]);
      TA_ASSERT(upper[d] <= range.upbound(d));
      extent[d] = upper[d] - lower[d];
      TA_ASSERT(extent[d] ==
                static_cast<index1_type>(upper_bound_d - lower_bound_d));
    }
 
    // Compute strides, volume, and offset, starting with last (least
    // significant) dimension
    for (int d = int(rank_) - 1; d >= 0; --d) {
      const auto range_stride_d = range_stride[d];
      stride[d] = range_stride_d;
      block_offset_ += lower[d] * range_stride_d;
      volume_ *= extent[d];
    }
  }
 
  template <typename PairRange,
            typename = std::enable_if_t<detail::is_gpair_range_v<PairRange>>>
  void init(const Range& range, const PairRange& bounds) {
    TA_ASSERT(range.rank());
 
    // Initialize the block range data members
    init_datavec(range.rank());
    offset_ = range.offset();
    volume_ = 1ul;
    rank_ = range.rank();
    block_offset_ = 0ul;
 
    // Construct temp pointers
    const auto* MADNESS_RESTRICT const range_stride = range.stride_data();
    auto* MADNESS_RESTRICT const lower = lobound_data_nc();
    auto* MADNESS_RESTRICT const upper = upbound_data_nc();
    auto* MADNESS_RESTRICT const extent = extent_data_nc();
    auto* MADNESS_RESTRICT const stride = stride_data_nc();
 
    // 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] >= range.lobound(d));
      TA_ASSERT(lower[d] < upper[d]);
      TA_ASSERT(upper[d] <= range.upbound(d));
      extent[d] = upper[d] - lower[d];
      TA_ASSERT(extent[d] ==
                static_cast<index1_type>(upper_bound_d - lower_bound_d));
      ++d;
    }
 
    // Compute strides, volume, and offset, starting with last (least
    // significant) dimension
    for (int d = int(rank_) - 1; d >= 0; --d) {
      const auto range_stride_d = range_stride[d];
      stride[d] = range_stride_d;
      block_offset_ += lower[d] * range_stride_d;
      volume_ *= extent[d];
    }
  }
 
 public:
  // Compiler generated functions
  BlockRange() = default;
  BlockRange(const BlockRange&) = default;
  BlockRange(BlockRange&&) = default;
  ~BlockRange() = default;
  BlockRange& operator=(const BlockRange&) = default;
  BlockRange& operator=(BlockRange&&) = default;
 
  // clang-format off
 
  // clang-format on
  template <typename Index1, typename Index2,
            typename = std::enable_if_t<detail::is_integral_range_v<Index1> &&
                                        detail::is_integral_range_v<Index2>>>
  BlockRange(const Range& range, const Index1& lower_bound,
             const Index2& upper_bound)
      : Range() {
    init(range, lower_bound, upper_bound);
  }
 
  // 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>>>
  BlockRange(const Range& range,
             const std::initializer_list<Index1>& lower_bound,
             const std::initializer_list<Index2>& upper_bound)
      : Range() {
    init(range, lower_bound, upper_bound);
  }
 
  // clang-format off
 
  // clang-format on
  template <typename PairRange,
            typename = std::enable_if_t<detail::is_gpair_range_v<PairRange>>>
  BlockRange(const Range& range, const PairRange& bounds) : Range() {
    init(range, bounds);
  }
 
  // clang-format off
 
  // clang-format on
  template <typename GPair>
  BlockRange(const Range& range, const std::initializer_list<GPair>& bounds,
             std::enable_if_t<detail::is_gpair_v<GPair>>* = nullptr)
      : Range() {
#ifndef NDEBUG
    if constexpr (detail::is_contiguous_range_v<GPair>) {
      for (auto&& bound_d : bounds) {
        TA_ASSERT(size(bound_d) == 2);
      }
    }
#endif
    init<std::initializer_list<GPair>>(range, bounds);
  }
 
  // clang-format off
 
  // clang-format on
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  BlockRange(const Range& range,
             const std::initializer_list<std::initializer_list<Index>>& bounds)
      : Range() {
#ifndef NDEBUG
    for (auto&& bound_d : bounds) {
      TA_ASSERT(size(bound_d) == 2);
    }
#endif
    init<std::initializer_list<std::initializer_list<Index>>>(range, bounds);
  }
 
 
  template <typename Index, typename std::enable_if_t<
                                detail::is_integral_range_v<Index>>* = nullptr>
  ordinal_type ordinal(const Index& index) const {
    return Range::ordinal(index);
  }
 
  template <typename... Index,
            typename std::enable_if<(sizeof...(Index) > 1ul)>::type* = nullptr>
  ordinal_type ordinal(const Index&... index) const {
    return Range::ordinal(index...);
  }
 
 
  ordinal_type ordinal(ordinal_type index) const {
    // Check that index is contained by range.
    TA_ASSERT(includes(index));
 
    // Construct result coordinate index object and allocate its memory.
    ordinal_type result = 0ul;
 
    // Get pointers to the data
    const auto* MADNESS_RESTRICT const size = extent_data();
    const auto* MADNESS_RESTRICT const stride = stride_data();
 
    // Compute the coordinate index of o in range.
    for (int i = int(rank_) - 1; i >= 0; --i) {
      const auto size_i = size[i];
      const auto stride_i = stride[i];
 
      // Compute result index element i
      result += (index % size_i) * stride_i;
      index /= size_i;
    }
 
    return result + block_offset_ - offset_;
  }
 
  template <typename Index>
  BlockRange& resize(const Index&, const Index&) {
    // This function is here to shadow the base class resize function
    TA_EXCEPTION("BlockRange::resize() is not supported");
    return *this;
  }
 
 
  template <typename Index>
  Range_& inplace_shift(const Index&) {
    TA_EXCEPTION("BlockRange::inplace_shift() is not supported");
    return *this;
  }
 
 
  template <typename Index>
  Range_ shift(const Index&) {
    TA_EXCEPTION("BlockRange::shift() is not supported");
    return *this;
  }
 
  void swap(BlockRange& other) {
    Range::swap(other);
    std::swap(block_offset_, other.block_offset_);
  }
 
  template <typename Archive>
  void serialize(const Archive& ar) const {
    Range::serialize(ar);
    ar& block_offset_;
  }
};  // BlockRange
 
// clang-format off
 
// clang-format on
inline bool is_congruent(const BlockRange& r1, const BlockRange& r2) {
  return is_congruent(static_cast<const Range&>(r1),
                      static_cast<const Range&>(r2));
}
 
// clang-format off
 
// clang-format on
inline bool is_congruent(const BlockRange& r1, const Range& r2) {
  return is_congruent(static_cast<const Range&>(r1), r2);
}
 
// clang-format off
 
// clang-format on
inline bool is_congruent(const Range& r1, const BlockRange& r2) {
  return is_congruent(r2, r1);
}
 
 
inline bool operator==(const BlockRange& r1, const BlockRange& r2) {
  return r1.block_offset_ == r2.block_offset_ &&
         static_cast<const Range&>(r1) == static_cast<const Range&>(r2);
}
 
 
inline bool is_contiguous(const BlockRange& range) { return true; }
 
}  // namespace TiledArray
 
#endif  // TILEDARRAY_BLOCK_RANGE_H__INCLUDED