SpECTRE Documentation Coverage Report
Current view: top level - DataStructures/Tensor/Expressions - Evaluate.hpp Hit Total Coverage
Commit: d0fc80462417e83e5cddfa1b9901bb4a9b6af4d6 Lines: 6 6 100.0 %
Date: 2024-03-29 00:33:31
Legend: Lines: hit not hit

          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

Generated by: LCOV version 1.14