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, 1e-14, abs(result)),
      96              :            "The magnitude squared is not real: " << result);
      97              :     return real(result);
      98              :   } else {
      99              :     return result;
     100              :   }
     101              : }
     102              : 
     103              : }  // namespace LinearSolver