tensor_interface.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
 *
 *  tensor_view.h
 *  May 29, 2015
 *
 */

#ifndef TILEDARRAY_TENSOR_TENSOR_VIEW_H__INCLUDED
#define TILEDARRAY_TENSOR_TENSOR_VIEW_H__INCLUDED

#include <TiledArray/tensor/kernels.h>
#include <TiledArray/tensor/complex.h>

namespace Eigen {

  // Forward declarations
  template <typename> class aligned_allocator;

} // namespace Eigen

namespace TiledArray {

  // Forward declarations
  class Permutation;
  template <typename, typename> class Tensor;
  class Range;
  namespace detail {
    template <typename, typename>
    class TensorInterface;
  }
  template <typename T, typename Index>
  void remap(detail::TensorInterface<T, Range> &, T* const, const Index&, const Index&);
  template <typename T, typename Index>
  void remap(detail::TensorInterface<const T, Range> &, T* const,
          const Index&, const Index&);
  template <typename T>
  void remap(detail::TensorInterface<T, Range> &, T* const,
          const std::initializer_list<std::size_t>&,
          const std::initializer_list<std::size_t>&);
  template <typename T>
  void remap(detail::TensorInterface<const T, Range> &, T* const,
      const std::initializer_list<std::size_t>&,
      const std::initializer_list<std::size_t>&);

  namespace detail {


    template <typename T, typename R>
    class TensorInterface {
    public:
      typedef TensorInterface<T, R> TensorInterface_; 
      typedef R range_type; 
      typedef typename range_type::size_type size_type; 
      typedef typename std::remove_const<T>::type
          value_type; 
      typedef typename std::add_lvalue_reference<T>::type
          reference; 
      typedef typename std::add_lvalue_reference<typename std::add_const<T>::type>::type
          const_reference; 
      typedef typename std::add_pointer<T>::type
          pointer; 
      typedef typename std::add_pointer<typename std::add_const<T>::type>::type
          const_pointer; 
      typedef typename std::ptrdiff_t difference_type; 
      typedef typename detail::numeric_type<value_type>::type
          numeric_type; 
      typedef typename detail::scalar_type<value_type>::type
          scalar_type; 

      typedef Tensor<T, Eigen::aligned_allocator<T> > result_tensor;

    private:

      template <typename X>
      using numeric_t = typename TiledArray::detail::numeric_type<X>::type;

      template <typename, typename>
      friend class TensorInterface;

      template <typename U, typename Index>
      friend void TiledArray::remap(detail::TensorInterface<U, Range>&,
                                    U* const, const Index&, const Index&);

      template <typename U, typename Index>
      friend void TiledArray::remap(detail::TensorInterface<const U, Range>&,
                                    U* const, const Index&, const Index&);

      template <typename U>
      friend void TiledArray::remap(detail::TensorInterface<U, Range> &,
              U* const,
              const std::initializer_list<std::size_t>&,
              const std::initializer_list<std::size_t>&);

      template <typename U>
      friend void TiledArray::remap(detail::TensorInterface<const U, Range>&,
                                    U* const,
                                    const std::initializer_list<std::size_t>&,
                                    const std::initializer_list<std::size_t>&);

      range_type range_;  
      pointer data_;  

    public:

      TensorInterface() = delete;
      ~TensorInterface() = default;
      TensorInterface(const TensorInterface_&) = default;
      TensorInterface(TensorInterface_&&) = default;
      TensorInterface_& operator=(const TensorInterface_&) = delete;
      TensorInterface_& operator=(TensorInterface_&&) = delete;


      template <typename U,
          typename std::enable_if<std::is_convertible<
              typename TensorInterface<U, R>::pointer, pointer>::value
          >::type* = nullptr>
      TensorInterface(const TensorInterface<U, R>& other) :
        range_(other.range_), data_(other.data_)
      { }


      template <typename U,
          typename std::enable_if<std::is_convertible<
              typename TensorInterface<U, R>::pointer, pointer>::value
          >::type* = nullptr>
      TensorInterface(TensorInterface<U, R>&& other) :
        range_(std::move(other.range_)), data_(other.data_)
      {
        other.data_ = nullptr;
      }


      TensorInterface(const range_type& range, pointer data) :
        range_(range), data_(data)
      {
        TA_ASSERT(data);
      }


