Line data Source code
1 1 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : /// \file
5 : /// Defines functions for evaluating `TensorExpression`s
6 :
7 : #pragma once
8 :
9 : #include <array>
10 : #include <complex>
11 : #include <cstddef>
12 : #include <type_traits>
13 :
14 : #include "DataStructures/ComplexDataVector.hpp"
15 : #include "DataStructures/DataVector.hpp"
16 : #include "DataStructures/Tensor/Expressions/DataTypeSupport.hpp"
17 : #include "DataStructures/Tensor/Expressions/IndexPropertyCheck.hpp"
18 : #include "DataStructures/Tensor/Expressions/LhsTensorSymmAndIndices.hpp"
19 : #include "DataStructures/Tensor/Expressions/SpatialSpacetimeIndex.hpp"
20 : #include "DataStructures/Tensor/Expressions/TensorExpression.hpp"
21 : #include "DataStructures/Tensor/Expressions/TensorIndex.hpp"
22 : #include "DataStructures/Tensor/Expressions/TensorIndexTransformation.hpp"
23 : #include "DataStructures/Tensor/Expressions/TimeIndex.hpp"
24 : #include "DataStructures/Tensor/Structure.hpp"
25 : #include "Utilities/Algorithm.hpp"
26 : #include "Utilities/ContainerHelpers.hpp"
27 : #include "Utilities/ErrorHandling/Assert.hpp"
28 : #include "Utilities/Gsl.hpp"
29 : #include "Utilities/Requires.hpp"
30 : #include "Utilities/TMPL.hpp"
31 :
32 : /// \cond
33 : template <typename X, typename Symm, typename IndexList>
34 : class Tensor;
35 : /// \endcond
36 :
37 : namespace tenex {
38 : namespace detail {
39 : template <size_t NumIndices>
40 : constexpr bool contains_indices_to_contract(
41 : const std::array<size_t, NumIndices>& tensorindices) {
42 : if constexpr (NumIndices < 2) {
43 : return false;
44 : } else {
45 : for (size_t i = 0; i < NumIndices - 1; i++) {
46 : for (size_t j = i + 1; j < NumIndices; j++) {
47 : const size_t current_tensorindex = gsl::at(tensorindices, i);
48 : // Concrete time indices are not contracted
49 : if ((not is_time_index_value(current_tensorindex)) and
50 : current_tensorindex == get_tensorindex_value_with_opposite_valence(
51 : gsl::at(tensorindices, j))) {
52 : return true;
53 : }
54 : }
55 : }
56 : return false;
57 : }
58 : }
59 :
60 : /// \brief Given the list of the positions of the LHS tensor's spacetime indices
61 : /// where a generic spatial index is used and the list of positions where a
62 : /// concrete time index is used, determine whether or not the component at the
63 : /// given LHS multi-index should be computed
64 : ///
65 : /// \details
66 : /// Not all of the LHS tensor's components may need to be computed. Cases when
67 : /// the component at a LHS multi-index should not be not evaluated:
68 : /// - If a generic spatial index is used for a spacetime index on the LHS,
69 : /// the components for which that index's concrete index is the time index
70 : /// should not be computed
71 : /// - If a concrete time index is used for a spacetime index on the LHS, the
72 : /// components for which that index's concrete index is a spatial index should
73 : /// not be computed
74 : ///
75 : /// \param lhs_multi_index the multi-index of the LHS tensor to check
76 : /// \param lhs_spatial_spacetime_index_positions the positions of the LHS
77 : /// tensor's spacetime indices where a generic spatial index is used
78 : /// \param lhs_time_index_positions the positions of the LHS tensor's spacetime
79 : /// indices where a concrete time index is used
80 : /// \return Whether or not `lhs_multi_index` is a multi-index of a component of
81 : /// the LHS tensor that should be computed
82 : template <size_t NumLhsIndices, size_t NumLhsSpatialSpacetimeIndices,
83 : size_t NumLhsConcreteTimeIndices>
84 : constexpr bool is_evaluated_lhs_multi_index(
85 : const std::array<size_t, NumLhsIndices>& lhs_multi_index,
86 : const std::array<size_t, NumLhsSpatialSpacetimeIndices>&
87 : lhs_spatial_spacetime_index_positions,
88 : const std::array<size_t, NumLhsConcreteTimeIndices>&
89 : lhs_time_index_positions) {
90 : for (size_t i = 0; i < lhs_spatial_spacetime_index_positions.size(); i++) {
91 : if (gsl::at(lhs_multi_index,
92 : gsl::at(lhs_spatial_spacetime_index_positions, i)) == 0) {
93 : return false;
94 : }
95 : }
96 : for (size_t i = 0; i < lhs_time_index_positions.size(); i++) {
97 : if (gsl::at(lhs_multi_index, gsl::at(lhs_time_index_positions, i)) != 0) {
98 : return false;
99 : }
100 : }
101 : return true;
102 : }
103 :
104 : template <typename SymmList>
105 : struct CheckNoLhsAntiSymmetries;
106 :
107 : template <template <typename...> class SymmList, typename... Symm>
108 : struct CheckNoLhsAntiSymmetries<SymmList<Symm...>> {
109 : static constexpr bool value = (... and (Symm::value > 0));
110 : };
111 :
112 : /// \brief Given a tensor and its list of tensor indices, return the
113 : /// canonicalized order of the tensor indices according to the tensor's symmetry
114 : ///
115 : /// \details
116 : /// The canonical ordering of a `Tensor`'s `TensorIndex`s
117 : /// (e.g. `ti::a`, `ti::b`, ti::c) is relevant to sets of indices that are
118 : /// symmetric. Within each set of symmetric indices, the `TensorIndex` used for
119 : /// each index can be freely reordered. Given a set of symmetric indices, this
120 : /// function defines the canonical order of the `TensorIndex`s assigned to them
121 : /// to be such that the lowest index positions take any generic spatial tensor
122 : /// indices (e.g. `ti::i`, `ti::j`), the next lowest index positions take any
123 : /// generic spacetime indices (e.g. `ti::a`, `ti::b`), and the highest index
124 : /// positions take any concrete time indices (e.g. `ti::t`, `ti::T`). We can
125 : /// imagine the canonical ordering of symmetric indices to look generally like:
126 : /// `[spatial indices ... | spacetime indices... | time indices...]`.
127 : ///
128 : /// Within the subsets of spatial, spacetime, and time indices, the
129 : /// `TensorIndex`s in each will be ordered such that lowercase indices come
130 : /// before uppercase, where both are ordered alphabetically. Another way of
131 : /// saying this is that if we had a rank N `Tensor` that was fully symmetric,
132 : /// its canonical ordering would take the following form:
133 : ///
134 : /// ```
135 : /// [ti::i, ti::j, ti::k, ..., ti::I, ti::J, ti::K, ..., // spatial
136 : /// ti::a, ti::b, ti::c, ..., ti::A, ti::B, ti::C, ..., // spacetime
137 : /// ti::t, ti::t, ti::t, ..., ti::T, ti::T, ti::T, ...] // time
138 : /// ```
139 : ///
140 : /// Here are some examples:
141 : ///
142 : /// ```
143 : /// symmetry: <1, 1, 1>
144 : /// set of `TensorIndex`s: {ti::t, ti::i, ti::a}
145 : /// canonical ordering: [ti::i, ti::a, ti::t]
146 : ///
147 : /// symmetry: <1, 1, 1>
148 : /// set of `TensorIndex`s: {ti::A, ti::a, ti::b}
149 : /// canonical ordering: [ti::a, ti::b, ti::A]
150 : /// ```
151 : ///
152 : /// When a `Tensor` is not fully symmetric, the `TensorIndex` labels for any
153 : /// indices that do not have symmetry with any other will simply keep their
154 : /// label because it cannot be swapped:
155 : ///
156 : /// ```
157 : /// symmetry: <1, 2, 1>
158 : /// set of `TensorIndex`s: {ti::a, ti::b, ti::i}
159 : /// canonical ordering: [ti::i, ti::b, ti::a]
160 : ///
161 : /// symmetry: <2, 1, 1>
162 : /// set of `TensorIndex`s: {ti::t, ti::k, ti::j}
163 : /// canonical ordering: [ti::t, ti::j, ti::k]
164 : /// ```
165 : ///
166 : /// If there is more than one set of symmetric indices, each of the subsets are
167 : /// individually reordered:
168 : ///
169 : /// ```
170 : /// symmetry: <1, 2, 2, 1>
171 : /// set of `TensorIndex`s: {ti::a, ti::b, ti::t, ti::i}
172 : /// canonical ordering: [ti::i, ti::b, ti::t, ti::a]
173 : /// ```
174 : ///
175 : /// The motivation for this specific canonical reordering is to quickly assess
176 : /// which components to assign to and which ones to skip when generic spatial
177 : /// and/or concrete time indices are used for symmetric spacetime indices in the
178 : /// resulting left hand side tensor when using `TensorExpression`s.
179 : ///
180 : /// Let's take the spacetime metric \f$g_{ab}\f$ as our motivating example. This
181 : /// tensor has symmetric spacetime indices, and let's say we only want to assign
182 : /// to \f$g_{ti}\f$. We want to loop over all 10 independent components of
183 : /// \f$g_{ab}\f$ and skip components outside of \f$g_{ti}\f$, e.g. \f$g_{xy}\f$
184 : /// (or \f$g_{12}\f$) and \f$g_{tt}\f$ (or \f$g_{00}\f$). To do so, when we see
185 : /// a multi-index like `{2, 1}` in our loop, we align `{2, 1}` with `{t, i}` and
186 : /// ask if the `2` is a valid index for `t` and if the `1` is a valid index for
187 : /// `i`. `1` is valid for `i`, but `2` is not valid for `t`, so we correctly
188 : /// skip over `{2, 1}` and don't assign to this component.
189 : ///
190 : /// However, this simple logic can lead to false positives or negatives when the
191 : /// indices are symmetric. What if the multi-index we're asking about is
192 : /// `{0, 1}` (\f$g_{01}\f$ or \f$g_{tx}\f$)? This logic would correctly
193 : /// determine that this is one of the components we want to assign to. But what
194 : /// if the multi-index was `{1, 0}`? The logic would incorrectly say to skip
195 : /// over and not assign to this multi-index because `1` is not a valid index for
196 : /// `t` and `0` is not valid for `i`. However, because \f$g_{ab}\f$ is
197 : /// symmetric, both `{0, 1}` and `{1, 0}` should give the same result, but
198 : /// `{1, 0}` gives us a false negative. Moreover and more generally, assigning
199 : /// to \f$g_{ti}\f$ and \f$g_{it}\f$ should yield the same behavior (assign to
200 : /// the same set of components).
201 : ///
202 : /// One way to address this would be to check the multi-indices for all
203 : /// permutations of symmetric index values, e.g. is `{0, 1}` *or* `{1, 0}`
204 : /// valid? And if so, then evaluate it. However, this adds work at runtime and
205 : /// the number of permutations to check increases as we increase the number of
206 : /// symmetric indices.
207 : ///
208 : /// The canonical reordering done by *this* function solves this problem by
209 : /// reordering the `TensorIndex`s to align nicely with the canonical multi-index
210 : /// ordering implemented by
211 : /// `::Tensor_detail::Structure::get_canonical_tensor_index`.
212 : /// `get_canonical_tensor_index` takes a storage index (which corresponds to an
213 : /// independent tensor component) and returns a canonical multi-index. This
214 : /// canonical multi-index is such that index values for symmetric indices will
215 : /// be ordered to increase from right to left. For example, for a rank 2
216 : /// symmetric tensor, `{1, 0}` is the canonical multi-index corresponding to the
217 : /// dependent multi-indices `{0, 1}` and `{1, 0}`. Therefore, `{1, 0}` would be
218 : /// the multi-index returned by `get_canonical_tensor_index` that corresponds to
219 : /// the single independent component. In other words, when we loop over the
220 : /// independent canonical multi-indices, we are looping over the multi-index
221 : /// permutations that are in the lower triangle of the N-dimensional matrix
222 : /// containing all multi-index permutations. The canonical reordering of LHS
223 : /// `TensorIndex`s for symmetric indices that is done by *this* function is
224 : /// implemented to match this: by making time indices the rightmost, then
225 : /// spacetime the next rightmost, and then spatial indices leftmost, we
226 : /// guarantee that looping over the lower triangle permutations given by
227 : /// `get_canonical_tensor_index` will not produce false positives or negatives
228 : /// using the earlier simple logic to check for valid multi-indices.
229 : ///
230 : /// We can use the spacetime metric as an example to demonstrate this. The
231 : /// lower triangle multi-indices that are looped over are ordered with index
232 : /// values increasing right to left, e.g. `{0, 0}`, `{1, 0}`, `{2, 0}`,
233 : /// `{2, 1}`, etc. If a user wants to assign to \f$g_{ti}\f$, then after this
234 : /// function internally reorders the LHS indices to \f$g_{it}\f$, when we loop
235 : /// and encounter `{1, 0}`, we correctly get that we should evaluate this
236 : /// component without having to check its other permutation. Likewise, if a user
237 : /// wants to assign to \f$g_{it}\f$, no reordering is done and we get the same
238 : /// correct behavior.
239 : ///
240 : /// This function works in general because:
241 : /// - any `0`s in the symmetric indices will "first be dealt" to any time
242 : /// `TensorIndex`s in the rightmost index positions and then any spacetime
243 : /// `TensorIndex`s, where `0` is correctly valid for both, but
244 : /// - if there are more `0`s than time + spacetime `TensorIndex`s, they will be
245 : /// "dealt" to spatial indices, which is always correctly invalid, and
246 : /// - if there are more time `TensorIndex`s than `0`s, values `> 0` will be
247 : /// "dealt" to time indices, which is also always correctly invalid.
248 : ///
249 : /// In this way, we don't ever have to check other permutations of sets of
250 : /// symmetric index values.
251 : ///
252 : /// \tparam LhsTensorIndices the `TensorIndex`s of the `Tensor`, e.g. `ti::a`,
253 : /// `ti::b`, `ti::c`
254 : /// \param canonical_symmetry the canonicalized symmetry values of the tensor
255 : /// (see `Symmetry` for definition of the canonical ordering of symmetry values)
256 : /// \return reordered values of `LhsTensorIndices::value...`
257 : template <typename... LhsTensorIndices, size_t NumIndices>
258 : constexpr std::array<size_t, NumIndices> get_reordered_tensorindex_values(
259 : const std::array<std::int32_t, NumIndices>& canonical_symmetry) {
260 : constexpr std::array<size_t, NumIndices> lhs_tensorindex_values = {
261 : {LhsTensorIndices::value...}};
262 : if constexpr (NumIndices < 2) {
263 : return lhs_tensorindex_values;
264 : } else {
265 : const auto compare = [](const size_t tensorindex_value1,
266 : const size_t tensorindex_value2) {
267 : // clang-tidy thinks these two branches are the same but they aren't:
268 : // if (tensorindex_value2 == ti::T.value)
269 : // if (tensorindex_value2 == ti::t.value)
270 : // NOLINTNEXTLINE (clang-tidy: bugprone-branch-clone)
271 : if (tensorindex_value2 == ti::T.value) {
272 : return false;
273 : } else if (is_time_index_value(tensorindex_value1)) {
274 : return true;
275 : } else if (tensorindex_value2 == ti::t.value) {
276 : return false;
277 : }
278 :
279 : return (is_generic_spacetime_index_value(tensorindex_value1) and
280 : is_generic_spatial_index_value(tensorindex_value2)) or
281 : (tensorindex_value1 > tensorindex_value2 and
282 : is_generic_spacetime_index_value(tensorindex_value1) ==
283 : is_generic_spacetime_index_value(tensorindex_value2));
284 : };
285 :
286 : std::array<size_t, NumIndices> reordered_lhs_tensorindex_values =
287 : lhs_tensorindex_values;
288 :
289 : std::int32_t max_symm_value = *alg::max_element(canonical_symmetry);
290 : std::int32_t symm_value_to_find = 1;
291 : while (symm_value_to_find <= max_symm_value) {
292 : size_t i = NumIndices - 1;
293 : while (true) {
294 : // skip forward until we get to the position with the value we want
295 : while (i > 0 and canonical_symmetry[i] != symm_value_to_find) {
296 : i--;
297 : }
298 : if (i == 0) {
299 : break;
300 : }
301 :
302 : size_t max_tensorindex_value = reordered_lhs_tensorindex_values[i];
303 : size_t max_index = i;
304 :
305 : size_t j = i - 1;
306 : // note: because we need to hit 0 and size_t wraps around to max size_t
307 : while (j < NumIndices) {
308 : const std::int32_t compare_symm_value = canonical_symmetry[j];
309 : const size_t compare_tensorindex_value =
310 : reordered_lhs_tensorindex_values[j];
311 : if (compare_symm_value == symm_value_to_find and
312 : compare(compare_tensorindex_value, max_tensorindex_value)) {
313 : max_tensorindex_value = compare_tensorindex_value;
314 : max_index = j;
315 : }
316 : j--;
317 : }
318 : reordered_lhs_tensorindex_values[max_index] =
319 : reordered_lhs_tensorindex_values[i];
320 : reordered_lhs_tensorindex_values[i] = max_tensorindex_value;
321 : i--;
322 : }
323 : symm_value_to_find++;
324 : }
325 :
326 : return reordered_lhs_tensorindex_values;
327 : }
328 : }
329 :
330 : /*!
331 : * \ingroup TensorExpressionsGroup
332 : * \brief Evaluate subtrees of the RHS expression or the RHS expression as a
333 : * whole and assign the result to the LHS tensor
334 : *
335 : * \details This is for internal use only and should never be directly called.
336 : * See `tenex::evaluate` and use it, instead.
337 : *
338 : * `EvaluateSubtrees` controls whether we wish to evaluate RHS subtrees or the
339 : * entire RHS expression as one expression. See`TensorExpression` documentation
340 : * on equation splitting for more details on what this means.
341 : *
342 : * If `EvaluateSubtrees == false`, then it's safe if the LHS tensor is used in
343 : * the RHS expression, so long as the generic index orders are the same. This
344 : * means that the callee of this function needs to first verify this is true
345 : * before calling this function. Under these conditions, this is a safe
346 : * operation because the implementation modifies each LHS component once and
347 : * does not revisit and access any LHS components after they've been updated.
348 : * For example, say we do `tenex::evaluate<ti_a, ti_b>(make_not_null(&L),
349 : * 5.0 * L(ti_a, ti_b));`. This function will first compute the RHS for some
350 : * concrete LHS, e.g. \f$L_{00}\f$. To compute this, it accesses \f$L_{00}\f$
351 : * in the RHS tree, multiplies it by `5.0`, then updates \f$L_{00}\f$ to be the
352 : * result of this multiplication. Next, it might compute \f$L_{01}\f$, where
353 : * only \f$L_{01}\f$ is accessed, and which hasn't yet been modified. Then the
354 : * next component is computed and updated, and so forth. These steps are
355 : * performed once for each unique LHS index. Therefore, it is important to note
356 : * that this kind of operation being safe to perform is
357 : * implementation-dependent. Specifically, the safety of the operation depends
358 : * on the order of LHS component access and assignment.
359 : *
360 : * \note `LhsTensorIndices` must be passed by reference because non-type
361 : * template parameters cannot be class types until C++20.
362 : *
363 : * @tparam EvaluateSubtrees whether or not to evaluate subtrees of RHS
364 : * expression
365 : * @tparam LhsTensorIndices the `TensorIndex`s of the `Tensor` on the LHS of the
366 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
367 : * @param lhs_tensor pointer to the resultant LHS `Tensor` to fill
368 : * @param rhs_tensorexpression the RHS TensorExpression to be evaluated
369 : */
370 : template <bool EvaluateSubtrees, typename... LhsTensorIndices,
371 : typename LhsDataType, typename LhsSymmetry, typename LhsIndexList,
372 : typename Derived, typename RhsDataType, typename RhsSymmetry,
373 : typename RhsIndexList, typename... RhsTensorIndices,
374 : size_t... LhsInts>
375 : void evaluate_impl(
376 : const gsl::not_null<Tensor<LhsDataType, LhsSymmetry, LhsIndexList>*>
377 : lhs_tensor,
378 : const TensorExpression<Derived, RhsDataType, RhsSymmetry, RhsIndexList,
379 : tmpl::list<RhsTensorIndices...>>&
380 : rhs_tensorexpression,
381 : const std::index_sequence<LhsInts...>& /*lhs_ints*/) {
382 : constexpr size_t num_lhs_indices = sizeof...(LhsTensorIndices);
383 : constexpr size_t num_rhs_indices = sizeof...(RhsTensorIndices);
384 :
385 : using lhs_tensorindex_list = tmpl::list<LhsTensorIndices...>;
386 : using rhs_tensorindex_list = tmpl::list<RhsTensorIndices...>;
387 :
388 : using lhs_tensor_type = typename std::decay_t<decltype(*lhs_tensor)>;
389 :
390 : static_assert(is_supported_tensor_datatype_v<LhsDataType> and
391 : is_supported_tensor_datatype_v<RhsDataType>,
392 : "TensorExpressions currently only support Tensors whose data "
393 : "type is double, std::complex<double>, DataVector, or "
394 : "ComplexDataVector. It is possible to add support for other "
395 : "data types that are supported by Tensor.");
396 : static_assert(
397 : is_assignable_v<LhsDataType, RhsDataType>,
398 : "Assignment of the LHS Tensor's data type to the RHS TensorExpression's "
399 : "data type is not supported. This happens from doing something like e.g. "
400 : "trying to assign a Tensor<double> to a Tensor<DataVector> or a "
401 : "Tensor<DataVector> to a Tensor<ComplexDataVector>.");
402 : // `Symmetry` currently prevents this because antisymmetries are not currently
403 : // supported for `Tensor`s. This check is repeated here because if
404 : // antisymmetries are later supported for `Tensor`, using antisymmetries in
405 : // `TensorExpression`s will not automatically work. The implementations of the
406 : // derived `TensorExpression` types assume no antisymmetries (assume positive
407 : // `Symmetry` values), so support for antisymmetries in `TensorExpression`s
408 : // will still need to be implemented.
409 : static_assert(CheckNoLhsAntiSymmetries<LhsSymmetry>::value,
410 : "Anti-symmetric Tensors are not currently supported by "
411 : "TensorExpressions.");
412 : static_assert(
413 : tmpl::equal_members<
414 : typename remove_time_indices<lhs_tensorindex_list>::type,
415 : typename remove_time_indices<rhs_tensorindex_list>::type>::value,
416 : "The generic indices on the LHS of a tensor equation (that is, the "
417 : "template parameters specified in evaluate<...>) must match the generic "
418 : "indices of the RHS TensorExpression. This error occurs as a result of a "
419 : "call like evaluate<ti::a, ti::b>(R(ti::A, ti::b) * S(ti::a, ti::c)), "
420 : "where the generic indices of the evaluated RHS expression are ti::b and "
421 : "ti::c, but the generic indices provided for the LHS are ti::a and "
422 : "ti::b.");
423 : static_assert(
424 : tensorindex_list_is_valid<lhs_tensorindex_list>::value,
425 : "Cannot assign a tensor expression to a LHS tensor with a repeated "
426 : "generic index, e.g. evaluate<ti::a, ti::a>. (Note that the concrete "
427 : "time indices (ti::T and ti::t) can be repeated.)");
428 : static_assert(
429 : not contains_indices_to_contract<num_lhs_indices>(
430 : {{LhsTensorIndices::value...}}),
431 : "Cannot assign a tensor expression to a LHS tensor with generic "
432 : "indices that would be contracted, e.g. evaluate<ti::A, ti::a>.");
433 : // `IndexPropertyCheck` does also check that valence (Up/Lo) of indices that
434 : // correspond in the RHS and LHS tensors are equal, but the assertion message
435 : // below does not mention this because a mismatch in valence should have been
436 : // caught due to the combination of (i) the Tensor::operator() assertion
437 : // checking that generic indices' valences match the tensor's indices'
438 : // valences and (ii) the above assertion that RHS and LHS generic indices
439 : // match
440 : static_assert(
441 : IndexPropertyCheck<LhsIndexList, RhsIndexList, lhs_tensorindex_list,
442 : rhs_tensorindex_list>::value,
443 : "At least one index of the tensor evaluated from the RHS expression "
444 : "cannot be evaluated to its corresponding index in the LHS tensor. This "
445 : "is due to a difference in number of spatial dimensions or Frame type "
446 : "between the index on the RHS and LHS. "
447 : "e.g. evaluate<ti::a, ti::b>(L, R(ti::b, ti::a));, where R's first "
448 : "index has 2 spatial dimensions but L's second index has 3 spatial "
449 : "dimensions. Check RHS and LHS indices that use the same generic index.");
450 : static_assert(Derived::height_relative_to_closest_tensor_leaf_in_subtree <
451 : std::numeric_limits<size_t>::max(),
452 : "Either no Tensors were found in the RHS TensorExpression or "
453 : "the depth of the tree exceeded the maximum size_t value (very "
454 : "unlikely). If there is indeed a Tensor in the RHS expression "
455 : "and assuming the tree's height is not actually the maximum "
456 : "size_t value, then there is a flaw in the logic for computing "
457 : "the derived TensorExpression types' member, "
458 : "height_relative_to_closest_tensor_leaf_in_subtree.");
459 :
460 : if constexpr (EvaluateSubtrees) {
461 : // Make sure the LHS tensor doesn't also appear in the RHS tensor expression
462 : (~rhs_tensorexpression).assert_lhs_tensor_not_in_rhs_expression(lhs_tensor);
463 : // If the LHS data type is a vector, size the LHS tensor components if their
464 : // size does not match the size from a `Tensor` in the RHS expression
465 : if constexpr (is_derived_of_vector_impl_v<LhsDataType>) {
466 : const size_t rhs_component_size =
467 : (~rhs_tensorexpression).get_rhs_tensor_component_size();
468 : if (rhs_component_size != (*lhs_tensor)[0].size()) {
469 : for (auto& lhs_component : *lhs_tensor) {
470 : lhs_component = LhsDataType(rhs_component_size);
471 : }
472 : }
473 : }
474 : }
475 :
476 : constexpr std::array<std::int32_t, num_lhs_indices> lhs_symmetry = {
477 : {tmpl::at_c<LhsSymmetry, LhsInts>::value...}};
478 : constexpr std::array<size_t, num_lhs_indices> reordered_tensorindex_values =
479 : get_reordered_tensorindex_values<LhsTensorIndices...>(lhs_symmetry);
480 : using reordered_lhs_tensorindex_list =
481 : tmpl::list<TensorIndex<reordered_tensorindex_values[LhsInts]>...>;
482 :
483 : constexpr std::array<size_t, num_rhs_indices> index_transformation =
484 : compute_tensorindex_transformation<num_lhs_indices, num_rhs_indices>(
485 : reordered_tensorindex_values, {{RhsTensorIndices::value...}});
486 :
487 : // positions of indices in LHS tensor where generic spatial indices are used
488 : // for spacetime indices
489 : constexpr auto lhs_spatial_spacetime_index_positions =
490 : get_spatial_spacetime_index_positions<LhsIndexList,
491 : reordered_lhs_tensorindex_list>();
492 : // positions of indices in RHS tensor where generic spatial indices are used
493 : // for spacetime indices
494 : constexpr auto rhs_spatial_spacetime_index_positions =
495 : get_spatial_spacetime_index_positions<RhsIndexList,
496 : rhs_tensorindex_list>();
497 :
498 : // positions of indices in LHS tensor where concrete time indices are used
499 : constexpr auto lhs_time_index_positions =
500 : get_time_index_positions<reordered_lhs_tensorindex_list>();
501 :
502 : using rhs_expression_type =
503 : typename std::decay_t<decltype(~rhs_tensorexpression)>;
504 :
505 : for (size_t i = 0; i < lhs_tensor_type::size(); i++) {
506 : auto lhs_multi_index =
507 : lhs_tensor_type::structure::get_canonical_tensor_index(i);
508 : if (is_evaluated_lhs_multi_index(lhs_multi_index,
509 : lhs_spatial_spacetime_index_positions,
510 : lhs_time_index_positions)) {
511 : for (size_t j = 0; j < lhs_spatial_spacetime_index_positions.size();
512 : j++) {
513 : gsl::at(lhs_multi_index,
514 : gsl::at(lhs_spatial_spacetime_index_positions, j)) -= 1;
515 : }
516 : auto rhs_multi_index =
517 : transform_multi_index(lhs_multi_index, index_transformation);
518 : for (size_t j = 0; j < rhs_spatial_spacetime_index_positions.size();
519 : j++) {
520 : gsl::at(rhs_multi_index,
521 : gsl::at(rhs_spatial_spacetime_index_positions, j)) += 1;
522 : }
523 :
524 : // The expression will either be evaluated as one whole expression
525 : // or it will be split up into subtrees that are evaluated one at a time.
526 : // See the section on splitting in the documentation for the
527 : // `TensorExpression` class to understand the logic and terminology used
528 : // in this control flow below.
529 : if constexpr (EvaluateSubtrees) {
530 : // the expression is split up, so evaluate subtrees at splits
531 : (~rhs_tensorexpression)
532 : .evaluate_primary_subtree((*lhs_tensor)[i], rhs_multi_index);
533 : if constexpr (not rhs_expression_type::is_primary_start) {
534 : // the root expression type is not the starting point of a leg, so it
535 : // has not yet been evaluated, so now we evaluate this last leg of the
536 : // expression at the root of the tree
537 : (*lhs_tensor)[i] =
538 : (~rhs_tensorexpression)
539 : .get_primary((*lhs_tensor)[i], rhs_multi_index);
540 : }
541 : } else {
542 : // the expression is not split up, so evaluate full expression
543 : (*lhs_tensor)[i] = (~rhs_tensorexpression).get(rhs_multi_index);
544 : }
545 : }
546 : }
547 : }
548 :
549 : /*!
550 : * \ingroup TensorExpressionsGroup
551 : * \brief Assign a value to components of the LHS tensor
552 : *
553 : * \details This is for internal use only and should never be directly called.
554 : * See `tenex::evaluate` and use it, instead.
555 : *
556 : * \note `LhsTensorIndices` must be passed by reference because non-type
557 : * template parameters cannot be class types until C++20.
558 : *
559 : * @tparam LhsTensorIndices the `TensorIndex`s of the `Tensor` on the LHS of the
560 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
561 : * @param lhs_tensor pointer to the resultant LHS `Tensor` to fill
562 : * @param rhs_value the RHS value to assigned
563 : */
564 : template <typename... LhsTensorIndices, typename X, typename LhsSymmetry,
565 : typename LhsIndexList, typename NumberType, size_t... LhsInts>
566 : void evaluate_impl(
567 : const gsl::not_null<Tensor<X, LhsSymmetry, LhsIndexList>*> lhs_tensor,
568 : const NumberType& rhs_value,
569 : const std::index_sequence<LhsInts...>& /*lhs_ints*/) {
570 : using lhs_tensor_type = typename std::decay_t<decltype(*lhs_tensor)>;
571 : constexpr size_t num_lhs_indices = sizeof...(LhsTensorIndices);
572 : using lhs_tensorindex_list = tmpl::list<LhsTensorIndices...>;
573 :
574 : static_assert(is_supported_tensor_datatype_v<X> and
575 : "TensorExpressions currently only support Tensors whose data "
576 : "type is double, std::complex<double>, DataVector, or "
577 : "ComplexDataVector. It is possible to add support for other "
578 : "data types that are supported by Tensor.");
579 : static_assert(
580 : is_assignable_v<X, NumberType>,
581 : "Assignment of the LHS Tensor's data type to the RHS number's data type "
582 : "is not supported within TensorExpressions. This happens from doing "
583 : "something like e.g. trying to assign a double to a DataVector or a "
584 : "DataVector to a ComplexDataVector.");
585 : // `Symmetry` currently prevents this because antisymmetries are not currently
586 : // supported for `Tensor`s. This check is repeated here because if
587 : // antisymmetries are later supported for `Tensor`, using antisymmetries in
588 : // `TensorExpression`s will not automatically work. The implementations of the
589 : // derived `TensorExpression` types assume no antisymmetries (assume positive
590 : // `Symmetry` values), so support for antisymmetries in `TensorExpression`s
591 : // will still need to be implemented.
592 : static_assert(CheckNoLhsAntiSymmetries<LhsSymmetry>::value,
593 : "Anti-symmetric Tensors are not currently supported by "
594 : "TensorExpressions.");
595 : static_assert(
596 : tensorindex_list_is_valid<lhs_tensorindex_list>::value,
597 : "Cannot assign a tensor expression to a LHS tensor with a repeated "
598 : "generic index, e.g. evaluate<ti::a, ti::a>. (Note that the concrete "
599 : "time indices (ti::T and ti::t) can be repeated.)");
600 : static_assert(
601 : not contains_indices_to_contract<num_lhs_indices>(
602 : {{LhsTensorIndices::value...}}),
603 : "Cannot assign a tensor expression to a LHS tensor with generic "
604 : "indices that would be contracted, e.g. evaluate<ti::A, ti::a>.");
605 :
606 : if constexpr (is_derived_of_vector_impl_v<X>) {
607 : ASSERT(get_size((*lhs_tensor)[0]) > 0,
608 : "Tensors with vector components must be sized before calling "
609 : "\ntenex::evaluate<...>("
610 : "\n\tgsl::not_null<Tensor<VectorType, ...>*>, number).");
611 : }
612 :
613 : constexpr std::array<std::int32_t, num_lhs_indices> lhs_symmetry = {
614 : {tmpl::at_c<LhsSymmetry, LhsInts>::value...}};
615 : constexpr std::array<size_t, num_lhs_indices> reordered_tensorindex_values =
616 : get_reordered_tensorindex_values<LhsTensorIndices...>(lhs_symmetry);
617 : (void)reordered_tensorindex_values; // silence false unused variable warning
618 : using reordered_lhs_tensorindex_list =
619 : tmpl::list<TensorIndex<reordered_tensorindex_values[LhsInts]>...>;
620 :
621 : // positions of indices in LHS tensor where generic spatial indices are used
622 : // for spacetime indices
623 : constexpr auto lhs_spatial_spacetime_index_positions =
624 : get_spatial_spacetime_index_positions<LhsIndexList,
625 : reordered_lhs_tensorindex_list>();
626 :
627 : // positions of indices in LHS tensor where concrete time indices are used
628 : constexpr auto lhs_time_index_positions =
629 : get_time_index_positions<reordered_lhs_tensorindex_list>();
630 :
631 : for (size_t i = 0; i < lhs_tensor_type::size(); i++) {
632 : auto lhs_multi_index =
633 : lhs_tensor_type::structure::get_canonical_tensor_index(i);
634 : if (is_evaluated_lhs_multi_index(lhs_multi_index,
635 : lhs_spatial_spacetime_index_positions,
636 : lhs_time_index_positions)) {
637 : (*lhs_tensor)[i] = rhs_value;
638 : }
639 : }
640 : }
641 : } // namespace detail
642 :
643 : /*!
644 : * \ingroup TensorExpressionsGroup
645 : * \brief Assign the result of a RHS tensor expression to a tensor with the LHS
646 : * index order set in the template parameters
647 : *
648 : * \details Uses the right hand side (RHS) TensorExpression's index ordering
649 : * (`RhsTE::args_list`) and the desired left hand side (LHS) tensor's index
650 : * ordering (`LhsTensorIndices`) to fill the provided LHS Tensor with that LHS
651 : * index ordering. This can carry out the evaluation of a RHS tensor expression
652 : * to a LHS tensor with the same index ordering, such as \f$L_{ab} = R_{ab}\f$,
653 : * or different ordering, such as \f$L_{ba} = R_{ab}\f$.
654 : *
655 : * The symmetry of the provided LHS Tensor need not match the symmetry
656 : * determined from evaluating the RHS TensorExpression according to its order of
657 : * operations. This allows one to specify LHS symmetries (via `lhs_tensor`) that
658 : * may not be preserved by the RHS expression's order of operations, which
659 : * depends on how the expression is written and implemented.
660 : *
661 : * The LHS `Tensor` cannot be part of the RHS expression, e.g.
662 : * `evaluate(make_not_null(&L), L() + R());`, because the LHS `Tensor` will
663 : * generally not be computed correctly when the RHS `TensorExpression` is split
664 : * up and the LHS tensor components are computed by accumulating the result of
665 : * subtrees (see the section on splitting in the documentation for the
666 : * `TensorExpression` class). If you need to use the LHS `Tensor` on the RHS,
667 : * use `tenex::update` instead.
668 : *
669 : * ### Example usage
670 : * Given `Tensor`s `R`, `S`, `T`, `G`, and `H`, we can compute the LHS tensor
671 : * \f$L\f$ in the equation \f$L_{a} = R_{ab} S^{b} + G_{a} - H_{ba}{}^{b} T\f$
672 : * by doing:
673 : *
674 : * \snippet Test_MixedOperations.cpp use_evaluate_with_result_as_arg
675 : *
676 : * \note `LhsTensorIndices` must be passed by reference because non-type
677 : * template parameters cannot be class types until C++20.
678 : *
679 : * @tparam LhsTensorIndices the `TensorIndex`s of the `Tensor` on the LHS of the
680 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
681 : * @param lhs_tensor pointer to the resultant LHS `Tensor` to fill
682 : * @param rhs_tensorexpression the RHS TensorExpression to be evaluated
683 : */
684 : template <auto&... LhsTensorIndices, typename LhsDataType, typename LhsSymmetry,
685 : typename LhsIndexList, typename Derived, typename RhsDataType,
686 : typename RhsSymmetry, typename RhsIndexList,
687 : typename... RhsTensorIndices>
688 1 : void evaluate(
689 : const gsl::not_null<Tensor<LhsDataType, LhsSymmetry, LhsIndexList>*>
690 : lhs_tensor,
691 : const TensorExpression<Derived, RhsDataType, RhsSymmetry, RhsIndexList,
692 : tmpl::list<RhsTensorIndices...>>&
693 : rhs_tensorexpression) {
694 : using rhs_expression_type =
695 : typename std::decay_t<decltype(~rhs_tensorexpression)>;
696 : constexpr bool evaluate_subtrees =
697 : rhs_expression_type::primary_subtree_contains_primary_start;
698 : detail::evaluate_impl<evaluate_subtrees,
699 : std::decay_t<decltype(LhsTensorIndices)>...>(
700 : lhs_tensor, rhs_tensorexpression,
701 : std::make_index_sequence<sizeof...(LhsTensorIndices)>{});
702 : }
703 :
704 : /// @{
705 : /*!
706 : * \ingroup TensorExpressionsGroup
707 : * \brief Assign a number to components of a tensor with the LHS index order
708 : * set in the template parameters
709 : *
710 : * \details
711 : * Example usage:
712 : * \snippet Test_MixedOperations.cpp assign_double_to_index_subsets
713 : *
714 : * \note The components of the LHS `Tensor` passed in must already be sized
715 : * because there is no way to infer component size from the RHS
716 : *
717 : * \note `LhsTensorIndices` must be passed by reference because non-type
718 : * template parameters cannot be class types until C++20.
719 : *
720 : * @tparam LhsTensorIndices the `TensorIndex`s of the `Tensor` on the LHS of the
721 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
722 : * @param lhs_tensor pointer to the resultant LHS `Tensor` to fill
723 : * @param rhs_value the RHS value to assign
724 : */
725 : template <auto&... LhsTensorIndices, typename X, typename LhsSymmetry,
726 : typename LhsIndexList, typename N,
727 : Requires<std::is_arithmetic_v<N>> = nullptr>
728 1 : void evaluate(
729 : const gsl::not_null<Tensor<X, LhsSymmetry, LhsIndexList>*> lhs_tensor,
730 : const N rhs_value) {
731 : detail::evaluate_impl<std::decay_t<decltype(LhsTensorIndices)>...>(
732 : lhs_tensor, rhs_value,
733 : std::make_index_sequence<sizeof...(LhsTensorIndices)>{});
734 : }
735 : template <auto&... LhsTensorIndices, typename X, typename LhsSymmetry,
736 : typename LhsIndexList, typename N>
737 1 : void evaluate(
738 : const gsl::not_null<Tensor<X, LhsSymmetry, LhsIndexList>*> lhs_tensor,
739 : const std::complex<N>& rhs_value) {
740 : detail::evaluate_impl<std::decay_t<decltype(LhsTensorIndices)>...>(
741 : lhs_tensor, rhs_value,
742 : std::make_index_sequence<sizeof...(LhsTensorIndices)>{});
743 : }
744 : /// @}
745 :
746 : /*!
747 : * \ingroup TensorExpressionsGroup
748 : * \brief Assign the result of a RHS tensor expression to a tensor with the LHS
749 : * index order set in the template parameters
750 : *
751 : * \details Uses the right hand side (RHS) TensorExpression's index ordering
752 : * (`RhsTE::args_list`) and the desired left hand side (LHS) tensor's index
753 : * ordering (`LhsTensorIndices`) to construct a LHS Tensor with that LHS index
754 : * ordering. This can carry out the evaluation of a RHS tensor expression to a
755 : * LHS tensor with the same index ordering, such as \f$L_{ab} = R_{ab}\f$, or
756 : * different ordering, such as \f$L_{ba} = R_{ab}\f$.
757 : *
758 : * The symmetry of the returned LHS Tensor depends on the order of operations in
759 : * the RHS TensorExpression, i.e. how the expression is written. If you would
760 : * like to specify the symmetry of the LHS Tensor instead of it being determined
761 : * by the order of operations in the RHS expression, please use the other
762 : * `tenex::evaluate` overload that takes an empty LHS Tensor as its first
763 : * argument.
764 : *
765 : * ### Example usage
766 : * Given `Tensor`s `R`, `S`, `T`, `G`, and `H`, we can compute the LHS tensor
767 : * \f$L\f$ in the equation \f$L_{a} = R_{ab} S^{b} + G_{a} - H_{ba}{}^{b} T\f$
768 : * by doing:
769 : *
770 : * \snippet Test_MixedOperations.cpp use_evaluate_to_return_result
771 : *
772 : * \parblock
773 : * \note If a generic spatial index is used for a spacetime index in the RHS
774 : * tensor, its corresponding index in the LHS tensor type will be a spatial
775 : * index with the same valence, frame, and number of spatial dimensions. If a
776 : * concrete time index is used for a spacetime index in the RHS tensor, the
777 : * index will not appear in the LHS tensor (i.e. there will NOT be a
778 : * corresponding LHS index where only the time index of that index has been
779 : * computed and its spatial indices are empty).
780 : * \endparblock
781 : *
782 : * \parblock
783 : * \note `LhsTensorIndices` must be passed by reference because non-type
784 : * template parameters cannot be class types until C++20.
785 : * \endparblock
786 : *
787 : * @tparam LhsTensorIndices the TensorIndexs of the Tensor on the LHS of the
788 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
789 : * @param rhs_tensorexpression the RHS TensorExpression to be evaluated
790 : * @return the resultant LHS Tensor with index order specified by
791 : * LhsTensorIndices
792 : */
793 : template <auto&... LhsTensorIndices, typename RhsTE,
794 : Requires<std::is_base_of_v<Expression, RhsTE>> = nullptr>
795 1 : auto evaluate(const RhsTE& rhs_tensorexpression) {
796 : using lhs_tensorindex_list =
797 : tmpl::list<std::decay_t<decltype(LhsTensorIndices)>...>;
798 : using rhs_tensorindex_list = typename RhsTE::args_list;
799 : using rhs_symmetry = typename RhsTE::symmetry;
800 : using rhs_tensorindextype_list = typename RhsTE::index_list;
801 :
802 : // Stores (potentially reordered) symmetry and indices needed for constructing
803 : // the LHS tensor, with index order specified by LhsTensorIndices
804 : using lhs_tensor_symm_and_indices =
805 : LhsTensorSymmAndIndices<rhs_tensorindex_list, lhs_tensorindex_list,
806 : rhs_symmetry, rhs_tensorindextype_list>;
807 :
808 : Tensor<typename RhsTE::type, typename lhs_tensor_symm_and_indices::symmetry,
809 : typename lhs_tensor_symm_and_indices::tensorindextype_list>
810 : lhs_tensor{};
811 :
812 : evaluate<LhsTensorIndices...>(make_not_null(&lhs_tensor),
813 : rhs_tensorexpression);
814 : return lhs_tensor;
815 : }
816 :
817 : /*!
818 : * \ingroup TensorExpressionsGroup
819 : * \brief If the LHS tensor is used in the RHS expression, this should be used
820 : * to assign a LHS tensor to the result of a RHS tensor expression that contains
821 : * it
822 : *
823 : * \details See documentation for `tenex::evaluate` for basic functionality.
824 : *
825 : * `tenex::update` differs from `tenex::evaluate` in that `tenex::update` should
826 : * be used when some LHS `Tensor` has been partially computed, and now we would
827 : * like to update it with a RHS expression that contains it. In other words,
828 : * this should be used when we would like to emulate assignment operations like
829 : * `LHS +=`, `LHS -=`, `LHS *=`, etc.
830 : *
831 : * One important difference to note with `tenex::update` is that it cannot split
832 : * up the RHS expression and evaluate subtrees, while `tenex::evaluate` can (see
833 : * `TensorExpression` documentation). From benchmarking, it was found that the
834 : * runtime of `DataVector` expressions scales poorly as we increase the number
835 : * of operations. For this reason, when the data type held by the tensors in the
836 : * expression is `DataVector`, it's best to avoid passing RHS expressions with a
837 : * large number of operations (e.g. an inner product that sums over many terms).
838 : *
839 : * ### Example usage
840 : * In implementing a large equation with many operations, we can manually break
841 : * up the equation and evaluate different subexpressions at a time by making one
842 : * initial call to `tenex::evaluate` followed by any number of calls to
843 : * `tenex::update` that use the LHS tensor in the RHS expression and will
844 : * compute the rest of the equation:
845 : *
846 : * \snippet Test_MixedOperations.cpp use_update
847 : *
848 : * \note `LhsTensorIndices` must be passed by reference because non-type
849 : * template parameters cannot be class types until C++20.
850 : *
851 : * @tparam LhsTensorIndices the TensorIndexs of the Tensor on the LHS of the
852 : * tensor expression, e.g. `ti_a`, `ti_b`, `ti_c`
853 : * @param lhs_tensor pointer to the resultant LHS Tensor to fill
854 : * @param rhs_tensorexpression the RHS TensorExpression to be evaluated
855 : */
856 : template <auto&... LhsTensorIndices, typename LhsDataType, typename RhsDataType,
857 : typename LhsSymmetry, typename LhsIndexList, typename Derived,
858 : typename RhsSymmetry, typename RhsIndexList,
859 : typename... RhsTensorIndices>
860 1 : void update(
861 : const gsl::not_null<Tensor<LhsDataType, LhsSymmetry, LhsIndexList>*>
862 : lhs_tensor,
863 : const TensorExpression<Derived, RhsDataType, RhsSymmetry, RhsIndexList,
864 : tmpl::list<RhsTensorIndices...>>&
865 : rhs_tensorexpression) {
866 : using lhs_tensorindex_list =
867 : tmpl::list<std::decay_t<decltype(LhsTensorIndices)>...>;
868 : // Assert that each instance of the LHS tensor in the RHS tensor expression
869 : // uses the same generic index order that the LHS uses
870 : (~rhs_tensorexpression)
871 : .template assert_lhs_tensorindices_same_in_rhs<lhs_tensorindex_list>(
872 : lhs_tensor);
873 :
874 : detail::evaluate_impl<false, std::decay_t<decltype(LhsTensorIndices)>...>(
875 : lhs_tensor, rhs_tensorexpression,
876 : std::make_index_sequence<sizeof...(LhsTensorIndices)>{});
877 : }
878 : } // namespace tenex
|
