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/tensor_impl.h>
#include <TiledArray/distributed_storage.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::size_type ordinal_type;

      Impl* tensor_; 
      ordinal_type index_; 

      // Not allowed
      TileReference<Impl>& operator=(const TileReference<Impl>&);
    public:

      TileReference(Impl* tensor, const typename Impl::size_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::size_type index_; 

      // Not allowed
      TileConstReference<Impl>& operator=(const TileConstReference<Impl>&);
    public:

      TileConstReference(const Impl* tensor, const typename Impl::size_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::size_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_::size_type size_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>
      future get(const Index& i) const {
        TA_ASSERT(! TensorImpl_::is_zero(i));
        return data_.get(TensorImpl_::trange().tiles_range().ordinal(i));
      }


      template <typename Integer>
      future get(const std::initializer_list<Integer>& i) const {
        return get<std::initializer_list<Integer>>(i);
      }


      template <typename Index, typename Value>
      void set(const Index& i, const Value& value) {
        TA_ASSERT(! TensorImpl_::is_zero(i));
        data_.set(TensorImpl_::trange().tiles_range().ordinal(i), value);
      }


      iterator begin() {
        // 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 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(); }

    }; // 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