array_impl.h

Go to the documentation of this file.

/*
 *  This file is a part of TiledArray.
 *  Copyright (C) 2014  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
 *
 *  array_impl.h
 *  Oct 24, 2014
 *
 */
 
#ifndef TILEDARRAY_ARRAY_IMPL_H__INCLUDED
#define TILEDARRAY_ARRAY_IMPL_H__INCLUDED
 
#include <TiledArray/distributed_storage.h>
#include <TiledArray/tensor_impl.h>
#include <TiledArray/transform_iterator.h>
#include <TiledArray/type_traits.h>
 
namespace TiledArray {
namespace detail {
 
// Forward declaration
template <typename>
class TensorReference;
template <typename>
class TensorConstReference;
template <typename, typename>
class ArrayIterator;
 
 
template <typename Impl>
class TileReference {
 private:
  template <typename, typename>
  friend class ArrayIterator;
 
  template <typename>
  friend class TileConstReference;
 
  typedef typename Impl::range_type range_type;
  typedef typename Impl::range_type::index index_type;
  typedef typename Impl::ordinal_type ordinal_type;
 
  Impl* tensor_;        
  ordinal_type index_;  
 
  // Not allowed
  TileReference<Impl>& operator=(const TileReference<Impl>&);
 
 public:
  TileReference(Impl* tensor, const typename Impl::ordinal_type index)
      : tensor_(tensor), index_(index) {}
 
  TileReference(const TileReference<Impl>& other)
      : tensor_(other.tensor_), index_(other.index_) {}
 
  template <typename Value>
  TileReference<Impl>& operator=(const Value& value) {
    tensor_->set(index_, value);
    return *this;
  }
 
  typename Impl::future future() const {
    TA_ASSERT(tensor_);
    return tensor_->get(index_);
  }
 
  typename Impl::value_type get() const {
    // NOTE: return by value to avoid lifetime issues.
    TA_ASSERT(tensor_);
    return future().get();
  }
 
  operator typename Impl::future() const { return this->future(); }
 
  operator typename Impl::value_type() const { return get(); }
 
 
  index_type index() const {
    TA_ASSERT(tensor_);
    return tensor_->tiles_range().idx(index_);
  }
 
 
  ordinal_type ordinal() const { return index_; }
 
 
  range_type make_range() const {
    TA_ASSERT(tensor_);
    return tensor_->trange().make_tile_range(index_);
  }
 
};  // class TileReference
 
template <typename Impl>
bool operator==(const TileReference<Impl>& a, const TileReference<Impl>& b) {
  return a.get() == b.get();
}
 
template <typename Impl>
bool operator!=(const TileReference<Impl>& a, const TileReference<Impl>& b) {
  return !(a == b);
}
 
template <typename Impl>
std::ostream& operator<<(std::ostream& os, const TileReference<Impl>& a) {
  os << a.get();
  return os;
}
 
 
template <typename Impl>
class TileConstReference {
 private:
  template <typename, typename>
  friend class ArrayIterator;
 
  const Impl* tensor_;  
  typename Impl::ordinal_type index_;  
 
  // Not allowed
  TileConstReference<Impl>& operator=(const TileConstReference<Impl>&);
 
 public:
  TileConstReference(const Impl* tensor,
                     const typename Impl::ordinal_type index)
      : tensor_(tensor), index_(index) {}
 
  TileConstReference(const TileConstReference<Impl>& other)
      : tensor_(other.tensor_), index_(other.index_) {}
 
  TileConstReference(const TileReference<Impl>& other)
      : tensor_(other.tensor_), index_(other.index_) {}
 
  typename Impl::future future() const {
    TA_ASSERT(tensor_);
    return tensor_->get(index_);
  }
 
  typename Impl::value_type get() const {
    // NOTE: return by value to avoid lifetime issues.
    TA_ASSERT(tensor_);
    return future().get();
  }
 
  operator typename Impl::future() const { return tensor_->get(index_); }
 
