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              : /// @}