Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <tuple>
9 : #include <utility>
10 :
11 : #include "DataStructures/Tensor/Tensor.hpp"
12 : #include "Utilities/MakeWithValue.hpp"
13 : #include "Utilities/TMPL.hpp"
14 :
15 : template <typename X, typename Symm, typename IndexList>
16 0 : class Tensor;
17 :
18 : namespace detail {
19 : // Get the values that encode the generic tensor indices for the first operand,
20 : // second operand, and result tensor so that the first NumIndicesToContract
21 : // indices will contract and the TensorExpression written with them will be
22 : // valid
23 : //
24 : // Note: Implementation assumes that the indices in the index pairs to contract
25 : // have the same type (spatial or spacetime)
26 : template <size_t NumIndicesToContract, size_t NumIndices1, size_t NumIndices2>
27 : constexpr auto
28 : get_tensor_index_values_for_tensors_to_contract_and_result_tensor(
29 : const std::array<bool, NumIndices1>& index_type_is_spacetime1,
30 : const std::array<bool, NumIndices1>& valence_is_lower1,
31 : const std::array<bool, NumIndices2>& index_type_is_spacetime2,
32 : const std::array<bool, NumIndices2>& valence_is_lower2) {
33 : constexpr size_t num_result_indices =
34 : NumIndices1 + NumIndices2 - 2 * NumIndicesToContract;
35 :
36 : // (first op index values, second op index values, result tensor index values)
37 : std::tuple<std::array<size_t, NumIndices1>, std::array<size_t, NumIndices2>,
38 : std::array<size_t, num_result_indices>>
39 : tensor_index_values{};
40 :
41 : // the next lower spacetime index value that we have not yet used
42 : size_t next_lower_spacetime_value = 0;
43 : // the next lower spatial index value that we have not yet used
44 : size_t next_lower_spatial_value = tenex::TensorIndex_detail::spatial_sentinel;
45 :
46 : // assign first operand's tensor index values
47 : for (size_t i = 0; i < NumIndices1; i++) {
48 : const bool index1_is_spacetime = gsl::at(index_type_is_spacetime1, i);
49 : const bool index1_is_lower = gsl::at(valence_is_lower1, i);
50 : if (index1_is_spacetime) {
51 : gsl::at(std::get<0>(tensor_index_values), i) =
52 : index1_is_lower ? next_lower_spacetime_value
53 : : tenex::get_tensorindex_value_with_opposite_valence(
54 : next_lower_spacetime_value);
55 : next_lower_spacetime_value++;
56 : } else {
57 : gsl::at(std::get<0>(tensor_index_values), i) =
58 : index1_is_lower ? next_lower_spatial_value
59 : : tenex::get_tensorindex_value_with_opposite_valence(
60 : next_lower_spatial_value);
61 : next_lower_spatial_value++;
62 : }
63 : }
64 :
65 : // assign the first NumIndicesToContract index values of the second operand so
66 : // that they will contract with the first NumIndicesToContract indices of the
67 : // first operand
68 : for (size_t i = 0; i < NumIndicesToContract; i++) {
69 : gsl::at(std::get<1>(tensor_index_values), i) =
70 : tenex::get_tensorindex_value_with_opposite_valence(
71 : gsl::at(std::get<0>(tensor_index_values), i));
72 : }
73 :
74 : // assign the remaining index values of the second operand so that they are
75 : // not duplicates of any previously-used indices
76 : for (size_t i = NumIndicesToContract; i < NumIndices2; i++) {
77 : const bool index2_is_spacetime = gsl::at(index_type_is_spacetime2, i);
78 : const bool index2_is_lower = gsl::at(valence_is_lower2, i);
79 : if (index2_is_spacetime) {
80 : gsl::at(std::get<1>(tensor_index_values), i) =
81 : index2_is_lower ? next_lower_spacetime_value
82 : : tenex::get_tensorindex_value_with_opposite_valence(
83 : next_lower_spacetime_value);
84 : next_lower_spacetime_value++;
85 : } else {
86 : gsl::at(std::get<1>(tensor_index_values), i) =
87 : index2_is_lower ? next_lower_spatial_value
88 : : tenex::get_tensorindex_value_with_opposite_valence(
89 : next_lower_spatial_value);
90 : next_lower_spatial_value++;
91 : }
92 : }
93 :
94 : // assign the free indices of the first operand to the result tensor indices
95 : for (size_t i = 0; i < NumIndices1 - NumIndicesToContract; i++) {
96 : gsl::at(std::get<2>(tensor_index_values), i) =
97 : gsl::at(std::get<0>(tensor_index_values), i + NumIndicesToContract);
98 : }
99 :
100 : // assign the free indices of the second operand to the result tensor indices
101 : for (size_t i = 0; i < NumIndices2 - NumIndicesToContract; i++) {
102 : gsl::at(std::get<2>(tensor_index_values),
103 : i + NumIndices1 - NumIndicesToContract) =
104 : gsl::at(std::get<1>(tensor_index_values), i + NumIndicesToContract);
105 : }
106 :
107 : return tensor_index_values;
108 : }
109 :
110 : // \brief Helper struct for contracting the first `NumIndicesToContract`
111 : // indices of two `Tensor`s
112 : //
113 : // \tparam NumIndicesToContract the number of indices to contract
114 : // \param T1 the type of the first `Tensor` of the two to contract
115 : // \param T2 the type of the second `Tensor` of the two to contract
116 : template <size_t NumIndicesToContract, typename T1, typename T2,
117 : size_t NumIndices1 = tmpl::size<typename T1::symmetry>::value,
118 : size_t NumIndices2 = tmpl::size<typename T2::symmetry>::value,
119 : size_t NumResultIndices =
120 : NumIndices1 + NumIndices2 - 2 * NumIndicesToContract,
121 : typename ContractedIndicesSeq =
122 : std::make_index_sequence<NumIndicesToContract>,
123 : typename TensorIndices1Seq = std::make_index_sequence<NumIndices1>,
124 : typename TensorIndices2Seq = std::make_index_sequence<NumIndices2>,
125 : typename ResultIndicesSeq =
126 : std::make_index_sequence<NumResultIndices>>
127 : struct contract_first_n_indices_impl;
128 :
129 : template <size_t NumIndicesToContract, typename X1, typename Symm1,
130 : typename... Indices1, typename X2, typename Symm2,
131 : typename... Indices2, size_t NumIndices1, size_t NumIndices2,
132 : size_t NumResultIndices, size_t... ContractedIndicesInts,
133 : size_t... Ints1, size_t... Ints2, size_t... ResultInts>
134 : struct contract_first_n_indices_impl<
135 : NumIndicesToContract, Tensor<X1, Symm1, tmpl::list<Indices1...>>,
136 : Tensor<X2, Symm2, tmpl::list<Indices2...>>, NumIndices1, NumIndices2,
137 : NumResultIndices, std::index_sequence<ContractedIndicesInts...>,
138 : std::index_sequence<Ints1...>, std::index_sequence<Ints2...>,
139 : std::index_sequence<ResultInts...>> {
140 : static_assert(NumIndicesToContract <= sizeof...(Indices1) and
141 : NumIndicesToContract <= sizeof...(Indices2),
142 : "Cannot request to contract more indices than indices in "
143 : "either of the two Tensors to contract.");
144 :
145 : static constexpr std::array<bool, NumIndices1> valence_is_lower1 = {
146 : {(Indices1::ul == UpLo::Lo)...}};
147 : static constexpr std::array<bool, NumIndices2> valence_is_lower2 = {
148 : {(Indices2::ul == UpLo::Lo)...}};
149 : static constexpr std::array<bool, NumIndices1> index_type_is_spacetime1 = {
150 : {(Indices1::index_type == IndexType::Spacetime)...}};
151 : static constexpr std::array<bool, NumIndices2> index_type_is_spacetime2 = {
152 : {(Indices2::index_type == IndexType::Spacetime)...}};
153 :
154 : // if this is removed and support for automatic contraction of a spatial index
155 : // with a spacetime index is wanted, the implementation of
156 : // get_tensor_index_values_for_tensors_to_contract_and_result_tensor() also
157 : // needs to be updated to support this
158 : static_assert(
159 : (... and (index_type_is_spacetime1[ContractedIndicesInts] ==
160 : index_type_is_spacetime2[ContractedIndicesInts])),
161 : "You are trying to automatically contract a spatial index with a "
162 : "spacetime index, but this is not supported for "
163 : "contract_first_n_indices().");
164 :
165 : // the values of the generic indices (i.e. `TensorIndex::value`s) that
166 : // uniquely identify different generic indices
167 : static constexpr auto tensor_index_values =
168 : get_tensor_index_values_for_tensors_to_contract_and_result_tensor<
169 : NumIndicesToContract>(index_type_is_spacetime1, valence_is_lower1,
170 : index_type_is_spacetime2, valence_is_lower2);
171 : static constexpr std::array<size_t, NumIndices1> tensor_index_values1 =
172 : std::get<0>(tensor_index_values);
173 : static constexpr std::array<size_t, NumIndices2> tensor_index_values2 =
174 : std::get<1>(tensor_index_values);
175 : static constexpr std::array<size_t, NumResultIndices>
176 : result_tensor_index_values = std::get<2>(tensor_index_values);
177 :
178 : using lhs_tensorindex_list =
179 : tmpl::list<TensorIndex<result_tensor_index_values[ResultInts]>...>;
180 :
181 : // \brief Contract first N indices by evaluating a `TensorExpression`
182 : //
183 : // \param lhs_tensor the result LHS `Tensor`
184 : // \param tensor1 the first `Tensor` of the two to contract
185 : // \param tensor2 the second `Tensor` of the two to contract
186 : template <typename LhsTensor>
187 : static void apply(const gsl::not_null<LhsTensor*> lhs_tensor,
188 : const Tensor<X1, Symm1, tmpl::list<Indices1...>>& tensor1,
189 : const Tensor<X2, Symm2, tmpl::list<Indices2...>>& tensor2) {
190 : constexpr bool evaluate_subtrees =
191 : decltype(tensor1(TensorIndex<tensor_index_values1[Ints1]>{}...) *
192 : tensor2(TensorIndex<tensor_index_values2[Ints2]>{}...))::
193 : primary_subtree_contains_primary_start;
194 : // Calls `evaluate_impl()` instead of `evaluate()` because `evaluate()`
195 : // takes `TensorIndex` lvalue references as template parameters, but here we
196 : // have `TensorIndex` types, and `evaluate()` cannot be overloaded to accept
197 : // both the former and the latter, so we simply circumvent the "top level"
198 : // `evaluate()` call that is normally used when evaluating the result of a
199 : // `TensorExpression`
200 : tenex::detail::evaluate_impl<
201 : evaluate_subtrees,
202 : TensorIndex<result_tensor_index_values[ResultInts]>...>(
203 : lhs_tensor,
204 : tensor1(TensorIndex<tensor_index_values1[Ints1]>{}...) *
205 : tensor2(TensorIndex<tensor_index_values2[Ints2]>{}...),
206 : std::make_index_sequence<NumResultIndices>{});
207 : }
208 :
209 : // \brief Contract first N indices by evaluating a `TensorExpression`
210 : //
211 : // \param tensor1 the first `Tensor` of the two to contract
212 : // \param tensor2 the second `Tensor` of the two to contract
213 : static auto apply(const Tensor<X1, Symm1, tmpl::list<Indices1...>>& tensor1,
214 : const Tensor<X2, Symm2, tmpl::list<Indices2...>>& tensor2) {
215 : using rhs_expression =
216 : decltype(tensor1(TensorIndex<tensor_index_values1[Ints1]>{}...) *
217 : tensor2(TensorIndex<tensor_index_values2[Ints2]>{}...));
218 : using rhs_tensorindex_list = typename rhs_expression::args_list;
219 : using rhs_symmetry = typename rhs_expression::symmetry;
220 : using rhs_tensorindextype_list = typename rhs_expression::index_list;
221 :
222 : using lhs_tensor_symm_and_indices =
223 : tenex::LhsTensorSymmAndIndices<rhs_tensorindex_list,
224 : lhs_tensorindex_list, rhs_symmetry,
225 : rhs_tensorindextype_list>;
226 :
227 : Tensor<typename rhs_expression::type,
228 : typename lhs_tensor_symm_and_indices::symmetry,
229 : typename lhs_tensor_symm_and_indices::tensorindextype_list>
230 : lhs_tensor{};
231 :
232 : apply(make_not_null(&lhs_tensor), tensor1, tensor2);
233 :
234 : return lhs_tensor;
235 : }
236 : };
237 : } // namespace detail
238 :
239 : /// \ingroup TensorGroup
240 : /// \brief Contract the first N indices of two `Tensor`s
241 : ///
242 : /// \details
243 : /// The indices of `lhs_tensor` should be the concatenation of the uncontracted
244 : /// indices of `tensor1` and the uncontracted indices of `tensor2`, in this
245 : /// order. For example, if `tensor1` is rank 3, `tensor2` is rank 4, and we want
246 : /// to contract the first two indices, the indices of `lhs_tensor` need to be
247 : /// the last index of `tensor1` followed by the 3rd index and then the 4th index
248 : /// of `tensor2`.
249 : ///
250 : /// The index types (spatial or spacetime) must be the same for the two indices
251 : /// in a pair of indices being contracted. Support can be added to this function
252 : /// to automatically contract the spatial indices of a spacetime index with a
253 : /// spatial index.
254 : ///
255 : /// \tparam NumIndicesToContract the number of indices to contract
256 : /// \param lhs_tensor the result LHS `Tensor`
257 : /// \param tensor1 the first `Tensor` of the two to contract
258 : /// \param tensor2 the second `Tensor` of the two to contract
259 : template <size_t NumIndicesToContract, typename LhsTensor, typename T1,
260 : typename T2>
261 1 : void contract_first_n_indices(const gsl::not_null<LhsTensor*> lhs_tensor,
262 : const T1& tensor1, const T2& tensor2) {
263 : detail::contract_first_n_indices_impl<NumIndicesToContract, T1, T2>::apply(
264 : lhs_tensor, tensor1, tensor2);
265 : }
266 :
267 : /// \ingroup TensorGroup
268 : /// \brief Contract the first N indices of two `Tensor`s
269 : ///
270 : /// \details
271 : /// The indices of the returned `Tensor` will be the concatenation of the
272 : /// uncontracted indices of `tensor1` and the uncontracted indices of `tensor2`,
273 : /// in this order. For example, if `tensor1` is rank 3, `tensor2` is rank 4, and
274 : /// we want to contract the first two indices, the indices of the returned
275 : /// `Tensor` will be the last index of `tensor1` followed by the 3rd index and
276 : /// then the 4th index of `tensor2`.
277 : ///
278 : /// The index types (spatial or spacetime) must be the same for the two indices
279 : /// in a pair of indices being contracted. Support can be added to this function
280 : /// to automatically contract the spatial indices of a spacetime index with a
281 : /// spatial index.
282 : ///
283 : /// \tparam NumIndicesToContract the number of indices to contract
284 : /// \param tensor1 the first `Tensor` of the two to contract
285 : /// \param tensor2 the second `Tensor` of the two to contract
286 : template <size_t NumIndicesToContract, typename T1, typename T2>
287 1 : auto contract_first_n_indices(const T1& tensor1, const T2& tensor2) {
288 : return detail::contract_first_n_indices_impl<NumIndicesToContract, T1,
289 : T2>::apply(tensor1, tensor2);
290 : }
|