      TensorInterface(range_type&& range, pointer data) :
        range_(std::move(range)), data_(data)
      {
        TA_ASSERT(data);
      }

      template <typename T1,
          typename std::enable_if<detail::is_tensor<T1>::value>::type* = nullptr>
      TensorInterface_& operator=(const T1& other) {
        TA_ASSERT(data_ != other.data());

        detail::inplace_tensor_op([] (numeric_type& MADNESS_RESTRICT result,
            const numeric_t<T1> arg)
            { result = arg; }, *this, other);

        return *this;
      }


      const range_type& range() const { return range_; }


      size_type size() const { return range_.volume(); }



      pointer data() const { return data_; }


      const_reference operator[](const size_type index) const {
        TA_ASSERT(range_.includes(index));
        return data_[range_.ordinal(index)];
      }


      reference operator[](const size_type index) {
        TA_ASSERT(range_.includes(index));
        return data_[range_.ordinal(index)];
      }



      template<typename... Index>
      reference operator()(const Index&... idx) {
        TA_ASSERT(range_.includes(idx...));
        return data_[range_.ordinal(idx...)];
      }


      template<typename... Index>
      const_reference operator()(const Index&... idx) const {
        TA_ASSERT(range_.includes(idx...));
        return data_[range_.ordinal(idx...)];
      }


      constexpr bool empty() const { return false; }


      void swap(TensorInterface_& other) {
        range_.swap(other.range_);
        std::swap(data_, other.data_);
      }



      template <typename Index>
      TensorInterface_& shift_to(const Index& bound_shift) {
        range_.inplace_shift(bound_shift);
        return *this;
      }


      template <typename Index>
      result_tensor shift(const Index& bound_shift) const {
        return result_tensor(range_.shift(bound_shift), data_);
      }

      // Generic vector operations


      template <typename Right, typename Op,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      result_tensor binary(const Right& right, Op&& op) const {
        return result_tensor(*this, right, op);
      }


      template <typename Right, typename Op,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      result_tensor binary(const Right& right, Op&& op, const Permutation& perm) const {
        return result_tensor(*this, right, op, perm);
      }


      template <typename Right, typename Op,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      TensorInterface_& inplace_binary(const Right& right, Op&& op) {
        detail::inplace_tensor_op(op, *this, right);
        return *this;
      }


      template <typename Op>
      result_tensor unary(Op&& op) const {
        return result_tensor(*this, op);
      }


      template <typename Op>
      result_tensor unary(Op&& op, const Permutation& perm) const {
        return result_tensor(*this, op, perm);
      }


      template <typename Op>
      TensorInterface_& inplace_unary(Op&& op) {
        detail::inplace_tensor_op(op, *this);
        return *this;
      }

      // Scale operation


      template <typename Scalar,
          typename std::enable_if<detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor scale(const Scalar factor) const {
        return unary([=] (const numeric_type a) -> numeric_type
            { return a * factor; });
      }


      template <typename Scalar,
          typename std::enable_if<detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor scale(const Scalar factor, const Permutation& perm) const {
        return unary([=] (const numeric_type a) -> numeric_type
            { return a * factor; }, perm);
      }


      template <typename Scalar,
          typename std::enable_if<detail::is_numeric<Scalar>::value>::type* = nullptr>
      TensorInterface_& scale_to(const Scalar factor) {
        return inplace_unary([=] (numeric_type& MADNESS_RESTRICT res) { res *= factor; });
      }


      // Addition operations


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      result_tensor add(const Right& right) const {
        return binary(right, [] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return l + r; });
      }


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      result_tensor add(const Right& right, const Permutation& perm) const {
        return binary(right, [] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return l + r; }, perm);
      }


