ALPSCore reference
tensor_base.hpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1998-2018 ALPS Collaboration. See COPYRIGHT.TXT
3  * All rights reserved. Use is subject to license terms. See LICENSE.TXT
4  * For use in publications, see ACKNOWLEDGE.TXT
5  */
6 
7 #ifndef ALPSCORE_GF_TENSOR_H
8 #define ALPSCORE_GF_TENSOR_H
9 
10 
11 #include <array>
12 #include <iostream>
13 #include <numeric>
14 #include <type_traits>
15 #include <vector>
16 
17 #include <Eigen/Core>
18 #include <Eigen/Dense>
19 
23 
24 
25 namespace alps {
26  namespace numerics {
27  namespace detail {
33  template<typename T, typename St>
34  struct is_storage {
35  static constexpr bool value = std::is_same < St, simple_storage < T > > ::value ||
36  std::is_same < St, data_view < T > > ::value;
37  };
38 
44  template <bool...> struct bool_pack;
45  template <bool... v>
46  using all_true = std::is_same<bool_pack<true, v...>, bool_pack<v..., true>>;
47 
55  template<typename T, size_t D, typename C>
56  class tensor_base;
57  }
61  template<typename T, size_t D>
62  using tensor = detail::tensor_base < T, D, detail::data_storage < T > >;
63 #ifdef ALPS_HAVE_SHARED_ALLOCATOR
64 
67  template<typename T, size_t D>
68  using shared_tensor = detail::tensor_base < T, D, shared_storage< T > >;
69 #endif
70 
73  template<typename T, size_t D>
74  using tensor_view = detail::tensor_base < T, D, detail::data_view < T > >;
75 
76  namespace detail {
77  template<typename X, int Rows = Eigen::Dynamic, int Cols = Eigen::Dynamic>
78  using MatrixMap = Eigen::Map < Eigen::Matrix < X, Rows, Cols, Eigen::RowMajor > >;
79  template<typename X, int Rows = Eigen::Dynamic, int Cols = Eigen::Dynamic>
80  using ConstMatrixMap = Eigen::Map <const Eigen::Matrix < X, Rows, Cols, Eigen::RowMajor > >;
81 
87  template<typename T, size_t Dim, typename Container>
88  class tensor_base {
89  public:
90  // types definitions
91  typedef T prec;
92  private:
93  typedef data_view < T > viewType;
94  typedef data_view < const typename std::remove_const<T>::type > constViewType;
95  typedef data_view < typename std::remove_const<T>::type > nonconstViewType;
96  typedef data_storage < T > storageType;
97  typedef data_storage < const typename std::remove_const<T>::type > constStorageType;
98  typedef data_storage < typename std::remove_const<T>::type > nonconstStorageType;
100  typedef tensor_base < T, Dim, Container > tType;
102  typedef tensor < T, Dim > tensorType;
104  typedef tensor_view < T, Dim > tensorViewType;
106  template<typename St>
107  using genericTensor = tensor_base < T, Dim, St >;
108 
109  private:
110  // fields definitions
112  static size_t constexpr dim = Dim;
114  Container storage_;
116  std::array < size_t, Dim > shape_;
118  std::array < size_t, Dim > acc_sizes_;
119 
120  public:
121 
125  tensor_base() : storage_(0) {};
126 
131  tensor_base(Container &&container, const std::array < size_t, Dim >& sizes) : storage_(container), shape_(sizes) {
132  static_assert(is_storage< T, Container>::value, "Should be either data_storage or data_view type");
133  assert(storage_.size() == size());
134  fill_acc_sizes();
135  }
136 
137  tensor_base(Container &container, const std::array < size_t, Dim >& sizes) : storage_(container), shape_(sizes) {
138  static_assert(is_storage< T, Container>::value, "Should be either data_storage or data_view type");
139  assert(storage_.size() == size());
140  fill_acc_sizes();
141  }
142 
143  template<typename Container2>
144  tensor_base(Container2 &container, size_t size, size_t offset, const std::array < size_t, Dim >& sizes) :
145  storage_(container, size, offset), shape_(sizes) {
146  assert(storage_.size() == this->size());
147  fill_acc_sizes();
148  }
149 
150  template<typename Container2>
151  tensor_base(const Container2 &container, size_t size, size_t offset, const std::array < size_t, Dim >& sizes) :
152  storage_(container, size, offset), shape_(sizes) {
153  assert(storage_.size() == this->size());
154  fill_acc_sizes();
155  }
156 
157  template<typename...Indices>
158  tensor_base(typename std::enable_if< all_true<std::is_convertible<Indices, std::size_t>::value...>::value, Container >::type &container,
159  size_t size1, Indices...sizes) : tensor_base(container, {{size1, size_t(sizes)...}}) {}
160 
168  tensor_base(T *data, const std::array < size_t, Dim > & sizes) : storage_(data, size(sizes)), shape_(sizes) {
169  fill_acc_sizes();
170  }
171 
172  template<typename...Indices>
173  tensor_base(T *data, size_t size1, Indices...sizes) : tensor_base(data, {{size1, size_t(sizes)...}}) {}
174 
181  template<typename X = Container>
182  tensor_base(typename std::enable_if < std::is_same < X, data_storage < T > >::value,
183  const std::array < size_t, Dim > & >::type sizes) : storage_(size(sizes)), shape_(sizes) {
184  fill_acc_sizes();
185  }
186 
187  template<typename X = Container, typename...Indices>
188  tensor_base(typename std::enable_if < std::is_same < X, data_storage < T > >::value,
189  size_t>::type size1, Indices...sizes) : storage_(size({{size1, size_t(sizes)...}})), shape_({{size1, size_t(sizes)...}}) {
190  static_assert(sizeof...(Indices) + 1 == Dim, "Wrong dimension");
191  fill_acc_sizes();
192  }
193 
194  // this constructor create a view of other tensor. that is why rhs is not const
195  template<typename St = Container>
196  tensor_base(typename std::enable_if<std::is_same<St, viewType>::value, tensorType>::type& rhs) :
197  storage_(rhs.storage()), shape_(rhs.shape()), acc_sizes_(rhs.acc_sizes()) {}
198 
200  tensor_base(const tType& rhs) = default;
202  tensor_base(tType &&rhs) = default;
204  template<typename T2, typename St, typename = std::enable_if<std::is_same<Container, storageType>::value, void >>
205  tensor_base(const tensor_base<T2, Dim, St> &rhs) : storage_(rhs.storage()), shape_(rhs.shape()), acc_sizes_(rhs.acc_sizes()) {}
207  template<typename T2, typename St, typename = std::enable_if<std::is_same<Container, storageType>::value, void >>
208  tensor_base(tensor_base<T2, Dim, St> &&rhs) noexcept: storage_(rhs.storage()), shape_(rhs.shape()), acc_sizes_(rhs.acc_sizes()) {}
209 
211  template<typename T2, typename St>
212  tensor_base < T, Dim, Container > &operator=(const tensor_base < T2, Dim, St> &rhs){
213  assert(size()==rhs.size());
214  storage_ = rhs.storage();
215  return *this;
216  };
218  tensor_base < T, Dim, Container > &operator=(const tensor_base < T, Dim, Container > &rhs) = default;
220  tensor_base < T, Dim, Container > &operator=(tensor_base < T, Dim, Container > &&rhs) = default;
222  template<typename T2, typename St>
223  bool operator==(const tensor_base<T2, Dim, St>& rhs) const {
224  return std::equal(shape_.begin(), shape_.end(), rhs.shape().begin()) && storage_ == rhs.storage();
225  }
226 
234  template<typename ...IndexTypes>
235  const T & operator()(typename std::enable_if < sizeof...(IndexTypes) == Dim - 1, size_t >::type t1, IndexTypes ... indices) const {
236  return storage_.data(index(t1, indices...));
237  }
238 
242  template<typename ...IndexTypes>
243  T &operator()(typename std::enable_if < sizeof...(IndexTypes) == Dim - 1, size_t >::type t1, IndexTypes ... indices) {
244  return storage_.data(index(t1, indices...));
245  }
246 
247 
258  template<typename ...IndexTypes>
259  tensor_view < T, Dim - (sizeof...(IndexTypes)) - 1 > operator()(typename std::enable_if < (sizeof...(IndexTypes) < Dim - 1),
260  size_t >::type t1, IndexTypes ... indices) {
261  std::array < size_t, Dim - (sizeof...(IndexTypes)) - 1 > sizes;
262  size_t s = new_size(sizes);
263  return tensor_view < T, Dim - (sizeof...(IndexTypes)) - 1 > ( storage_, s, index(t1, indices...), sizes);
264  }
265 
266  template<typename ...IndexTypes>
267  tensor_view <const typename std::remove_const<T>::type, Dim - (sizeof...(IndexTypes)) - 1 > operator()
268  (typename std::enable_if < (sizeof...(IndexTypes) < Dim - 1), size_t >::type t1, IndexTypes ... indices) const {
269  std::array < size_t, Dim - (sizeof...(IndexTypes)) - 1 > sizes;
270  size_t s = new_size(sizes);
271  return tensor_view <const typename std::remove_const<T>::type, Dim - (sizeof...(IndexTypes)) - 1 >
272  (storage_, s, index(t1, indices...), sizes );
273  }
274 
275  /*
276  * Basic arithmetic operations
277  */
285  template<typename S>
286  typename std::enable_if < !std::is_same < S, tensorType >::value, tensor_base< decltype(S{} + T{}), Dim,
287  data_storage< decltype(S{} * T{}) > > >::type operator*(S scalar) const {
288  tensor_base< decltype(S{} + T{}), Dim, data_storage< decltype(S{} * T{}) > > x(*this);
289  return (x *= static_cast<decltype(S{} + T{})>(scalar));
290  };
291 
299  template<typename S>
300  typename std::enable_if < std::is_same < S, tensorType >::value, tensorType >::type operator*(const S& rhs) const {
301  tensorType x(*this);
302  return x*=rhs;
303  };
304 
308  template<typename S>
309  typename std::enable_if < !std::is_same < S, tensorType >::value, tType & >::type operator*=(S scalar) {
310  static_assert(std::is_convertible<S, T>::value, "Can't perform inplace multiplication: S can be casted into T");
311  Eigen::Map < Eigen::Matrix < T, 1, Eigen::Dynamic > > M(&storage_.data(0), storage_.size());
312  M *= T(scalar);
313  return *this;
314  };
315 
319  template<typename S>
320  typename std::enable_if < std::is_same < S, tensorType >::value, tensorType & >::type operator*=(const S& rhs) {
321  Eigen::Map < Eigen::Array < T, 1, Eigen::Dynamic > > M1(&storage_.data(0), storage_.size());
322  Eigen::Map < const Eigen::Array < T, 1, Eigen::Dynamic > > M2(&rhs.storage().data(0), rhs.storage().size());
323  M1*=M2;
324  return *this;
325  };
326 
330  template<typename S>
331  typename std::enable_if < !std::is_same < S, tensorType >::value, tensor_base< decltype(S{} + T{}), Dim,
332  data_storage< decltype(S{} * T{}) > > >::type operator/(S scalar) const {
333  tensor_base< decltype(S{} + T{}), Dim, data_storage< decltype(S{} * T{}) > > x(*this);
334  return (x /= scalar);
335  };
336 
340  template<typename S>
341  typename std::enable_if < !std::is_same < S, tensorType >::value, tType & >::type operator/=(S scalar) {
342  static_assert(std::is_convertible<S, T>::value, "Can not perform inplace division: S can be casted into T");
343  Eigen::Map < Eigen::Matrix < T, 1, Eigen::Dynamic > > M(&storage_.data(0), storage_.size());
344  M *= T(1.0)/T(scalar);
345  return *this;
346  };
347 
351  tensor_base < T, Dim, Container > operator-() const {
352  return *this * T(-1.0);
353  };
354 
358  void set_zero() {
359  set_number(T(0));
360  }
361 
368  template<typename num_type>
369  void set_number(num_type value) {
370  static_assert(std::is_convertible<num_type, T>::value, "Can not assign value to the tensor. Value can not be cast into the tensor value type.");
371  std::fill_n(data(), size(), T(value));
372  }
373 
377  size_t size() const {
378  return storage_.size();
379  }
380 
384  size_t num_elements() const {
385  return storage_.size();
386  }
387 
392  template<size_t M = Dim>
393  typename std::enable_if < M == 2, tensorType >::type
394  inverse() {
395  if (shape_[0] != shape_[1]) {
396  throw std::invalid_argument("Can not do inversion of the non-square matrix.");
397  }
398  tensorType x(*this);
399  Eigen::Map < Eigen::Matrix < T, Eigen::Dynamic, Eigen::Dynamic > > Mt(&(x.storage().data(0)), shape_[0], shape_[1]);
400  Mt = Mt.inverse();
401  return x;
402  };
403 
412  template<typename S, typename Ct>
413  tensor < decltype(S{} + T{}), Dim > operator+(const tensor_base < S, Dim, Ct > &rhs) {
414  tensor < decltype(S{} + T{}), Dim > x(*this);
415  return (x += rhs);
416  };
417 
421  template<typename S, typename Ct>
422  tensor < decltype(S{} - T{}), Dim > operator-(const tensor_base < S, Dim, Ct > &rhs) {
423  tensor < decltype(S{} - T{}), Dim > x(*this);
424  return (x -= rhs);
425  };
426 
430  template<typename S, typename Ct>
431  typename std::enable_if <
432  std::is_same < S, T >::value || std::is_same < T, std::complex < double>>::value
433  || std::is_same < T, std::complex < float>>::value, tType & >::type operator+=(const tensor_base < S, Dim, Ct > &y) {
434  MatrixMap < T, 1, Eigen::Dynamic > M1(&storage_.data(0), storage_.size());
435  ConstMatrixMap < S, 1, Eigen::Dynamic > M2(&y.storage().data(0), y.storage().size());
436  M1.noalias() += M2;
437  return (*this);
438  };
439 
443  template<typename S>
444  typename std::enable_if < std::is_same < S, tensorType >::value ||
445  std::is_same < S, tensorViewType >::value, tType & >::type operator-=(const S &y) {
446  MatrixMap < T, 1, Eigen::Dynamic > M1(&storage_.data(0), storage_.size());
447  ConstMatrixMap < T, 1, Eigen::Dynamic> M2(&y.storage().data(0), y.storage().size());
448  M1.noalias() -= M2;
449  return (*this);
450  };
451 
455  tensorType dot(const tensorType &y) const {
456  static_assert(Dim == 2, "Can not do inversion for not 2D tensor.");
457  if (shape_[1] != y.shape_[0]) {
458  throw std::invalid_argument("Can not do multiplication. Dimensions missmatches.");
459  }
460  tensorType x({{shape_[0], y.shape_[1]}});
461  ConstMatrixMap < T, Eigen::Dynamic, Eigen::Dynamic > M1(&storage().data(0), shape_[0], shape_[1]);
462  ConstMatrixMap < T, Eigen::Dynamic, Eigen::Dynamic > M2(&y.storage().data(0), y.shape()[0], y.shape()[1]);
463  MatrixMap < T > M3(&x.storage().data(0), x.shape()[0], x.shape()[1]);
464  M3 = M1*M2;
465  return x;
466  };
467 
471  MatrixMap < T > matrix() {
472  static_assert(Dim == 2, "Can not return Eigen matrix view for not 2D tensor.");
473  return MatrixMap < T >(&storage().data(0), shape_[0], shape_[1]);
474  };
475 
477  const std::array < size_t, Dim > &shape() const { return shape_; };
478 
480  template<typename ...Inds>
481  typename std::enable_if<are_all_integrals<Inds...>::value>::type reshape(Inds...inds) {
482  static_assert(sizeof...(Inds) == Dim, "New shape should have the same dimension.");
483  std::array<size_t, Dim> shape = {{size_t(inds)...}};
484  reshape(shape);
485  }
486 
488  template<typename X = Container>
489  typename std::enable_if<std::is_same < X, data_storage < T > >::value, void>::type reshape(const std::array<size_t, Dim>& shape) {
490  size_t new_size = size(shape);
491  storage_.resize(new_size);
492  shape_ = shape;
493  fill_acc_sizes();
494  }
495 
497  template<typename X = Container>
498  typename std::enable_if<std::is_same < X, data_view < T > >::value, void>::type reshape(const std::array<size_t, Dim>& shape) {
499  size_t new_size = size(shape);
500  if(new_size != size()) {
501  throw std::invalid_argument("Wrong size. Can't reshape tensor.");
502  }
503  shape_ = shape;
504  fill_acc_sizes();
505  }
506 
508  const std::array < size_t, Dim > &acc_sizes() const { return acc_sizes_; };
509 
511  static size_t dimension() { return Dim; }
512 
514  const prec *data() const { return storage_.data(); }
515 
517  prec *data() { return storage_.data(); }
518 
520  const Container &storage() const { return storage_; }
521 
523  Container &storage() { return storage_; }
524 
526  template<typename ...Indices>
527  inline size_t index(const Indices&...indices) const {
528  return index_impl(make_index_sequence<sizeof...(Indices)>(), size_t(indices)...);
529  }
530 
531  private:
535  template<size_t... I, typename ...Indices>
536  inline size_t index_impl(index_sequence<I...>, const Indices&... indices) const {
537  std::array<size_t, sizeof...(Indices)> a{{indices * acc_sizes_[I] ...}};
538  return std::accumulate(a.begin(), a.end(), size_t(0), std::plus<size_t>());
539  }
543  void fill_acc_sizes() {
544  acc_sizes_[shape_.size()-1] = 1;
545  for(int k = int(shape_.size()) - 2; k >= 0; --k)
546  acc_sizes_[k] = acc_sizes_[k+1] * shape_[k+1];
547  }
549  size_t size(const std::array < size_t, Dim > &shape) {
550  return std::accumulate(shape.begin(), shape.end(), size_t(1), std::multiplies<size_t>());
551  }
553  template<size_t M>
554  inline size_t new_size(std::array < size_t, M >& sizes) const {
555  size_t s = 1;
556  for (size_t i = 0; i < sizes.size(); ++i) {
557  sizes[i] = shape_[i + Dim - M];
558  s *= sizes[i];
559  }
560  return s;
561  }
562  };
563 
564  }
565  }
566 }
567 
568 #endif //ALPSCORE_GF_TENSOR_H
index_mesh::index_type index
Definition: mesh.hpp:1247
make_integer_sequence< size_t, N > make_index_sequence
detail::tensor_base< T, D, detail::data_storage< T > > tensor
Definition: tensor_base.hpp:62
boost::array< T, N > operator*(boost::array< T, N > lhs, boost::array< U, N > const &rhs)
std::enable_if<!is_sequence< T >::value, std::size_t >::type size(T const &)
Definition: size.hpp:20
bool operator==(const dictionary &lhs, const dictionary &rhs)
Definition: dictionary.hpp:96
boost::array< T, N > & operator+=(boost::array< T, N > &lhs, boost::array< T, N > const &rhs)
boost::array< T, N > & operator/=(boost::array< T, N > &lhs, boost::array< T, N > const &rhs)
boost::array< T, N > & operator*=(boost::array< T, N > &lhs, boost::array< T, N > const &rhs)
detail::tensor_base< T, D, detail::data_view< T > > tensor_view
Definition: tensor_base.hpp:74
boost::array< T, N > & operator-=(boost::array< T, N > &lhs, boost::array< T, N > const &rhs)
boost::array< T, N > operator-(boost::array< T, N > lhs)