  operator typename Impl::value_type() const { return get(); }
};  // class TileConstReference
 
template <typename Impl>
bool operator==(const TileConstReference<Impl>& a,
                const TileConstReference<Impl>& b) {
  return a.get() == b.get();
}
 
template <typename Impl>
bool operator!=(const TileConstReference<Impl>& a,
                const TileConstReference<Impl>& b) {
  return !(a == b);
}
 
template <typename Impl>
std::ostream& operator<<(std::ostream& os, const TileConstReference<Impl>& a) {
  os << a.get();
  return os;
}
 
}  // namespace detail
}  // namespace TiledArray
 
namespace madness {
namespace detail {
 
// The following class specializations are required so MADNESS will do the
// right thing when given a TileReference or TileConstReference object
// as an input for task functions.
 
template <typename Impl>
struct task_arg<TiledArray::detail::TileReference<Impl>> {
  typedef typename Impl::value_type type;
  typedef typename Impl::future holderT;
};  // struct task_arg<TiledArray::detail::TileReference<Impl> >
 
template <typename Impl>
struct task_arg<TiledArray::detail::TileConstReference<Impl>> {
  typedef typename Impl::value_type type;
  typedef typename Impl::future holderT;
};  // struct task_arg<TiledArray::detail::TileConstReference<Impl> >
 
}  // namespace detail
}  // namespace madness
 
namespace TiledArray {
namespace detail {
 
 
template <typename Impl, typename Reference>
class ArrayIterator {
 private:
  // Give access to other iterator types.
  template <typename, typename>
  friend class ArrayIterator;
 
  Impl* array_;
  typename Impl::pmap_interface::const_iterator it_;
 
 public:
  typedef ptrdiff_t difference_type;  
  typedef
      typename Impl::future value_type;  
  typedef PointerProxy<value_type> pointer;  
  typedef Reference reference;  
  typedef std::forward_iterator_tag
      iterator_category;  
  typedef ArrayIterator<Impl, Reference> ArrayIterator_;  
  typedef typename Impl::range_type::index index_type;
  typedef typename Impl::ordinal_type ordinal_type;
  typedef typename Impl::range_type range_type;
  typedef typename Impl::value_type tile_type;
 
 private:
  void advance() {
    TA_ASSERT(array_);
    const typename Impl::pmap_interface::const_iterator end =
        array_->pmap()->end();
    do {
      ++it_;
    } while ((it_ != end) && array_->is_zero(*it_));
  }
 
 public:
  ArrayIterator() : array_(NULL), it_() {}
 
