Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : ///\file
       5             : /// Defines an inner product for the linear solver
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <array>
      10             : 
      11             : #include "DataStructures/Variables.hpp"
      12             : #include "Utilities/Blas.hpp"
      13             : #include "Utilities/EqualWithinRoundoff.hpp"
      14             : #include "Utilities/ErrorHandling/Assert.hpp"
      15             : #include "Utilities/ForceInline.hpp"
      16             : 
      17             : namespace LinearSolver {
      18             : 
      19             : /// \ingroup LinearSolverGroup
      20             : /// Implementations of LinearSolver::inner_product.
      21           1 : namespace InnerProductImpls {
      22             : 
      23             : /// The inner product between any types that have a `dot` product
      24             : template <typename Lhs, typename Rhs>
      25           1 : struct InnerProductImpl {
      26           0 :   static auto apply(const Lhs& lhs, const Rhs& rhs) {
      27             :     return dot(conj(lhs), rhs);
      28             :   }
      29             : };
      30             : 
      31             : /// The inner product between `Variables`
      32             : template <typename LhsTagsList, typename RhsTagsList>
      33           1 : struct InnerProductImpl<Variables<LhsTagsList>, Variables<RhsTagsList>> {
      34           0 :   using ValueType = typename Variables<LhsTagsList>::value_type;
      35           0 :   static ValueType apply(const Variables<LhsTagsList>& lhs,
      36             :                          const Variables<RhsTagsList>& rhs) {
      37             :     const auto size = lhs.size();
      38             :     ASSERT(size == rhs.size(),
      39             :            "The Variables must be of the same size to take an inner product");
      40             :     if constexpr (std::is_same_v<ValueType, std::complex<double>>) {
      41             :       return zdotc_(size, lhs.data(), 1, rhs.data(), 1);
      42             :     } else {
      43             :       return ddot_(size, lhs.data(), 1, rhs.data(), 1);
      44             :     }
      45             :   }
      46             : };
      47             : 
      48             : }  // namespace InnerProductImpls
      49             : 
      50             : /*!
      51             :  * \ingroup LinearSolverGroup
      52             :  * \brief The local part of the Euclidean inner product on the vector space
      53             :  * w.r.t. which the addition and scalar multiplication of both `Lhs` and `Rhs`
      54             :  * is defined.
      55             :  *
      56             :  * \details The linear solver works under the following assumptions:
      57             :  * - The data represented by \p lhs and \p rhs can each be interpreted as the
      58             :  * local chunk of a vector of the same vector space \f$V\f$. _Local_ means there
      59             :  * are vectors \f$q, p\in V\f$ such that \p lhs and \p rhs represent the
      60             :  * components of these vectors w.r.t. a subset \f$B_i\f$ of a basis
      61             :  * \f$B\subset V\f$.
      62             :  * - The `*` and `+` operators of `Lhs` and `Rhs` implement the scalar
      63             :  * multiplication and addition in the vector space _locally_, i.e. restricted to
      64             :  * \f$B_i\f$ in the above sense.
      65             :  * - The inner product is the local part \f$\langle p,q\rangle|_{B_i}\f$ of the
      66             :  * standard Euclidean dot product in the vector space so that globally it is
      67             :  * \f$\langle p,q\rangle=\sum_{i}\langle p,q\rangle|_{B_i}\f$ for
      68             :  * \f$B=\mathop{\dot{\bigcup}}_i B_i\f$.
      69             :  *
      70             :  * In practice this means that the full vectors \f$p\f$ and \f$q\f$ can be
      71             :  * distributed on many elements, where each only holds local chunks \p lhs and
      72             :  * \p rhs of the components. Scalar multiplication and addition can be performed
      73             :  * locally as expected, but computing the full inner product requires a global
      74             :  * reduction over all elements that sums their local `inner_product`s.
      75             :  */
      76             : template <typename Lhs, typename Rhs>
      77           1 : SPECTRE_ALWAYS_INLINE auto inner_product(const Lhs& lhs, const Rhs& rhs) {
      78             :   return InnerProductImpls::InnerProductImpl<Lhs, Rhs>::apply(lhs, rhs);
      79             : }
      80             : 
      81             : /*!
      82             :  * \ingroup LinearSolverGroup
      83             :  * \brief The local part of the Euclidean inner product of a vector with itself.
      84             :  *
      85             :  * \see LinearSolver::inner_product for details on the inner product. This
      86             :  * function just returns the inner product of a vector with itself, and ensures
      87             :  * that the result is real.
      88             :  */
      89             : template <typename T>
      90           1 : SPECTRE_ALWAYS_INLINE double magnitude_square(const T& vector) {
      91             :   auto result =
      92             :       InnerProductImpls::InnerProductImpl<T, T>::apply(vector, vector);
      93             :   if constexpr (std::is_same_v<std::decay_t<decltype(result)>,
      94             :                                std::complex<double>>) {
      95             :     ASSERT(equal_within_roundoff(imag(result), 0.0),
      96             :            "The magnitude squared is not real. The imaginary part is: "
      97             :                << imag(result));
      98             :     return real(result);
      99             :   } else {
     100             :     return result;
     101             :   }
     102             : }
     103             : 
     104             : }  // namespace LinearSolver