Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : 
       8             : #include "DataStructures/Tensor/Tensor.hpp"
       9             : #include "DataStructures/Variables.hpp"
      10             : #include "Utilities/Gsl.hpp"
      11             : #include "Utilities/TMPL.hpp"
      12             : 
      13             : /*!
      14             :  * \brief The `Tensor` type resulting from the outer product of two tensors
      15             :  *
      16             :  * \see outer_product
      17             :  */
      18             : template <typename DataType, typename SymmLhs, typename IndicesLhs,
      19             :           typename SymmRhs, typename IndicesRhs>
      20           1 : using OuterProductResultTensor = Tensor<
      21             :     DataType,
      22             :     tmpl::append<
      23             :         tmpl::transform<
      24             :             SymmLhs,
      25             :             tmpl::plus<tmpl::_1,
      26             :                        tmpl::fold<SymmRhs, tmpl::int32_t<0>,
      27             :                                   tmpl::max<tmpl::_state, tmpl::_element>>>>,
      28             :         SymmRhs>,
      29             :     tmpl::append<IndicesLhs, IndicesRhs>>;
      30             : 
      31             : /// @{
      32             : /*!
      33             :  * \ingroup TensorGroup
      34             :  * \brief The outer product (or tensor product) of two tensors
      35             :  *
      36             :  * \details
      37             :  * Computes $A_{i\ldots j\ldots} = B_{i\ldots} C_{j\ldots}$ for two tensors
      38             :  * $B_{i\ldots}$ and $C_{j\ldots}$. Both tensors can have arbitrary indices and
      39             :  * symmetries.
      40             :  */
      41             : template <typename DataTypeLhs, typename SymmLhs, typename IndicesLhs,
      42             :           typename DataTypeRhs, typename SymmRhs, typename IndicesRhs,
      43             :           typename DataTypeResult = decltype(blaze::evaluate(DataTypeLhs() *
      44             :                                                              DataTypeRhs()))>
      45           1 : void outer_product(
      46             :     const gsl::not_null<OuterProductResultTensor<
      47             :         DataTypeResult, SymmLhs, IndicesLhs, SymmRhs, IndicesRhs>*>
      48             :         result,
      49             :     const Tensor<DataTypeLhs, SymmLhs, IndicesLhs>& lhs,
      50             :     const Tensor<DataTypeRhs, SymmRhs, IndicesRhs>& rhs) {
      51             :   for (auto it_lhs = lhs.begin(); it_lhs != lhs.end(); ++it_lhs) {
      52             :     for (auto it_rhs = rhs.begin(); it_rhs != rhs.end(); ++it_rhs) {
      53             :       const auto result_indices = concatenate(lhs.get_tensor_index(it_lhs),
      54             :                                               rhs.get_tensor_index(it_rhs));
      55             :       result->get(result_indices) = *it_lhs * *it_rhs;
      56             :     }
      57             :   }
      58             : }
      59             : 
      60             : template <typename DataTypeLhs, typename SymmLhs, typename IndicesLhs,
      61             :           typename DataTypeRhs, typename SymmRhs, typename IndicesRhs,
      62             :           typename DataTypeResult = decltype(blaze::evaluate(DataTypeLhs() *
      63             :                                                              DataTypeRhs()))>
      64           1 : auto outer_product(const Tensor<DataTypeLhs, SymmLhs, IndicesLhs>& lhs,
      65             :                    const Tensor<DataTypeRhs, SymmRhs, IndicesRhs>& rhs)
      66             :     -> OuterProductResultTensor<DataTypeResult, SymmLhs, IndicesLhs, SymmRhs,
      67             :                                 IndicesRhs> {
      68             :   OuterProductResultTensor<DataTypeResult, SymmLhs, IndicesLhs, SymmRhs,
      69             :                            IndicesRhs>
      70             :       result{};
      71             :   ::outer_product(make_not_null(&result), lhs, rhs);
      72             :   return result;
      73             : }
      74             : /// @}