  ArrayIterator(Impl* tensor, typename Impl::pmap_interface::const_iterator it)
      : array_(tensor), it_(it) {}
 
 
  ArrayIterator(const ArrayIterator_& other)
      : array_(other.array_), it_(other.it_) {}
 
 
  template <
      typename I, typename R,
      typename std::enable_if<!((!std::is_const<Impl>::value) &&
                                std::is_const<I>::value)>::type* = nullptr>
  ArrayIterator(const ArrayIterator<I, R>& other)
      : array_(other.array_), it_(other.it_) {}
 
 
  ArrayIterator_& operator=(const ArrayIterator_& other) {
    array_ = other.array_;
    it_ = other.it_;
 
    return *this;
  }
 
 
  template <typename R>
  ArrayIterator_& operator=(const ArrayIterator<Impl, R>& other) {
    array_ = other.array_;
    it_ = other.it_;
 
    return *this;
  }
 
 
  ArrayIterator_& operator++() {
    advance();
    return *this;
  }
 
 
  ArrayIterator_ operator++(int) {
    ArrayIterator_ tmp(*this);
    advance();
    return tmp;
  }
 
 
  template <typename I, typename R>
  bool operator==(const ArrayIterator<I, R>& other) const {
    return (array_ == other.array_) && (it_ == other.it_);
  }
 
 
  template <typename I, typename R>
  bool operator!=(const ArrayIterator<I, R>& other) const {
    return (array_ != other.array_) || (it_ != other.it_);
  }
 
 
  reference operator*() const {
    TA_ASSERT(array_);
    return reference(array_, *it_);
  }
 
 
  pointer operator->() const {
    TA_ASSERT(array_);
    return pointer(array_->get(*it_));
  }
 
 
  index_type index() const {
    TA_ASSERT(array_);
    return array_->tiles_range().idx(*it_);
  }
 
 
  ordinal_type ordinal() const {
    TA_ASSERT(array_);
    return *it_;
  }
 
 
  range_type make_range() const {
    TA_ASSERT(array_);
    TA_ASSERT(it_ != array_->pmap()->end());
    return array_->trange().make_tile_range(*it_);
  }
 
};  // class TensorIterator
 
 
template <typename Tile, typename Policy>
class ArrayImpl : public TensorImpl<Policy> {
 public:
  typedef ArrayImpl<Tile, Policy> ArrayImpl_;  
  typedef TensorImpl<Policy> TensorImpl_;  
  typedef typename TensorImpl_::index1_type index1_type;    
  typedef typename TensorImpl_::ordinal_type ordinal_type;  
  typedef typename TensorImpl_::policy_type
      policy_type;  
  typedef typename TensorImpl_::trange_type
      trange_type;  
  typedef typename TensorImpl_::range_type
      range_type;  
  typedef typename TensorImpl_::shape_type shape_type;  
  typedef typename TensorImpl_::pmap_interface
      pmap_interface;       
  typedef Tile value_type;  
  typedef
      typename eval_trait<Tile>::type eval_type;  
  typedef typename numeric_type<value_type>::type
      numeric_type;  
  typedef DistributedStorage<value_type>
      storage_type;                              
  typedef typename storage_type::future future;  
  typedef TileReference<ArrayImpl_> reference;   
  typedef TileConstReference<ArrayImpl_>
      const_reference;  
  typedef ArrayIterator<ArrayImpl_, reference> iterator;  
  typedef ArrayIterator<const ArrayImpl_, const_reference>
      const_iterator;  
 
 private:
  storage_type data_;  
 
 public:
 
  ArrayImpl(World& world, const trange_type& trange, const shape_type& shape,
            const std::shared_ptr<pmap_interface>& pmap)
      : TensorImpl_(world, trange, shape, pmap),
        data_(world, trange.tiles_range().volume(), pmap) {}
 
