Contract.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 Expression Templates for contracting tensor indices on a single
6 /// Tensor
7 
8 #pragma once
9 
12 #include "Utilities/Requires.hpp"
13 
14 /*!
15  * \ingroup TensorExpressionsGroup
16  * Holds all possible TensorExpressions currently implemented
17  */
18 namespace TensorExpressions {
19 
20 namespace detail {
21 
22 template <typename I1, typename I2>
23 using indices_contractible =
25  I1::dim == I2::dim and I1::ul != I2::ul and
26  I1::fr == I2::fr and I1::index == I2::index>;
27 
28 template <typename T, typename X, typename SymmList, typename IndexList,
29  typename Args>
30 struct ComputeContractedTypeImpl;
31 
32 template <typename T, typename X, template <typename...> class SymmList,
33  typename IndexList, typename Args, typename... Symm>
34 struct ComputeContractedTypeImpl<T, X, SymmList<Symm...>, IndexList, Args> {
35  using type =
36  TensorExpression<T, X, Symmetry<Symm::value...>, IndexList, Args>;
37 };
38 
39 template <typename Index1, typename Index2, typename T, typename X,
40  typename Symm, typename IndexList, typename Args>
41 using ComputeContractedType = typename ComputeContractedTypeImpl<
42  T, X, tmpl::erase<tmpl::erase<Symm, Index2>, Index1>,
43  tmpl::erase<tmpl::erase<IndexList, Index2>, Index1>,
44  tmpl::erase<tmpl::erase<Args, Index2>, Index1>>::type;
45 
46 template <int I, typename Index1, typename Index2>
47 struct ComputeContractionImpl {
48  template <typename... LhsIndices, typename T, typename S>
49  static SPECTRE_ALWAYS_INLINE typename T::type apply(S tensor_index,
50  const T& t1) {
51  tensor_index[Index1::value] = I;
52  tensor_index[Index2::value] = I;
53  return t1.template get<LhsIndices...>(tensor_index) +
54  ComputeContractionImpl<I - 1, Index1, Index2>::template apply<
55  LhsIndices...>(tensor_index, t1);
56  }
57 };
58 
59 template <typename Index1, typename Index2>
60 struct ComputeContractionImpl<0, Index1, Index2> {
61  template <typename... LhsIndices, typename T, typename S>
62  static SPECTRE_ALWAYS_INLINE typename T::type apply(S tensor_index,
63  const T& t1) {
64  tensor_index[Index1::value] = 0;
65  tensor_index[Index2::value] = 0;
66  return t1.template get<LhsIndices...>(tensor_index);
67  }
68 };
69 } // namespace detail
70 
71 /*!
72  * \ingroup TensorExpressionsGroup
73  */
74 template <typename Index1, typename Index2, typename T, typename X,
75  typename Symm, typename IndexList, typename ArgsList>
77  : public TensorExpression<
78  TensorContract<Index1, Index2, T, X, Symm, IndexList, ArgsList>, X,
79  typename detail::ComputeContractedType<Index1, Index2, T, X, Symm,
80  IndexList, ArgsList>::symmetry,
81  typename detail::ComputeContractedType<
82  Index1, Index2, T, X, Symm, IndexList, ArgsList>::index_list,
83  typename detail::ComputeContractedType<
84  Index1, Index2, T, X, Symm, IndexList, ArgsList>::args_list>,
85  public Expression {
86  using CI1 = tmpl::at<IndexList, Index1>;
87  using CI2 = tmpl::at<IndexList, Index2>;
88  static_assert(tmpl::size<Symm>::value > 1 and
89  tmpl::size<IndexList>::value > 1,
90  "Cannot contract indices on a Tensor with rank less than 2");
92  "Cannot contract the requested indices.");
93 
94  using new_type = detail::ComputeContractedType<Index1, Index2, T, X, Symm,
95  IndexList, ArgsList>;
96 
97  using type = X;
98  using symmetry = typename new_type::symmetry;
99  using index_list = typename new_type::index_list;
100  static constexpr auto num_tensor_indices =
101  tmpl::size<index_list>::value == 0 ? 1 : tmpl::size<index_list>::value;
102  using args_list = tmpl::sort<typename new_type::args_list>;
103 
104  explicit TensorContract(
105  const TensorExpression<T, X, Symm, IndexList, ArgsList>& t)
106  : t_(~t) {}
107 
108  template <size_t I, size_t Rank, Requires<(I <= Index1::value)> = nullptr>
109  SPECTRE_ALWAYS_INLINE void fill_contracting_tensor_index(
110  std::array<int, Rank>& tensor_index_in,
111  const std::array<int, num_tensor_indices>& tensor_index) const {
112  // -100 is for the slot that will be set later. Easy to debug.
113  tensor_index_in[I] = I == Index1::value ? -100 : tensor_index[I];
114  fill_contracting_tensor_index<I + 1>(tensor_index_in, tensor_index);
115  }
116 
117  template <size_t I, size_t Rank,
118  Requires<(I > Index1::value and I <= Index2::value and
119  I < Rank - 1)> = nullptr>
120  SPECTRE_ALWAYS_INLINE void fill_contracting_tensor_index(
121  std::array<int, Rank>& tensor_index_in,
122  const std::array<int, Rank - 2>& tensor_index) const {
123  // tensor_index is Rank - 2 since it shouldn't be called for Rank 2 case
124  // -200 is for the slot that will be set later. Easy to debug.
125  tensor_index_in[I] = I == Index2::value ? -200 : tensor_index[I - 1];
126  fill_contracting_tensor_index<I + 1>(tensor_index_in, tensor_index);
127  }
128 
129  template <size_t I, size_t Rank,
130  Requires<(I > Index2::value and I < Rank - 1)> = nullptr>
131  SPECTRE_ALWAYS_INLINE void fill_contracting_tensor_index(
132  std::array<int, Rank>& tensor_index_in,
133  const std::array<int, Rank - 2>& tensor_index) const {
134  // Left as Rank - 2 since it should never be called for the Rank 2 case
135  tensor_index_in[I] = tensor_index[I - 2];
136  fill_contracting_tensor_index<I + 1>(tensor_index_in, tensor_index);
137  }
138 
139  template <size_t I, size_t Rank, Requires<(I == Rank - 1)> = nullptr>
140  SPECTRE_ALWAYS_INLINE void fill_contracting_tensor_index(
141  std::array<int, Rank>& tensor_index_in,
142  const std::array<int, num_tensor_indices>& tensor_index) const {
143  tensor_index_in[I] = tensor_index[I - 2];
144  }
145 
146  template <typename... LhsIndices, typename U>
148  get(const std::array<U, num_tensor_indices>& new_tensor_index) const {
149  // new_tensor_index is the one with _fewer_ components, ie post-contraction
151  // Manually unrolled for loops to compute the tensor_index from the
152  // new_tensor_index
153  fill_contracting_tensor_index<0>(tensor_index, new_tensor_index);
154  return detail::ComputeContractionImpl<CI1::dim - 1, Index1, Index2>::
155  template apply<LhsIndices...>(tensor_index, t_);
156  }
157 
158  private:
160  TensorExpression<T, X, Symm, IndexList, ArgsList>>
161  t_;
162 };
163 
164 /*!
165  * \ingroup TensorExpressionsGroup
166  */
167 template <int Index1, int Index2, typename T, typename X, typename Symm,
168  typename IndexList, typename Args>
169 SPECTRE_ALWAYS_INLINE auto contract(
170  const TensorExpression<T, X, Symm, IndexList, Args>& t) {
171  return TensorContract<tmpl::int32_t<Index1>, tmpl::int32_t<Index2>, T, X,
172  Symm, IndexList, Args>(~t);
173 }
174 
175 namespace detail {
176 // Helper struct to allow contractions by using repeated indices in operator()
177 // calls to tensor.
178 template <template <typename> class TE, typename ReplacedArgList, typename I,
179  typename TotalContracted>
180 struct fully_contract_helper {
181  template <typename T>
182  SPECTRE_ALWAYS_INLINE static constexpr auto apply(const T& t) -> decltype(
183  contract<
184  tmpl::index_of<ReplacedArgList, ti_contracted_t<I::value>>::value,
185  tmpl::index_of<ReplacedArgList,
186  ti_contracted_t<I::value + 1>>::value>(
187  fully_contract_helper<TE, ReplacedArgList, tmpl::next<I>,
188  TotalContracted>::apply(t))) {
189  return contract<
190  tmpl::index_of<ReplacedArgList, ti_contracted_t<I::value>>::value,
191  tmpl::index_of<ReplacedArgList, ti_contracted_t<I::value + 1>>::value>(
192  fully_contract_helper<TE, ReplacedArgList, tmpl::next<I>,
193  TotalContracted>::apply(t));
194  }
195 };
196 
197 template <template <typename> class TE, typename ReplacedArgList,
198  typename TotalContracted>
199 struct fully_contract_helper<TE, ReplacedArgList,
200  tmpl::int32_t<TotalContracted::value - 1>,
201  TotalContracted> {
202  using I = tmpl::int32_t<TotalContracted::value - 1>;
203  template <typename T>
204  SPECTRE_ALWAYS_INLINE static constexpr auto apply(const T& t) -> decltype(
205  contract<
206  tmpl::index_of<ReplacedArgList, ti_contracted_t<I::value>>::value,
207  tmpl::index_of<ReplacedArgList,
208  ti_contracted_t<I::value + 1>>::value>(
209  TE<ReplacedArgList>(t))) {
210  return contract<
211  tmpl::index_of<ReplacedArgList, ti_contracted_t<I::value>>::value,
212  tmpl::index_of<ReplacedArgList, ti_contracted_t<I::value + 1>>::value>(
213  TE<ReplacedArgList>(t));
214  }
215 };
216 } // namespace detail
217 
218 /*!
219  * \ingroup TensorExpressionsGroup
220  * \brief Represents a fully contracted Tensor
221  */
222 template <template <typename> class TE, typename ReplacedArgList, typename I,
223  typename TotalContracted>
224 using fully_contracted =
225  detail::fully_contract_helper<TE, ReplacedArgList, I, TotalContracted>;
226 } // namespace TensorExpressions
Marks a class as being a TensorExpression.
Definition: TensorExpression.hpp:271
Defines the type alias Requires.
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1595
Definition: Determinant.hpp:11
typename detail::SymmetryImpl< std::make_index_sequence< sizeof...(T)>, tmpl::integral_list< std::int32_t, T... > >::type Symmetry
Computes the canonical symmetry from the integers T
Definition: Symmetry.hpp:81
#define SPECTRE_ALWAYS_INLINE
Always inline a function. Only use this if you benchmarked the code.
Definition: ForceInline.hpp:20
Defines class.
Defines metafunctions used to comute the Symmetry<...> for a Tensor.
Definition: AddSubtract.hpp:28
detail::fully_contract_helper< TE, ReplacedArgList, I, TotalContracted > fully_contracted
Represents a fully contracted Tensor.
Definition: Contract.hpp:225
typename Requires_detail::requires_impl< B >::template_error_type_failed_to_meet_requirements_on_template_parameters Requires
Express requirements on the template parameters of a function or class, replaces std::enable_if_t ...
Definition: Requires.hpp:67
Definition: Contract.hpp:76