7 #ifndef ALPSCORE_GF_TENSOR_H 8 #define ALPSCORE_GF_TENSOR_H 14 #include <type_traits> 18 #include <Eigen/Dense> 33 template<
typename T,
typename St>
35 static constexpr
bool value = std::is_same < St, simple_storage < T > > ::value ||
36 std::is_same < St, data_view < T > > ::value;
44 template <
bool...>
struct bool_pack;
46 using all_true = std::is_same<bool_pack<
true, v...>, bool_pack<v...,
true>>;
55 template<
typename T,
size_t D,
typename C>
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 67 template<
typename T,
size_t D>
68 using shared_tensor = detail::tensor_base < T, D, shared_storage< T > >;
73 template<
typename T,
size_t D>
74 using tensor_view = detail::tensor_base < T, D, detail::data_view < T > >;
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 > >;
87 template<
typename T,
size_t Dim,
typename Container>
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;
106 template<
typename St>
107 using genericTensor = tensor_base < T, Dim, St >;
112 static size_t constexpr dim = Dim;
116 std::array < size_t, Dim > shape_;
118 std::array < size_t, Dim > acc_sizes_;
125 tensor_base() : storage_(0) {};
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());
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());
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());
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());
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)...}}) {}
168 tensor_base(T *data,
const std::array < size_t, Dim > & sizes) : storage_(data,
size(sizes)), shape_(sizes) {
172 template<
typename...Indices>
173 tensor_base(T *data,
size_t size1, Indices...sizes) : tensor_base(data, {{size1, size_t(sizes)...}}) {}
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) {
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");
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()) {}
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()) {}
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();
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();
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...));
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...));
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);
266 template<
typename ...IndexTypes>
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);
272 (storage_, s,
index(t1, indices...), sizes );
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));
300 typename std::enable_if < std::is_same < S, tensorType >::value, tensorType >::type
operator*(
const S& rhs)
const {
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());
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());
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);
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);
351 tensor_base < T, Dim, Container >
operator-()
const {
352 return *
this * T(-1.0);
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));
377 size_t size()
const {
378 return storage_.size();
384 size_t num_elements()
const {
385 return storage_.size();
392 template<
size_t M = Dim>
393 typename std::enable_if < M == 2, tensorType >::type
395 if (shape_[0] != shape_[1]) {
396 throw std::invalid_argument(
"Can not do inversion of the non-square matrix.");
399 Eigen::Map < Eigen::Matrix < T, Eigen::Dynamic, Eigen::Dynamic > > Mt(&(x.storage().data(0)), shape_[0], shape_[1]);
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);
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);
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());
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());
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.");
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]);
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]);
477 const std::array < size_t, Dim > &shape()
const {
return shape_; };
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)...}};
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);
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.");
508 const std::array < size_t, Dim > &acc_sizes()
const {
return acc_sizes_; };
511 static size_t dimension() {
return Dim; }
514 const prec *data()
const {
return storage_.data(); }
517 prec *data() {
return storage_.data(); }
520 const Container &storage()
const {
return storage_; }
523 Container &storage() {
return storage_; }
526 template<
typename ...Indices>
527 inline size_t index(
const Indices&...indices)
const {
535 template<
size_t... I,
typename ...Indices>
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>());
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];
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>());
554 inline size_t new_size(std::array < size_t, M >& sizes)
const {
556 for (
size_t i = 0; i < sizes.size(); ++i) {
557 sizes[i] = shape_[i + Dim - M];
568 #endif //ALPSCORE_GF_TENSOR_H index_mesh::index_type index
make_integer_sequence< size_t, N > make_index_sequence
detail::tensor_base< T, D, detail::data_storage< T > > tensor
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 &)
bool operator==(const dictionary &lhs, const dictionary &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)
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
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)