  virtual ~ArrayImpl() {}
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index> ||
                                        detail::is_integral_range_v<Index>>>
  future get(const Index& i) const {
    TA_ASSERT(!TensorImpl_::is_zero(i));
    return data_.get(TensorImpl_::trange().tiles_range().ordinal(i));
  }
 
 
  template <typename Integer,
            typename = std::enable_if_t<std::is_integral_v<Integer>>>
  future get(const std::initializer_list<Integer>& i) const {
    return get<std::initializer_list<Integer>>(i);
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index> ||
                                        detail::is_integral_range_v<Index>>>
  const future& get_local(const Index& i) const {
    TA_ASSERT(!TensorImpl_::is_zero(i) && TensorImpl_::is_local(i));
    return data_.get_local(TensorImpl_::trange().tiles_range().ordinal(i));
  }
 
 
  template <typename Integer,
            typename = std::enable_if_t<std::is_integral_v<Integer>>>
  const future& get_local(const std::initializer_list<Integer>& i) const {
    return get_local<std::initializer_list<Integer>>(i);
  }
 
 
  template <typename Index,
            typename = std::enable_if_t<std::is_integral_v<Index> ||
                                        detail::is_integral_range_v<Index>>>
  future& get_local(const Index& i) {
    TA_ASSERT(!TensorImpl_::is_zero(i) && TensorImpl_::is_local(i));
    return data_.get_local(TensorImpl_::trange().tiles_range().ordinal(i));
  }
 
 
  template <typename Integer,
            typename = std::enable_if_t<std::is_integral_v<Integer>>>
  future& get_local(const std::initializer_list<Integer>& i) {
    return get_local<std::initializer_list<Integer>>(i);
  }
 
 
  template <typename Index, typename Value,
            typename = std::enable_if_t<std::is_integral_v<Index> ||
                                        detail::is_integral_range_v<Index>>>
  void set(const Index& i, Value&& value) {
    TA_ASSERT(!TensorImpl_::is_zero(i));
    const auto ord = TensorImpl_::trange().tiles_range().ordinal(i);
    data_.set(ord, std::forward<Value>(value));
    if (set_notifier_accessor()) {
      set_notifier_accessor()(*this, ord);
    }
  }
 
 
  template <typename Index, typename Value,
            typename = std::enable_if_t<std::is_integral_v<Index>>>
  void set(const std::initializer_list<Index>& i, Value&& value) {
    TA_ASSERT(!TensorImpl_::is_zero(i));
    const auto ord = TensorImpl_::trange().tiles_range().ordinal(i);
    data_.set(ord, std::forward<Value>(value));
    if (set_notifier_accessor()) {
      set_notifier_accessor()(*this, ord);
    }
  }
 
 
  iterator begin() {
    // Get the pmap iterator
    typename pmap_interface::const_iterator it = TensorImpl_::pmap()->begin();
 
    // Find the first non-zero iterator
    const typename pmap_interface::const_iterator end =
        TensorImpl_::pmap()->end();
    while ((it != end) && TensorImpl_::is_zero(*it)) ++it;
 
    // Construct and return the iterator
    return iterator(this, it);
  }
 
 
  const_iterator cbegin() const {
    // Get the pmap iterator
    typename pmap_interface::const_iterator it = TensorImpl_::pmap()->begin();
 
    // Find the fist non-zero iterator
    const typename pmap_interface::const_iterator end =
        TensorImpl_::pmap()->end();
    while ((it != end) && TensorImpl_::is_zero(*it)) ++it;
 
    // Construct and return the iterator
    return const_iterator(this, it);
  }
 
 
  iterator end() { return iterator(this, TensorImpl_::pmap()->end()); }
 
 
  const_iterator cend() const {
    return const_iterator(this, TensorImpl_::pmap()->end());
  }
 
 
  const madness::uniqueidT& id() const { return data_.id(); }
 
  static std::function<void(const ArrayImpl_&, int64_t)>&
  set_notifier_accessor() {
    static std::function<void(const ArrayImpl_&, int64_t)> value;
    return value;
  }
 
};  // class ArrayImpl
 
#ifndef TILEDARRAY_HEADER_ONLY
 
extern template class ArrayImpl<
    Tensor<double, Eigen::aligned_allocator<double>>, DensePolicy>;
extern template class ArrayImpl<Tensor<float, Eigen::aligned_allocator<float>>,
                                DensePolicy>;
extern template class ArrayImpl<Tensor<int, Eigen::aligned_allocator<int>>,
                                DensePolicy>;
extern template class ArrayImpl<Tensor<long, Eigen::aligned_allocator<long>>,
                                DensePolicy>;
//    extern template
//    class ArrayImpl<Tensor<std::complex<double>,
//    Eigen::aligned_allocator<std::complex<double> > >, DensePolicy>; extern
//    template class ArrayImpl<Tensor<std::complex<float>,
//    Eigen::aligned_allocator<std::complex<float> > >, DensePolicy>;
 
extern template class ArrayImpl<
    Tensor<double, Eigen::aligned_allocator<double>>, SparsePolicy>;
extern template class ArrayImpl<Tensor<float, Eigen::aligned_allocator<float>>,
                                SparsePolicy>;
extern template class ArrayImpl<Tensor<int, Eigen::aligned_allocator<int>>,
                                SparsePolicy>;
extern template class ArrayImpl<Tensor<long, Eigen::aligned_allocator<long>>,
                                SparsePolicy>;
//    extern template
//    class ArrayImpl<Tensor<std::complex<double>,
//    Eigen::aligned_allocator<std::complex<double> > >, SparsePolicy>; extern
//    template class ArrayImpl<Tensor<std::complex<float>,
//    Eigen::aligned_allocator<std::complex<float> > >, SparsePolicy>;
 
#endif  // TILEDARRAY_HEADER_ONLY
 
}  // namespace detail
}  // namespace TiledArray
 
#endif  // TILEDARRAY_ARRAY_IMPL_H__INCLUDED