      template <typename Right, typename Scalar,
          typename std::enable_if<is_tensor<Right>::value &&
          detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor add(const Right& right, const Scalar factor) const {
        return binary(right, [=] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return (l + r) * factor; });
      }


      template <typename Right, typename Scalar,
          typename std::enable_if<is_tensor<Right>::value &&
          detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor add(const Right& right, const Scalar factor,
          const Permutation& perm) const
      {
        return binary(right,  [=] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return (l + r) * factor; }, perm);
      }


      result_tensor add(const numeric_type value) const {
        return unary([=] (const numeric_type a) -> numeric_type
            { return a + value; });
      }


      result_tensor add(const numeric_type value, const Permutation& perm) const {
        return unary([=] (const numeric_type a) -> numeric_type
            { return a + value; }, perm);
      }


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      TensorInterface_& add_to(const Right& right) {
        return inplace_binary(right, [] (numeric_type& MADNESS_RESTRICT l,
            const numeric_t<Right> r)
            { l += r; });
      }


      template <typename Right, typename Scalar,
          typename std::enable_if<is_tensor<Right>::value &&
          detail::is_numeric<Scalar>::value>::type* = nullptr>
      TensorInterface_& add_to(const Right& right, const Scalar factor) {
        return inplace_binary(right, [=] (numeric_type& MADNESS_RESTRICT l,
            const numeric_t<Right> r)
            { (l += r) *= factor; });
      }


      TensorInterface_& add_to(const numeric_type value) {
        return inplace_unary([=] (numeric_type& MADNESS_RESTRICT res) { res += value; });
      }

      // Subtraction operations


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      result_tensor subt(const Right& right) const {
        return binary(right, [] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return l - r; });
      }


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      result_tensor subt(const Right& right, const Permutation& perm) const {
        return binary(right, [] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return l - r; }, perm);
      }


      template <typename Right, typename Scalar,
          typename std::enable_if<is_tensor<Right>::value &&
          detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor subt(const Right& right, const numeric_type factor) const {
        return binary(right, [=] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return (l - r) * factor; });
      }


      template <typename Right, typename Scalar,
          typename std::enable_if<is_tensor<Right>::value &&
          detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor subt(const Right& right, const Scalar factor,
          const Permutation& perm) const
      {
        return binary(right, [=] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return (l - r) * factor; }, perm);
      }


      result_tensor subt(const numeric_type value) const {
        return add(-value);
      }


      result_tensor subt(const numeric_type value, const Permutation& perm) const {
        return add(-value, perm);
      }


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      TensorInterface_& subt_to(const Right& right) {
        return inplace_binary(right, [] (numeric_type& MADNESS_RESTRICT l,
            const numeric_t<Right> r)
            { l -= r; });
      }


      template <typename Right, typename Scalar,
          typename std::enable_if<is_tensor<Right>::value &&
          detail::is_numeric<Scalar>::value>::type* = nullptr>
      TensorInterface_& subt_to(const Right& right, const Scalar factor) {
        return inplace_binary(right, [=] (numeric_type& MADNESS_RESTRICT l,
            const numeric_t<Right> r)
            { (l -= r) *= factor; });
      }


      TensorInterface_& subt_to(const numeric_type value) {
        return add_to(-value);
      }

      // Multiplication operations


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      result_tensor mult(const Right& right) const {
        return binary(right, [] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return l * r; });
      }


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      result_tensor mult(const Right& right, const Permutation& perm) const {
        return binary(right, [] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return l * r; }, perm);
      }


      template <typename Right, typename Scalar,
          typename std::enable_if<is_tensor<Right>::value &&
          detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor mult(const Right& right, const Scalar factor) const {
        return binary(right, [=] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return (l * r) * factor; });
      }


      template <typename Right, typename Scalar,
          typename std::enable_if<is_tensor<Right>::value &&
          detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor mult(const Right& right, const Scalar factor,
          const Permutation& perm) const
      {
        return binary(right,  [=] (const numeric_type l,
            const numeric_t<Right> r)
            -> numeric_type { return (l * r) * factor; }, perm);
      }


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      TensorInterface_& mult_to(const Right& right) {
        return inplace_binary(right, [] (numeric_type& MADNESS_RESTRICT l,
            const numeric_t<Right> r)
            { l *= r; });
      }


      template <typename Right, typename Scalar,
          typename std::enable_if<is_tensor<Right>::value &&
          detail::is_numeric<Scalar>::value>::type* = nullptr>
      TensorInterface_& mult_to(const Right& right, const Scalar factor) {
        return inplace_binary(right, [=] (numeric_type& MADNESS_RESTRICT l,
            const numeric_t<Right> r)
            { (l *= r) *= factor; });
      }

      // Negation operations


      result_tensor neg() const {
        return unary([] (const numeric_type r) -> numeric_type { return -r; });
      }


      result_tensor neg(const Permutation& perm) const {
        return unary([] (const numeric_type l) -> numeric_type { return -l; },
            perm);
      }


      TensorInterface_& neg_to() {
        return inplace_unary([] (numeric_type& MADNESS_RESTRICT l) { l = -l; });
      }



      result_tensor conj() const { return scale(conj_op()); }


      template <typename Scalar,
          typename std::enable_if<detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor conj(const Scalar factor) const {
        return scale(conj_op(factor));
      }


      result_tensor conj(const Permutation& perm) const {
        return scale(conj_op(), perm);
      }


      template <typename Scalar,
          typename std::enable_if<detail::is_numeric<Scalar>::value>::type* = nullptr>
      result_tensor conj(const Scalar factor, const Permutation& perm) const {
        return scale(conj_op(factor), perm);
      }


      TensorInterface_& conj_to() {
        return scale_to(conj_op());
      }


      template <typename Scalar,
          typename std::enable_if<detail::is_numeric<Scalar>::value>::type* = nullptr>
      TensorInterface_& conj_to(const Scalar factor) {
        return scale_to(conj_op(factor));
      }



      template <typename ReduceOp, typename JoinOp>
      numeric_type reduce(ReduceOp&& reduce_op, JoinOp&& join_op,
          const numeric_type identity) const
      {
        return detail::tensor_reduce(reduce_op, join_op, identity, *this);
      }


      template <typename Right, typename ReduceOp, typename JoinOp,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      numeric_type reduce(const Right& other, ReduceOp&& reduce_op, JoinOp&& join_op,
          const numeric_type identity) const
      {
        return detail::tensor_reduce(reduce_op, join_op, identity, *this, other);
      }


      numeric_type sum() const {
        auto sum_op = [] (numeric_type& MADNESS_RESTRICT res, const numeric_type arg)
                { res += arg; };
        return reduce(sum_op, sum_op, numeric_type(0));
      }


      numeric_type product() const {
        auto mult_op = [] (numeric_type& MADNESS_RESTRICT res, const numeric_type arg)
                { res *= arg; };
        return reduce(mult_op, mult_op, numeric_type(1));
      }


      scalar_type squared_norm() const {
        auto square_op = [] (scalar_type& MADNESS_RESTRICT res, const numeric_type arg)
                { res += TiledArray::detail::norm(arg); };
        auto sum_op = [] (scalar_type& MADNESS_RESTRICT res, const scalar_type arg)
                { res += arg; };
        return reduce(square_op, sum_op, numeric_type(0));
      }


      numeric_type norm() const {
        return std::sqrt(squared_norm());
      }


      numeric_type min() const {
        auto min_op = [] (numeric_type& MADNESS_RESTRICT res, const numeric_type arg)
                { res = std::min(res, arg); };
        return reduce(min_op, min_op, std::numeric_limits<numeric_type>::max());
      }


      numeric_type max() const {
        auto max_op = [] (numeric_type& MADNESS_RESTRICT res, const numeric_type arg)
                { res = std::max(res, arg); };
        return reduce(max_op, max_op, std::numeric_limits<numeric_type>::min());
      }


      numeric_type abs_min() const {
        auto abs_min_op = [] (numeric_type& MADNESS_RESTRICT res, const numeric_type arg)
                { res = std::min(res, std::abs(arg)); };
        auto min_op = [] (numeric_type& MADNESS_RESTRICT res, const numeric_type arg)
                { res = std::min(res, arg); };
        return reduce(abs_min_op, min_op, std::numeric_limits<numeric_type>::max());
      }


      numeric_type abs_max() const {
        auto abs_max_op = [] (numeric_type& MADNESS_RESTRICT res, const numeric_type arg)
                { res = std::max(res, std::abs(arg)); };
        auto max_op = [] (numeric_type& MADNESS_RESTRICT res, const numeric_type arg)
                { res = std::max(res, arg); };
        return reduce(abs_max_op, max_op, numeric_type(0));
      }


      template <typename Right,
          typename std::enable_if<is_tensor<Right>::value>::type* = nullptr>
      numeric_type dot(const Right& other) const {
        auto mult_add_op = [] (numeric_type& res, const numeric_type l,
                  const numeric_t<Right> r)
                  { res += l * r; };
        auto add_op = [] (numeric_type& MADNESS_RESTRICT res, const numeric_type value)
              { res += value; };
        return reduce(other, mult_add_op, add_op, numeric_type(0));
      }

    }; // class TensorInterface

  } // namespace detail
} // namespace TiledArray

#endif // TILEDARRAY_TENSOR_TENSOR_VIEW_H__INCLUDED