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, tensor1(TensorIndex<tensor_index_values1[Ints1]>{}...) *
204 : tensor2(TensorIndex<tensor_index_values2[Ints2]>{}...));
205 : }
206 :
207 : // \brief Contract first N indices by evaluating a `TensorExpression`
208 : //
209 : // \param tensor1 the first `Tensor` of the two to contract
210 : // \param tensor2 the second `Tensor` of the two to contract
211 : static auto apply(const Tensor<X1, Symm1, tmpl::list<Indices1...>>& tensor1,
212 : const Tensor<X2, Symm2, tmpl::list<Indices2...>>& tensor2) {
213 : using rhs_expression =
214 : decltype(tensor1(TensorIndex<tensor_index_values1[Ints1]>{}...) *
215 : tensor2(TensorIndex<tensor_index_values2[Ints2]>{}...));
216 : using rhs_tensorindex_list = typename rhs_expression::args_list;
217 : using rhs_symmetry = typename rhs_expression::symmetry;
218 : using rhs_tensorindextype_list = typename rhs_expression::index_list;
219 :
220 : using lhs_tensor_symm_and_indices =
221 : tenex::LhsTensorSymmAndIndices<rhs_tensorindex_list,
222 : lhs_tensorindex_list, rhs_symmetry,
223 : rhs_tensorindextype_list>;
224 :
225 : Tensor<typename rhs_expression::type,
226 : typename lhs_tensor_symm_and_indices::symmetry,
227 : typename lhs_tensor_symm_and_indices::tensorindextype_list>
228 : lhs_tensor{};
229 :
230 : apply(make_not_null(&lhs_tensor), tensor1, tensor2);
231 :
232 : return lhs_tensor;
233 : }
234 : };
235 : } // namespace detail
236 :
237 : /// \ingroup TensorGroup
238 : /// \brief Contract the first N indices of two `Tensor`s
239 : ///
240 : /// \details
241 : /// The indices of `lhs_tensor` should be the concatenation of the uncontracted
242 : /// indices of `tensor1` and the uncontracted indices of `tensor2`, in this
243 : /// order. For example, if `tensor1` is rank 3, `tensor2` is rank 4, and we want
244 : /// to contract the first two indices, the indices of `lhs_tensor` need to be
245 : /// the last index of `tensor1` followed by the 3rd index and then the 4th index
246 : /// of `tensor2`.
247 : ///
248 : /// The index types (spatial or spacetime) must be the same for the two indices
249 : /// in a pair of indices being contracted. Support can be added to this function
250 : /// to automatically contract the spatial indices of a spacetime index with a
251 : /// spatial index.
252 : ///
253 : /// \tparam NumIndicesToContract the number of indices to contract
254 : /// \param lhs_tensor the result LHS `Tensor`
255 : /// \param tensor1 the first `Tensor` of the two to contract
256 : /// \param tensor2 the second `Tensor` of the two to contract
257 : template <size_t NumIndicesToContract, typename LhsTensor, typename T1,
258 : typename T2>
259 1 : void contract_first_n_indices(const gsl::not_null<LhsTensor*> lhs_tensor,
260 : const T1& tensor1, const T2& tensor2) {
261 : detail::contract_first_n_indices_impl<NumIndicesToContract, T1, T2>::apply(
262 : lhs_tensor, tensor1, tensor2);
263 : }
264 :
265 : /// \ingroup TensorGroup
266 : /// \brief Contract the first N indices of two `Tensor`s
267 : ///
268 : /// \details
269 : /// The indices of the returned `Tensor` will be the concatenation of the
270 : /// uncontracted indices of `tensor1` and the uncontracted indices of `tensor2`,
271 : /// in this order. For example, if `tensor1` is rank 3, `tensor2` is rank 4, and
272 : /// we want to contract the first two indices, the indices of the returned
273 : /// `Tensor` will be the last index of `tensor1` followed by the 3rd index and
274 : /// then the 4th index of `tensor2`.
275 : ///
276 : /// The index types (spatial or spacetime) must be the same for the two indices
277 : /// in a pair of indices being contracted. Support can be added to this function
278 : /// to automatically contract the spatial indices of a spacetime index with a
279 : /// spatial index.
280 : ///
281 : /// \tparam NumIndicesToContract the number of indices to contract
282 : /// \param tensor1 the first `Tensor` of the two to contract
283 : /// \param tensor2 the second `Tensor` of the two to contract
284 : template <size_t NumIndicesToContract, typename T1, typename T2>
285 1 : auto contract_first_n_indices(const T1& tensor1, const T2& tensor2) {
286 : return detail::contract_first_n_indices_impl<NumIndicesToContract, T1,
287 : T2>::apply(tensor1, tensor2);
288 : }
|
