Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <ostream>
       9             : 
      10             : #include "DataStructures/Variables.hpp"
      11             : #include "Utilities/DereferenceWrapper.hpp"
      12             : #include "Utilities/ErrorHandling/Assert.hpp"
      13             : #include "Utilities/Gsl.hpp"
      14             : 
      15             : /// \cond
      16             : template <size_t Dim>
      17             : class Index;
      18             : // IWYU pragma: no_forward_declare Variables
      19             : /// \endcond
      20             : 
      21             : namespace apply_matrices_detail {
      22             : template <typename ElementType, size_t Dim, bool... DimensionIsIdentity>
      23             : struct Impl {
      24             :   template <typename MatrixType>
      25             :   static void apply(gsl::not_null<ElementType*> result,
      26             :                     const std::array<MatrixType, Dim>& matrices,
      27             :                     const ElementType* data, const Index<Dim>& extents,
      28             :                     size_t number_of_independent_components);
      29             : };
      30             : 
      31             : template <typename MatrixType, size_t Dim>
      32             : size_t result_size(const std::array<MatrixType, Dim>& matrices,
      33             :                    const Index<Dim>& extents) {
      34             :   size_t num_points_result = 1;
      35             :   for (size_t d = 0; d < Dim; ++d) {
      36             :     const size_t cols = dereference_wrapper(gsl::at(matrices, d)).columns();
      37             :     if (cols == 0) {
      38             :       // An empty matrix is treated as the identity.
      39             :       num_points_result *= extents[d];
      40             :     } else {
      41             :       ASSERT(cols == extents[d],
      42             :              "Matrix " << d << " has wrong number of columns: " << cols
      43             :                        << " (expected " << extents[d] << ")");
      44             :       num_points_result *= dereference_wrapper(gsl::at(matrices, d)).rows();
      45             :     }
      46             :   }
      47             :   return num_points_result;
      48             : }
      49             : }  // namespace apply_matrices_detail
      50             : 
      51             : /// @{
      52             : /// \ingroup NumericalAlgorithmsGroup
      53             : /// \brief Multiply by matrices in each dimension
      54             : ///
      55             : /// Multiplies each stripe in the first dimension of `u` by
      56             : /// `matrices[0]`, each stripe in the second dimension of `u` by
      57             : /// `matrices[1]`, and so on.  If any of the matrices are empty they
      58             : /// will be treated as the identity, but the matrix multiplications
      59             : /// will be skipped for increased efficiency.
      60             : ///
      61             : /// \note The element type stored in the vectors to be transformed may be either
      62             : /// `double` or `std::complex<double>`. The matrix, however, must be real. In
      63             : /// the case of acting on a vector of complex values, the matrix is treated as
      64             : /// having zero imaginary part. This is chosen for efficiency in all
      65             : /// use-cases for spectral matrix arithmetic so far encountered.
      66             : template <typename VariableTags, typename MatrixType, size_t Dim>
      67           1 : void apply_matrices(const gsl::not_null<Variables<VariableTags>*> result,
      68             :                     const std::array<MatrixType, Dim>& matrices,
      69             :                     const Variables<VariableTags>& u,
      70             :                     const Index<Dim>& extents) {
      71             :   ASSERT(u.number_of_grid_points() == extents.product(),
      72             :          "Mismatch between extents (" << extents.product()
      73             :                                       << ") and variables ("
      74             :                                       << u.number_of_grid_points() << ").");
      75             :   ASSERT(result->number_of_grid_points() ==
      76             :              apply_matrices_detail::result_size(matrices, extents),
      77             :          "result has wrong size.  Expected "
      78             :              << apply_matrices_detail::result_size(matrices, extents)
      79             :              << ", received " << result->number_of_grid_points());
      80             :   apply_matrices_detail::Impl<typename Variables<VariableTags>::value_type,
      81             :                               Dim>::apply(result->data(), matrices, u.data(),
      82             :                                           extents,
      83             :                                           u.number_of_independent_components);
      84             : }
      85             : 
      86             : template <typename VariableTags, typename MatrixType, size_t Dim>
      87           1 : Variables<VariableTags> apply_matrices(
      88             :     const std::array<MatrixType, Dim>& matrices,
      89             :     const Variables<VariableTags>& u, const Index<Dim>& extents) {
      90             :   Variables<VariableTags> result(
      91             :       apply_matrices_detail::result_size(matrices, extents));
      92             :   apply_matrices(make_not_null(&result), matrices, u, extents);
      93             :   return result;
      94             : }
      95             : 
      96             : // clang tidy mistakenly fails to identify this as a function definition
      97             : template <typename ResultType, typename MatrixType, typename VectorType,
      98             :           size_t Dim>
      99           1 : void apply_matrices(const gsl::not_null<ResultType*> result,  // NOLINT
     100             :                     const std::array<MatrixType, Dim>& matrices,
     101             :                     const VectorType& u, const Index<Dim>& extents) {
     102             :   const size_t number_of_independent_components = u.size() / extents.product();
     103             :   ASSERT(u.size() == number_of_independent_components * extents.product(),
     104             :          "The size of the vector u ("
     105             :              << u.size()
     106             :              << ") must be a multiple of the number of grid points ("
     107             :              << extents.product() << ").");
     108             :   ASSERT(result->size() ==
     109             :              number_of_independent_components *
     110             :                  apply_matrices_detail::result_size(matrices, extents),
     111             :          "result has wrong size.  Expected "
     112             :              << number_of_independent_components *
     113             :                     apply_matrices_detail::result_size(matrices, extents)
     114             :              << ", received " << result->size());
     115             :   apply_matrices_detail::Impl<typename VectorType::ElementType, Dim>::apply(
     116             :       result->data(), matrices, u.data(), extents,
     117             :       number_of_independent_components);
     118             : }
     119             : 
     120             : template <typename MatrixType, typename VectorType, size_t Dim>
     121           1 : VectorType apply_matrices(const std::array<MatrixType, Dim>& matrices,
     122             :                           const VectorType& u, const Index<Dim>& extents) {
     123             :   const size_t number_of_independent_components = u.size() / extents.product();
     124             :   VectorType result(number_of_independent_components *
     125             :                     apply_matrices_detail::result_size(matrices, extents));
     126             :   apply_matrices(make_not_null(&result), matrices, u, extents);
     127             :   return result;
     128             : }
     129             : 
     130             : template <typename ResultType, typename MatrixType, typename VectorType,
     131             :           size_t Dim>
     132           1 : ResultType apply_matrices(const std::array<MatrixType, Dim>& matrices,
     133             :                           const VectorType& u, const Index<Dim>& extents) {
     134             :   const size_t number_of_independent_components = u.size() / extents.product();
     135             :   ResultType result(number_of_independent_components *
     136             :                     apply_matrices_detail::result_size(matrices, extents));
     137             :   apply_matrices(make_not_null(&result), matrices, u, extents);
     138             :   return result;
     139             : }
     140             : /// @}