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/TensorExpression.hpp"
20 : #include "DataStructures/Tensor/Expressions/TensorIndex.hpp"
21 : #include "DataStructures/Tensor/Expressions/TensorIndexTransformation.hpp"
22 : #include "DataStructures/Tensor/Expressions/TimeIndex.hpp"
23 : #include "DataStructures/Tensor/Structure.hpp"
24 : #include "DataStructures/Tensor/Tensor.hpp"
25 : #include "Utilities/ContainerHelpers.hpp"
26 : #include "Utilities/ErrorHandling/Assert.hpp"
27 : #include "Utilities/Gsl.hpp"
28 : #include "Utilities/Requires.hpp"
29 : #include "Utilities/TMPL.hpp"
30 :
31 : namespace tenex {
32 : namespace detail {
33 : template <size_t NumIndices>
34 : constexpr bool contains_indices_to_contract(
35 : const std::array<size_t, NumIndices>& tensorindices) {
36 : if constexpr (NumIndices < 2) {
37 : return false;
38 : } else {
39 : for (size_t i = 0; i < NumIndices - 1; i++) {
40 : for (size_t j = i + 1; j < NumIndices; j++) {
41 : const size_t current_tensorindex = gsl::at(tensorindices, i);
42 : // Concrete time indices are not contracted
43 : if ((not is_time_index_value(current_tensorindex)) and
44 : current_tensorindex == get_tensorindex_value_with_opposite_valence(
45 : gsl::at(tensorindices, j))) {
46 : return true;
47 : }
48 : }
49 : }
50 : return false;
51 : }
52 : }
53 :
54 : /// \brief Given the list of the positions of the LHS tensor's spacetime indices
55 : /// where a generic spatial index is used and the list of positions where a
56 : /// concrete time index is used, determine whether or not the component at the
57 : /// given LHS multi-index should be computed
58 : ///
59 : /// \details
60 : /// Not all of the LHS tensor's components may need to be computed. Cases when
61 : /// the component at a LHS multi-index should not be not evaluated:
62 : /// - If a generic spatial index is used for a spacetime index on the LHS,
63 : /// the components for which that index's concrete index is the time index
64 : /// should not be computed
65 : /// - If a concrete time index is used for a spacetime index on the LHS, the
66 : /// components for which that index's concrete index is a spatial index should
67 : /// not be computed
68 : ///
69 : /// \param lhs_multi_index the multi-index of the LHS tensor to check
70 : /// \param lhs_spatial_spacetime_index_positions the positions of the LHS
71 : /// tensor's spacetime indices where a generic spatial index is used
72 : /// \param lhs_time_index_positions the positions of the LHS tensor's spacetime
73 : /// indices where a concrete time index is used
74 : /// \return Whether or not `lhs_multi_index` is a multi-index of a component of
75 : /// the LHS tensor that should be computed
76 : template <size_t NumLhsIndices, size_t NumLhsSpatialSpacetimeIndices,
77 : size_t NumLhsConcreteTimeIndices>
78 : constexpr bool is_evaluated_lhs_multi_index(
79 : const std::array<size_t, NumLhsIndices>& lhs_multi_index,
80 : const std::array<size_t, NumLhsSpatialSpacetimeIndices>&
81 : lhs_spatial_spacetime_index_positions,
82 : const std::array<size_t, NumLhsConcreteTimeIndices>&
83 : lhs_time_index_positions) {
84 : for (size_t i = 0; i < lhs_spatial_spacetime_index_positions.size(); i++) {
85 : if (gsl::at(lhs_multi_index,
86 : gsl::at(lhs_spatial_spacetime_index_positions, i)) == 0) {
87 : return false;
88 : }
89 : }
90 : for (size_t i = 0; i < lhs_time_index_positions.size(); i++) {
91 : if (gsl::at(lhs_multi_index, gsl::at(lhs_time_index_positions, i)) != 0) {
92 : return false;
93 : }
94 : }
95 : return true;
96 : }
97 :
98 : template <typename SymmList>
99 : struct CheckNoLhsAntiSymmetries;
100 :
101 : template <template <typename...> class SymmList, typename... Symm>
102 : struct CheckNoLhsAntiSymmetries<SymmList<Symm...>> {
103 : static constexpr bool value = (... and (Symm::value > 0));
104 : };
105 :
106 : /*!
107 : * \ingroup TensorExpressionsGroup
108 : * \brief Evaluate subtrees of the RHS expression or the RHS expression as a
109 : * whole and assign the result to the LHS tensor
110 : *
111 : * \details This is for internal use only and should never be directly called.
112 : * See `tenex::evaluate` and use it, instead.
113 : *
114 : * `EvaluateSubtrees` controls whether we wish to evaluate RHS subtrees or the
115 : * entire RHS expression as one expression. See`TensorExpression` documentation
116 : * on equation splitting for more details on what this means.
117 : *
118 : * If `EvaluateSubtrees == false`, then it's safe if the LHS tensor is used in
119 : * the RHS expression, so long as the generic index orders are the same. This
120 : * means that the callee of this function needs to first verify this is true
121 : * before calling this function. Under these conditions, this is a safe
122 : * operation because the implementation modifies each LHS component once and
123 : * does not revisit and access any LHS components after they've been updated.
124 : * For example, say we do `tenex::evaluate<ti_a, ti_b>(make_not_null(&L),
125 : * 5.0 * L(ti_a, ti_b));`. This function will first compute the RHS for some
126 : * concrete LHS, e.g. \f$L_{00}\f$. To compute this, it accesses \f$L_{00}\f$
127 : * in the RHS tree, multiplies it by `5.0`, then updates \f$L_{00}\f$ to be the
128 : * result of this multiplication. Next, it might compute \f$L_{01}\f$, where
129 : * only \f$L_{01}\f$ is accessed, and which hasn't yet been modified. Then the
130 : * next component is computed and updated, and so forth. These steps are
131 : * performed once for each unique LHS index. Therefore, it is important to note
132 : * that this kind of operation being safe to perform is
133 : * implementation-dependent. Specifically, the safety of the operation depends
134 : * on the order of LHS component access and assignment.
135 : *
136 : * \note `LhsTensorIndices` must be passed by reference because non-type
137 : * template parameters cannot be class types until C++20.
138 : *
139 : * @tparam EvaluateSubtrees whether or not to evaluate subtrees of RHS
140 : * expression
141 : * @tparam LhsTensorIndices the `TensorIndex`s of the `Tensor` on the LHS of the
142 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
143 : * @param lhs_tensor pointer to the resultant LHS `Tensor` to fill
144 : * @param rhs_tensorexpression the RHS TensorExpression to be evaluated
145 : */
146 : template <bool EvaluateSubtrees, typename... LhsTensorIndices,
147 : typename LhsDataType, typename LhsSymmetry, typename LhsIndexList,
148 : typename Derived, typename RhsDataType, typename RhsSymmetry,
149 : typename RhsIndexList, typename... RhsTensorIndices>
150 : void evaluate_impl(
151 : const gsl::not_null<Tensor<LhsDataType, LhsSymmetry, LhsIndexList>*>
152 : lhs_tensor,
153 : const TensorExpression<Derived, RhsDataType, RhsSymmetry, RhsIndexList,
154 : tmpl::list<RhsTensorIndices...>>&
155 : rhs_tensorexpression) {
156 : constexpr size_t num_lhs_indices = sizeof...(LhsTensorIndices);
157 : constexpr size_t num_rhs_indices = sizeof...(RhsTensorIndices);
158 :
159 : using lhs_tensorindex_list = tmpl::list<LhsTensorIndices...>;
160 : using rhs_tensorindex_list = tmpl::list<RhsTensorIndices...>;
161 :
162 : using lhs_tensor_type = typename std::decay_t<decltype(*lhs_tensor)>;
163 :
164 : static_assert(is_supported_tensor_datatype_v<LhsDataType> and
165 : is_supported_tensor_datatype_v<RhsDataType>,
166 : "TensorExpressions currently only support Tensors whose data "
167 : "type is double, std::complex<double>, DataVector, or "
168 : "ComplexDataVector. It is possible to add support for other "
169 : "data types that are supported by Tensor.");
170 : static_assert(
171 : is_assignable_v<LhsDataType, RhsDataType>,
172 : "Assignment of the LHS Tensor's data type to the RHS TensorExpression's "
173 : "data type is not supported. This happens from doing something like e.g. "
174 : "trying to assign a Tensor<double> to a Tensor<DataVector> or a "
175 : "Tensor<DataVector> to a Tensor<ComplexDataVector>.");
176 : // `Symmetry` currently prevents this because antisymmetries are not currently
177 : // supported for `Tensor`s. This check is repeated here because if
178 : // antisymmetries are later supported for `Tensor`, using antisymmetries in
179 : // `TensorExpression`s will not automatically work. The implementations of the
180 : // derived `TensorExpression` types assume no antisymmetries (assume positive
181 : // `Symmetry` values), so support for antisymmetries in `TensorExpression`s
182 : // will still need to be implemented.
183 : static_assert(CheckNoLhsAntiSymmetries<LhsSymmetry>::value,
184 : "Anti-symmetric Tensors are not currently supported by "
185 : "TensorExpressions.");
186 : static_assert(
187 : tmpl::equal_members<
188 : typename remove_time_indices<lhs_tensorindex_list>::type,
189 : typename remove_time_indices<rhs_tensorindex_list>::type>::value,
190 : "The generic indices on the LHS of a tensor equation (that is, the "
191 : "template parameters specified in evaluate<...>) must match the generic "
192 : "indices of the RHS TensorExpression. This error occurs as a result of a "
193 : "call like evaluate<ti::a, ti::b>(R(ti::A, ti::b) * S(ti::a, ti::c)), "
194 : "where the generic indices of the evaluated RHS expression are ti::b and "
195 : "ti::c, but the generic indices provided for the LHS are ti::a and "
196 : "ti::b.");
197 : static_assert(
198 : tensorindex_list_is_valid<lhs_tensorindex_list>::value,
199 : "Cannot assign a tensor expression to a LHS tensor with a repeated "
200 : "generic index, e.g. evaluate<ti::a, ti::a>. (Note that the concrete "
201 : "time indices (ti::T and ti::t) can be repeated.)");
202 : static_assert(
203 : not contains_indices_to_contract<num_lhs_indices>(
204 : {{LhsTensorIndices::value...}}),
205 : "Cannot assign a tensor expression to a LHS tensor with generic "
206 : "indices that would be contracted, e.g. evaluate<ti::A, ti::a>.");
207 : // `IndexPropertyCheck` does also check that valence (Up/Lo) of indices that
208 : // correspond in the RHS and LHS tensors are equal, but the assertion message
209 : // below does not mention this because a mismatch in valence should have been
210 : // caught due to the combination of (i) the Tensor::operator() assertion
211 : // checking that generic indices' valences match the tensor's indices'
212 : // valences and (ii) the above assertion that RHS and LHS generic indices
213 : // match
214 : static_assert(
215 : IndexPropertyCheck<LhsIndexList, RhsIndexList, lhs_tensorindex_list,
216 : rhs_tensorindex_list>::value,
217 : "At least one index of the tensor evaluated from the RHS expression "
218 : "cannot be evaluated to its corresponding index in the LHS tensor. This "
219 : "is due to a difference in number of spatial dimensions or Frame type "
220 : "between the index on the RHS and LHS. "
221 : "e.g. evaluate<ti::a, ti::b>(L, R(ti::b, ti::a));, where R's first "
222 : "index has 2 spatial dimensions but L's second index has 3 spatial "
223 : "dimensions. Check RHS and LHS indices that use the same generic index.");
224 : static_assert(Derived::height_relative_to_closest_tensor_leaf_in_subtree <
225 : std::numeric_limits<size_t>::max(),
226 : "Either no Tensors were found in the RHS TensorExpression or "
227 : "the depth of the tree exceeded the maximum size_t value (very "
228 : "unlikely). If there is indeed a Tensor in the RHS expression "
229 : "and assuming the tree's height is not actually the maximum "
230 : "size_t value, then there is a flaw in the logic for computing "
231 : "the derived TensorExpression types' member, "
232 : "height_relative_to_closest_tensor_leaf_in_subtree.");
233 :
234 : if constexpr (EvaluateSubtrees) {
235 : // Make sure the LHS tensor doesn't also appear in the RHS tensor expression
236 : (~rhs_tensorexpression).assert_lhs_tensor_not_in_rhs_expression(lhs_tensor);
237 : // If the LHS data type is a vector, size the LHS tensor components if their
238 : // size does not match the size from a `Tensor` in the RHS expression
239 : if constexpr (is_derived_of_vector_impl_v<LhsDataType>) {
240 : const size_t rhs_component_size =
241 : (~rhs_tensorexpression).get_rhs_tensor_component_size();
242 : if (rhs_component_size != (*lhs_tensor)[0].size()) {
243 : for (auto& lhs_component : *lhs_tensor) {
244 : lhs_component = LhsDataType(rhs_component_size);
245 : }
246 : }
247 : }
248 : }
249 :
250 : constexpr std::array<size_t, num_rhs_indices> index_transformation =
251 : compute_tensorindex_transformation<num_lhs_indices, num_rhs_indices>(
252 : {{LhsTensorIndices::value...}}, {{RhsTensorIndices::value...}});
253 :
254 : // positions of indices in LHS tensor where generic spatial indices are used
255 : // for spacetime indices
256 : constexpr auto lhs_spatial_spacetime_index_positions =
257 : get_spatial_spacetime_index_positions<LhsIndexList,
258 : lhs_tensorindex_list>();
259 : // positions of indices in RHS tensor where generic spatial indices are used
260 : // for spacetime indices
261 : constexpr auto rhs_spatial_spacetime_index_positions =
262 : get_spatial_spacetime_index_positions<RhsIndexList,
263 : rhs_tensorindex_list>();
264 :
265 : // positions of indices in LHS tensor where concrete time indices are used
266 : constexpr auto lhs_time_index_positions =
267 : get_time_index_positions<lhs_tensorindex_list>();
268 :
269 : using rhs_expression_type =
270 : typename std::decay_t<decltype(~rhs_tensorexpression)>;
271 :
272 : for (size_t i = 0; i < lhs_tensor_type::size(); i++) {
273 : auto lhs_multi_index =
274 : lhs_tensor_type::structure::get_canonical_tensor_index(i);
275 : if (is_evaluated_lhs_multi_index(lhs_multi_index,
276 : lhs_spatial_spacetime_index_positions,
277 : lhs_time_index_positions)) {
278 : for (size_t j = 0; j < lhs_spatial_spacetime_index_positions.size();
279 : j++) {
280 : gsl::at(lhs_multi_index,
281 : gsl::at(lhs_spatial_spacetime_index_positions, j)) -= 1;
282 : }
283 : auto rhs_multi_index =
284 : transform_multi_index(lhs_multi_index, index_transformation);
285 : for (size_t j = 0; j < rhs_spatial_spacetime_index_positions.size();
286 : j++) {
287 : gsl::at(rhs_multi_index,
288 : gsl::at(rhs_spatial_spacetime_index_positions, j)) += 1;
289 : }
290 :
291 : // The expression will either be evaluated as one whole expression
292 : // or it will be split up into subtrees that are evaluated one at a time.
293 : // See the section on splitting in the documentation for the
294 : // `TensorExpression` class to understand the logic and terminology used
295 : // in this control flow below.
296 : if constexpr (EvaluateSubtrees) {
297 : // the expression is split up, so evaluate subtrees at splits
298 : (~rhs_tensorexpression)
299 : .evaluate_primary_subtree((*lhs_tensor)[i], rhs_multi_index);
300 : if constexpr (not rhs_expression_type::is_primary_start) {
301 : // the root expression type is not the starting point of a leg, so it
302 : // has not yet been evaluated, so now we evaluate this last leg of the
303 : // expression at the root of the tree
304 : (*lhs_tensor)[i] =
305 : (~rhs_tensorexpression)
306 : .get_primary((*lhs_tensor)[i], rhs_multi_index);
307 : }
308 : } else {
309 : // the expression is not split up, so evaluate full expression
310 : (*lhs_tensor)[i] = (~rhs_tensorexpression).get(rhs_multi_index);
311 : }
312 : }
313 : }
314 : }
315 :
316 : /*!
317 : * \ingroup TensorExpressionsGroup
318 : * \brief Assign a value to components of the LHS tensor
319 : *
320 : * \details This is for internal use only and should never be directly called.
321 : * See `tenex::evaluate` and use it, instead.
322 : *
323 : * \note `LhsTensorIndices` must be passed by reference because non-type
324 : * template parameters cannot be class types until C++20.
325 : *
326 : * @tparam LhsTensorIndices the `TensorIndex`s of the `Tensor` on the LHS of the
327 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
328 : * @param lhs_tensor pointer to the resultant LHS `Tensor` to fill
329 : * @param rhs_value the RHS value to assigned
330 : */
331 : template <typename... LhsTensorIndices, typename X, typename LhsSymmetry,
332 : typename LhsIndexList, typename NumberType>
333 : void evaluate_impl(
334 : const gsl::not_null<Tensor<X, LhsSymmetry, LhsIndexList>*> lhs_tensor,
335 : const NumberType& rhs_value) {
336 : using lhs_tensor_type = typename std::decay_t<decltype(*lhs_tensor)>;
337 : constexpr size_t num_lhs_indices = sizeof...(LhsTensorIndices);
338 : using lhs_tensorindex_list = tmpl::list<LhsTensorIndices...>;
339 :
340 : static_assert(is_supported_tensor_datatype_v<X> and
341 : "TensorExpressions currently only support Tensors whose data "
342 : "type is double, std::complex<double>, DataVector, or "
343 : "ComplexDataVector. It is possible to add support for other "
344 : "data types that are supported by Tensor.");
345 : static_assert(
346 : is_assignable_v<X, NumberType>,
347 : "Assignment of the LHS Tensor's data type to the RHS number's data type "
348 : "is not supported within TensorExpressions. This happens from doing "
349 : "something like e.g. trying to assign a double to a DataVector or a "
350 : "DataVector to a ComplexDataVector.");
351 : // `Symmetry` currently prevents this because antisymmetries are not currently
352 : // supported for `Tensor`s. This check is repeated here because if
353 : // antisymmetries are later supported for `Tensor`, using antisymmetries in
354 : // `TensorExpression`s will not automatically work. The implementations of the
355 : // derived `TensorExpression` types assume no antisymmetries (assume positive
356 : // `Symmetry` values), so support for antisymmetries in `TensorExpression`s
357 : // will still need to be implemented.
358 : static_assert(CheckNoLhsAntiSymmetries<LhsSymmetry>::value,
359 : "Anti-symmetric Tensors are not currently supported by "
360 : "TensorExpressions.");
361 : static_assert(
362 : tensorindex_list_is_valid<lhs_tensorindex_list>::value,
363 : "Cannot assign a tensor expression to a LHS tensor with a repeated "
364 : "generic index, e.g. evaluate<ti::a, ti::a>. (Note that the concrete "
365 : "time indices (ti::T and ti::t) can be repeated.)");
366 : static_assert(
367 : not contains_indices_to_contract<num_lhs_indices>(
368 : {{LhsTensorIndices::value...}}),
369 : "Cannot assign a tensor expression to a LHS tensor with generic "
370 : "indices that would be contracted, e.g. evaluate<ti::A, ti::a>.");
371 :
372 : if constexpr (is_derived_of_vector_impl_v<X>) {
373 : ASSERT(get_size((*lhs_tensor)[0]) > 0,
374 : "Tensors with vector components must be sized before calling "
375 : "\ntenex::evaluate<...>("
376 : "\n\tgsl::not_null<Tensor<VectorType, ...>*>, number).");
377 : }
378 :
379 : // positions of indices in LHS tensor where generic spatial indices are used
380 : // for spacetime indices
381 : constexpr auto lhs_spatial_spacetime_index_positions =
382 : get_spatial_spacetime_index_positions<LhsIndexList,
383 : lhs_tensorindex_list>();
384 :
385 : // positions of indices in LHS tensor where concrete time indices are used
386 : constexpr auto lhs_time_index_positions =
387 : get_time_index_positions<lhs_tensorindex_list>();
388 :
389 : for (size_t i = 0; i < lhs_tensor_type::size(); i++) {
390 : auto lhs_multi_index =
391 : lhs_tensor_type::structure::get_canonical_tensor_index(i);
392 : if (is_evaluated_lhs_multi_index(lhs_multi_index,
393 : lhs_spatial_spacetime_index_positions,
394 : lhs_time_index_positions)) {
395 : (*lhs_tensor)[i] = rhs_value;
396 : }
397 : }
398 : }
399 : } // namespace detail
400 :
401 : /*!
402 : * \ingroup TensorExpressionsGroup
403 : * \brief Assign the result of a RHS tensor expression to a tensor with the LHS
404 : * index order set in the template parameters
405 : *
406 : * \details Uses the right hand side (RHS) TensorExpression's index ordering
407 : * (`RhsTE::args_list`) and the desired left hand side (LHS) tensor's index
408 : * ordering (`LhsTensorIndices`) to fill the provided LHS Tensor with that LHS
409 : * index ordering. This can carry out the evaluation of a RHS tensor expression
410 : * to a LHS tensor with the same index ordering, such as \f$L_{ab} = R_{ab}\f$,
411 : * or different ordering, such as \f$L_{ba} = R_{ab}\f$.
412 : *
413 : * The symmetry of the provided LHS Tensor need not match the symmetry
414 : * determined from evaluating the RHS TensorExpression according to its order of
415 : * operations. This allows one to specify LHS symmetries (via `lhs_tensor`) that
416 : * may not be preserved by the RHS expression's order of operations, which
417 : * depends on how the expression is written and implemented.
418 : *
419 : * The LHS `Tensor` cannot be part of the RHS expression, e.g.
420 : * `evaluate(make_not_null(&L), L() + R());`, because the LHS `Tensor` will
421 : * generally not be computed correctly when the RHS `TensorExpression` is split
422 : * up and the LHS tensor components are computed by accumulating the result of
423 : * subtrees (see the section on splitting in the documentation for the
424 : * `TensorExpression` class). If you need to use the LHS `Tensor` on the RHS,
425 : * use `tenex::update` instead.
426 : *
427 : * ### Example usage
428 : * Given `Tensor`s `R`, `S`, `T`, `G`, and `H`, we can compute the LHS tensor
429 : * \f$L\f$ in the equation \f$L_{a} = R_{ab} S^{b} + G_{a} - H_{ba}{}^{b} T\f$
430 : * by doing:
431 : *
432 : * \snippet Test_MixedOperations.cpp use_evaluate_with_result_as_arg
433 : *
434 : * \note `LhsTensorIndices` must be passed by reference because non-type
435 : * template parameters cannot be class types until C++20.
436 : *
437 : * @tparam LhsTensorIndices the `TensorIndex`s of the `Tensor` on the LHS of the
438 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
439 : * @param lhs_tensor pointer to the resultant LHS `Tensor` to fill
440 : * @param rhs_tensorexpression the RHS TensorExpression to be evaluated
441 : */
442 : template <auto&... LhsTensorIndices, typename LhsDataType, typename LhsSymmetry,
443 : typename LhsIndexList, typename Derived, typename RhsDataType,
444 : typename RhsSymmetry, typename RhsIndexList,
445 : typename... RhsTensorIndices>
446 1 : void evaluate(
447 : const gsl::not_null<Tensor<LhsDataType, LhsSymmetry, LhsIndexList>*>
448 : lhs_tensor,
449 : const TensorExpression<Derived, RhsDataType, RhsSymmetry, RhsIndexList,
450 : tmpl::list<RhsTensorIndices...>>&
451 : rhs_tensorexpression) {
452 : using rhs_expression_type =
453 : typename std::decay_t<decltype(~rhs_tensorexpression)>;
454 : constexpr bool evaluate_subtrees =
455 : rhs_expression_type::primary_subtree_contains_primary_start;
456 : detail::evaluate_impl<evaluate_subtrees,
457 : std::decay_t<decltype(LhsTensorIndices)>...>(
458 : lhs_tensor, rhs_tensorexpression);
459 : }
460 :
461 : /// @{
462 : /*!
463 : * \ingroup TensorExpressionsGroup
464 : * \brief Assign a number to components of a tensor with the LHS index order
465 : * set in the template parameters
466 : *
467 : * \details
468 : * Example usage:
469 : * \snippet Test_MixedOperations.cpp assign_double_to_index_subsets
470 : *
471 : * \note The components of the LHS `Tensor` passed in must already be sized
472 : * because there is no way to infer component size from the RHS
473 : *
474 : * \note `LhsTensorIndices` must be passed by reference because non-type
475 : * template parameters cannot be class types until C++20.
476 : *
477 : * @tparam LhsTensorIndices the `TensorIndex`s of the `Tensor` on the LHS of the
478 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
479 : * @param lhs_tensor pointer to the resultant LHS `Tensor` to fill
480 : * @param rhs_value the RHS value to assign
481 : */
482 : template <auto&... LhsTensorIndices, typename X, typename LhsSymmetry,
483 : typename LhsIndexList, typename N,
484 : Requires<std::is_arithmetic_v<N>> = nullptr>
485 1 : void evaluate(
486 : const gsl::not_null<Tensor<X, LhsSymmetry, LhsIndexList>*> lhs_tensor,
487 : const N rhs_value) {
488 : detail::evaluate_impl<std::decay_t<decltype(LhsTensorIndices)>...>(lhs_tensor,
489 : rhs_value);
490 : }
491 : template <auto&... LhsTensorIndices, typename X, typename LhsSymmetry,
492 : typename LhsIndexList, typename N>
493 1 : void evaluate(
494 : const gsl::not_null<Tensor<X, LhsSymmetry, LhsIndexList>*> lhs_tensor,
495 : const std::complex<N>& rhs_value) {
496 : detail::evaluate_impl<std::decay_t<decltype(LhsTensorIndices)>...>(lhs_tensor,
497 : rhs_value);
498 : }
499 : /// @}
500 :
501 : /*!
502 : * \ingroup TensorExpressionsGroup
503 : * \brief Assign the result of a RHS tensor expression to a tensor with the LHS
504 : * index order set in the template parameters
505 : *
506 : * \details Uses the right hand side (RHS) TensorExpression's index ordering
507 : * (`RhsTE::args_list`) and the desired left hand side (LHS) tensor's index
508 : * ordering (`LhsTensorIndices`) to construct a LHS Tensor with that LHS index
509 : * ordering. This can carry out the evaluation of a RHS tensor expression to a
510 : * LHS tensor with the same index ordering, such as \f$L_{ab} = R_{ab}\f$, or
511 : * different ordering, such as \f$L_{ba} = R_{ab}\f$.
512 : *
513 : * The symmetry of the returned LHS Tensor depends on the order of operations in
514 : * the RHS TensorExpression, i.e. how the expression is written. If you would
515 : * like to specify the symmetry of the LHS Tensor instead of it being determined
516 : * by the order of operations in the RHS expression, please use the other
517 : * `tenex::evaluate` overload that takes an empty LHS Tensor as its first
518 : * argument.
519 : *
520 : * ### Example usage
521 : * Given `Tensor`s `R`, `S`, `T`, `G`, and `H`, we can compute the LHS tensor
522 : * \f$L\f$ in the equation \f$L_{a} = R_{ab} S^{b} + G_{a} - H_{ba}{}^{b} T\f$
523 : * by doing:
524 : *
525 : * \snippet Test_MixedOperations.cpp use_evaluate_to_return_result
526 : *
527 : * \parblock
528 : * \note If a generic spatial index is used for a spacetime index in the RHS
529 : * tensor, its corresponding index in the LHS tensor type will be a spatial
530 : * index with the same valence, frame, and number of spatial dimensions. If a
531 : * concrete time index is used for a spacetime index in the RHS tensor, the
532 : * index will not appear in the LHS tensor (i.e. there will NOT be a
533 : * corresponding LHS index where only the time index of that index has been
534 : * computed and its spatial indices are empty).
535 : * \endparblock
536 : *
537 : * \parblock
538 : * \note `LhsTensorIndices` must be passed by reference because non-type
539 : * template parameters cannot be class types until C++20.
540 : * \endparblock
541 : *
542 : * @tparam LhsTensorIndices the TensorIndexs of the Tensor on the LHS of the
543 : * tensor expression, e.g. `ti::a`, `ti::b`, `ti::c`
544 : * @param rhs_tensorexpression the RHS TensorExpression to be evaluated
545 : * @return the resultant LHS Tensor with index order specified by
546 : * LhsTensorIndices
547 : */
548 : template <auto&... LhsTensorIndices, typename RhsTE,
549 : Requires<std::is_base_of_v<Expression, RhsTE>> = nullptr>
550 1 : auto evaluate(const RhsTE& rhs_tensorexpression) {
551 : using lhs_tensorindex_list =
552 : tmpl::list<std::decay_t<decltype(LhsTensorIndices)>...>;
553 : using rhs_tensorindex_list = typename RhsTE::args_list;
554 : using rhs_symmetry = typename RhsTE::symmetry;
555 : using rhs_tensorindextype_list = typename RhsTE::index_list;
556 :
557 : // Stores (potentially reordered) symmetry and indices needed for constructing
558 : // the LHS tensor, with index order specified by LhsTensorIndices
559 : using lhs_tensor_symm_and_indices =
560 : LhsTensorSymmAndIndices<rhs_tensorindex_list, lhs_tensorindex_list,
561 : rhs_symmetry, rhs_tensorindextype_list>;
562 :
563 : Tensor<typename RhsTE::type, typename lhs_tensor_symm_and_indices::symmetry,
564 : typename lhs_tensor_symm_and_indices::tensorindextype_list>
565 : lhs_tensor{};
566 : evaluate<LhsTensorIndices...>(make_not_null(&lhs_tensor),
567 : rhs_tensorexpression);
568 : return lhs_tensor;
569 : }
570 :
571 : /*!
572 : * \ingroup TensorExpressionsGroup
573 : * \brief If the LHS tensor is used in the RHS expression, this should be used
574 : * to assign a LHS tensor to the result of a RHS tensor expression that contains
575 : * it
576 : *
577 : * \details See documentation for `tenex::evaluate` for basic functionality.
578 : *
579 : * `tenex::update` differs from `tenex::evaluate` in that `tenex::update` should
580 : * be used when some LHS `Tensor` has been partially computed, and now we would
581 : * like to update it with a RHS expression that contains it. In other words,
582 : * this should be used when we would like to emulate assignment operations like
583 : * `LHS +=`, `LHS -=`, `LHS *=`, etc.
584 : *
585 : * One important difference to note with `tenex::update` is that it cannot split
586 : * up the RHS expression and evaluate subtrees, while `tenex::evaluate` can (see
587 : * `TensorExpression` documentation). From benchmarking, it was found that the
588 : * runtime of `DataVector` expressions scales poorly as we increase the number
589 : * of operations. For this reason, when the data type held by the tensors in the
590 : * expression is `DataVector`, it's best to avoid passing RHS expressions with a
591 : * large number of operations (e.g. an inner product that sums over many terms).
592 : *
593 : * ### Example usage
594 : * In implementing a large equation with many operations, we can manually break
595 : * up the equation and evaluate different subexpressions at a time by making one
596 : * initial call to `tenex::evaluate` followed by any number of calls to
597 : * `tenex::update` that use the LHS tensor in the RHS expression and will
598 : * compute the rest of the equation:
599 : *
600 : * \snippet Test_MixedOperations.cpp use_update
601 : *
602 : * \note `LhsTensorIndices` must be passed by reference because non-type
603 : * template parameters cannot be class types until C++20.
604 : *
605 : * @tparam LhsTensorIndices the TensorIndexs of the Tensor on the LHS of the
606 : * tensor expression, e.g. `ti_a`, `ti_b`, `ti_c`
607 : * @param lhs_tensor pointer to the resultant LHS Tensor to fill
608 : * @param rhs_tensorexpression the RHS TensorExpression to be evaluated
609 : */
610 : template <auto&... LhsTensorIndices, typename LhsDataType, typename RhsDataType,
611 : typename LhsSymmetry, typename LhsIndexList, typename Derived,
612 : typename RhsSymmetry, typename RhsIndexList,
613 : typename... RhsTensorIndices>
614 1 : void update(
615 : const gsl::not_null<Tensor<LhsDataType, LhsSymmetry, LhsIndexList>*>
616 : lhs_tensor,
617 : const TensorExpression<Derived, RhsDataType, RhsSymmetry, RhsIndexList,
618 : tmpl::list<RhsTensorIndices...>>&
619 : rhs_tensorexpression) {
620 : using lhs_tensorindex_list =
621 : tmpl::list<std::decay_t<decltype(LhsTensorIndices)>...>;
622 : // Assert that each instance of the LHS tensor in the RHS tensor expression
623 : // uses the same generic index order that the LHS uses
624 : (~rhs_tensorexpression)
625 : .template assert_lhs_tensorindices_same_in_rhs<lhs_tensorindex_list>(
626 : lhs_tensor);
627 :
628 : detail::evaluate_impl<false, std::decay_t<decltype(LhsTensorIndices)>...>(
629 : lhs_tensor, rhs_tensorexpression);
630 : }
631 : } // namespace tenex
|