TensorAsExpression.hpp
Go to the documentation of this file.
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 /// \file
5 /// Defines expressions that represent tensors
6 
7 #pragma once
8 
9 #include <array>
10 #include <cstddef>
11 
12 #include "DataStructures/Tensor/Expressions/LhsTensorSymmAndIndices.hpp"
17 #include "Utilities/Algorithm.hpp"
20 #include "Utilities/TMPL.hpp"
21 #include "Utilities/TypeTraits/IsA.hpp"
22 
23 namespace TensorExpressions {
24 /// \ingroup TensorExpressionsGroup
25 /// \brief Defines an expression representing a Tensor
26 ///
27 /// \details
28 /// In order to represent a tensor as an expression, instead of having Tensor
29 /// derive off of TensorExpression, a TensorAsExpression derives off of
30 /// TensorExpression and contains a pointer to a Tensor. The reason having
31 /// Tensor derive off of TensorExpression is problematic is that the index
32 /// structure is part of the type of the TensorExpression, so every possible
33 /// permutation and combination of indices must be derived from. For a rank 3
34 /// tensor, this is already over 500 base classes, which the Intel compiler
35 /// takes too long to compile.
36 ///
37 /// \tparam T the type of Tensor being represented as an expression
38 /// \tparam ArgsList the tensor indices, e.g. `_a` and `_b` in `F(_a, _b)`
39 template <typename T, typename ArgsList>
41 
42 template <typename X, typename Symm, template <typename...> class IndexList,
43  typename... Indices, template <typename...> class ArgsList,
44  typename... Args>
45 struct TensorAsExpression<Tensor<X, Symm, IndexList<Indices...>>,
46  ArgsList<Args...>>
47  : public TensorExpression<
48  TensorAsExpression<Tensor<X, Symm, IndexList<Indices...>>,
49  ArgsList<Args...>>,
50  X, Symm, IndexList<Indices...>, ArgsList<Args...>> {
51  using type = X;
52  using symmetry = Symm;
53  using index_list = IndexList<Indices...>;
54  static constexpr auto num_tensor_indices = tmpl::size<index_list>::value;
55  using args_list = ArgsList<Args...>;
56  using structure = Tensor_detail::Structure<symmetry, Indices...>;
57 
58  /// Construct an expression from a Tensor
59  explicit TensorAsExpression(const Tensor<X, Symm, IndexList<Indices...>>& t)
60  : t_(&t) {}
61  ~TensorAsExpression() override = default;
62 
63  /// \brief Computes the right hand side tensor multi-index that corresponds to
64  /// the left hand side tensor multi-index, according to their generic indices
65  ///
66  /// \details
67  /// Given the order of the generic indices for the left hand side (LHS) and
68  /// right hand side (RHS) and a specific LHS tensor multi-index, the
69  /// computation of the equivalent multi-index for the RHS tensor accounts for
70  /// differences in the ordering of the generic indices on the LHS and RHS.
71  ///
72  /// Here, the elements of `lhs_index_order` and `rhs_index_order` refer to
73  /// TensorIndex::values that correspond to generic tensor indices,
74  /// `lhs_tensor_multi_index` is a multi-index for the LHS tensor, and the
75  /// equivalent RHS tensor multi-index is returned. If we have LHS tensor
76  /// \f$L_{ab}\f$, RHS tensor \f$R_{ba}\f$, and the LHS component \f$L_{31}\f$,
77  /// the corresponding RHS component is \f$R_{13}\f$.
78  ///
79  /// Here is an example of what the algorithm does:
80  ///
81  /// `lhs_index_order`:
82  /// \code
83  /// [0, 1, 2] // i.e. abc
84  /// \endcode
85  /// `rhs_index_order`:
86  /// \code
87  /// [1, 2, 0] // i.e. bca
88  /// \endcode
89  /// `lhs_tensor_multi_index`:
90  /// \code
91  /// [4, 0, 3] // i.e. a = 4, b = 0, c = 3
92  /// \endcode
93  /// returned RHS tensor multi-index:
94  /// \code
95  /// [0, 3, 4] // i.e. b = 0, c = 3, a = 4
96  /// \endcode
97  ///
98  /// \param lhs_index_order the generic index order of the LHS tensor
99  /// \param rhs_index_order the generic index order of the RHS tensor
100  /// \param lhs_tensor_multi_index the specific LHS tensor multi-index
101  /// \return the RHS tensor multi-index that corresponds to
102  /// `lhs_tensor_multi_index`, according to the index orders in
103  /// `lhs_index_order` and `rhs_index_order`
104  SPECTRE_ALWAYS_INLINE static constexpr std::array<size_t, sizeof...(Indices)>
106  const std::array<size_t, sizeof...(Indices)>& lhs_index_order,
107  const std::array<size_t, sizeof...(Indices)>& rhs_index_order,
108  const std::array<size_t, sizeof...(Indices)>&
109  lhs_tensor_multi_index) noexcept {
110  std::array<size_t, sizeof...(Indices)> rhs_tensor_multi_index{};
111  for (size_t i = 0; i < sizeof...(Indices); ++i) {
112  gsl::at(rhs_tensor_multi_index,
113  static_cast<unsigned long>(std::distance(
114  rhs_index_order.begin(),
115  alg::find(rhs_index_order, gsl::at(lhs_index_order, i))))) =
116  gsl::at(lhs_tensor_multi_index, i);
117  }
118  return rhs_tensor_multi_index;
119  }
120 
121  /// \brief Computes a mapping from the storage indices of the left hand side
122  /// tensor to the right hand side tensor
123  ///
124  /// \tparam LhsStructure the Structure of the Tensor on the left hand side of
125  /// the TensorExpression
126  /// \tparam LhsIndices the TensorIndexs of the Tensor on the left hand side
127  /// \return the mapping from the left hand side to the right hand side storage
128  /// indices
129  template <typename LhsStructure, typename... LhsIndices>
130  SPECTRE_ALWAYS_INLINE static constexpr std::array<size_t,
131  LhsStructure::size()>
133  constexpr size_t num_components = LhsStructure::size();
134  std::array<size_t, num_components> lhs_to_rhs_map{};
135  const auto lhs_storage_to_tensor_indices =
136  LhsStructure::storage_to_tensor_index();
137  for (size_t lhs_storage_index = 0; lhs_storage_index < num_components;
138  ++lhs_storage_index) {
139  // `compute_rhs_tensor_index` will return the RHS tensor multi-index that
140  // corresponds to the LHS tensor multi-index, according to the order of
141  // the generic indices for the LHS and RHS. structure::get_storage_index
142  // will then get the RHS storage index that corresponds to this RHS
143  // tensor multi-index.
144  gsl::at(lhs_to_rhs_map, lhs_storage_index) =
145  structure::get_storage_index(compute_rhs_tensor_index(
146  {{LhsIndices::value...}}, {{Args::value...}},
147  lhs_storage_to_tensor_indices[lhs_storage_index]));
148  }
149  return lhs_to_rhs_map;
150  }
151 
152  /// \brief Returns the value at a left hand side tensor's storage index
153  ///
154  /// \details
155  /// One big challenge with TensorExpression implementation is the reordering
156  /// of the indices on the left hand side (LHS) and right hand side (RHS) of
157  /// the expression. The algorithms implemented in `compute_lhs_to_rhs_map` and
158  /// `compute_rhs_tensor_index` handle the index sorting by mapping between the
159  /// generic index orders of the LHS and RHS.
160  ///
161  /// \tparam LhsStructure the Structure of the Tensor on the LHS of the
162  /// TensorExpression
163  /// \tparam LhsIndices the TensorIndexs of the Tensor on the LHS of the tensor
164  /// expression
165  /// \param lhs_storage_index the storage index of the LHS tensor component to
166  /// retrieve
167  /// \return the value of the DataType of the component at `lhs_storage_index`
168  /// in the LHS tensor
169  template <typename LhsStructure, typename... LhsIndices>
170  SPECTRE_ALWAYS_INLINE decltype(auto) get(
171  const size_t lhs_storage_index) const noexcept {
172  if constexpr (std::is_same_v<LhsStructure, structure> and
173  std::is_same_v<tmpl::list<LhsIndices...>,
174  tmpl::list<Args...>>) {
175  return (*t_)[lhs_storage_index];
176  } else {
177  constexpr std::array<size_t, LhsStructure::size()> lhs_to_rhs_map =
178  compute_lhs_to_rhs_map<LhsStructure, LhsIndices...>();
179  return (*t_)[gsl::at(lhs_to_rhs_map, lhs_storage_index)];
180  }
181  }
182 
183  /// Retrieve the i'th entry of the Tensor being held
184  SPECTRE_ALWAYS_INLINE type operator[](const size_t i) const {
185  return t_->operator[](i);
186  }
187 
188  private:
189  const Tensor<X, Symm, IndexList<Indices...>>* t_ = nullptr;
190 };
191 } // namespace TensorExpressions
gsl::at
constexpr T & at(std::array< T, N > &arr, Size index)
Retrieve a entry from a container, with checks in Debug mode that the index being retrieved is valid.
Definition: Gsl.hpp:125
TensorExpressions::TensorAsExpression< Tensor< X, Symm, IndexList< Indices... > >, ArgsList< Args... > >::operator[]
type operator[](const size_t i) const
Retrieve the i'th entry of the Tensor being held.
Definition: TensorAsExpression.hpp:184
TensorExpression.hpp
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:660
Symmetry.hpp
Structure.hpp
TensorExpressions::TensorAsExpression< Tensor< X, Symm, IndexList< Indices... > >, ArgsList< Args... > >::TensorAsExpression
TensorAsExpression(const Tensor< X, Symm, IndexList< Indices... >> &t)
Construct an expression from a Tensor.
Definition: TensorAsExpression.hpp:59
TensorExpressions::TensorAsExpression
Defines an expression representing a Tensor.
Definition: TensorAsExpression.hpp:40
TensorExpressions
Definition: AddSubtract.hpp:30
SPECTRE_ALWAYS_INLINE
#define SPECTRE_ALWAYS_INLINE
Definition: ForceInline.hpp:16
TensorExpressions::TensorAsExpression< Tensor< X, Symm, IndexList< Indices... > >, ArgsList< Args... > >::compute_rhs_tensor_index
static constexpr std::array< size_t, sizeof...(Indices)> compute_rhs_tensor_index(const std::array< size_t, sizeof...(Indices)> &lhs_index_order, const std::array< size_t, sizeof...(Indices)> &rhs_index_order, const std::array< size_t, sizeof...(Indices)> &lhs_tensor_multi_index) noexcept
Computes the right hand side tensor multi-index that corresponds to the left hand side tensor multi-i...
Definition: TensorAsExpression.hpp:105
cstddef
Assert.hpp
array
ForceInline.hpp
alg::find
constexpr decltype(auto) find(Container &&c, const T &value)
Convenience wrapper around constexpr reimplementation of std::find.
Definition: Algorithm.hpp:225
Tensor.hpp
TensorExpressions::TensorAsExpression< Tensor< X, Symm, IndexList< Indices... > >, ArgsList< Args... > >::compute_lhs_to_rhs_map
static constexpr std::array< size_t, LhsStructure::size()> compute_lhs_to_rhs_map() noexcept
Computes a mapping from the storage indices of the left hand side tensor to the right hand side tenso...
Definition: TensorAsExpression.hpp:132
TMPL.hpp