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/ForceInline.hpp" 14 : 15 : namespace LinearSolver { 16 : 17 : /// \ingroup LinearSolverGroup 18 : /// Implementations of LinearSolver::inner_product. 19 1 : namespace InnerProductImpls { 20 : 21 : /// The inner product between any types that have a `dot` product 22 : template <typename Lhs, typename Rhs> 23 1 : struct InnerProductImpl { 24 0 : static double apply(const Lhs& lhs, const Rhs& rhs) { return dot(lhs, rhs); } 25 : }; 26 : 27 : /// The inner product between `Variables` 28 : template <typename LhsTagsList, typename RhsTagsList> 29 1 : struct InnerProductImpl<Variables<LhsTagsList>, Variables<RhsTagsList>> { 30 0 : static double apply(const Variables<LhsTagsList>& lhs, 31 : const Variables<RhsTagsList>& rhs) { 32 : const auto size = lhs.size(); 33 : ASSERT(size == rhs.size(), 34 : "The Variables must be of the same size to take an inner product"); 35 : return ddot_(size, lhs.data(), 1, rhs.data(), 1); 36 : } 37 : }; 38 : 39 : } // namespace InnerProductImpls 40 : 41 : /*! 42 : * \ingroup LinearSolverGroup 43 : * \brief The local part of the Euclidean inner product on the vector space 44 : * w.r.t. which the addition and scalar multiplication of both `Lhs` and `Rhs` 45 : * is defined. 46 : * 47 : * \details The linear solver works under the following assumptions: 48 : * - The data represented by \p lhs and \p rhs can each be interpreted as the 49 : * local chunk of a vector of the same vector space \f$V\f$. _Local_ means there 50 : * are vectors \f$q, p\in V\f$ such that \p lhs and \p rhs represent the 51 : * components of these vectors w.r.t. a subset \f$B_i\f$ of a basis 52 : * \f$B\subset V\f$. 53 : * - The `*` and `+` operators of `Lhs` and `Rhs` implement the scalar 54 : * multiplication and addition in the vector space _locally_, i.e. restricted to 55 : * \f$B_i\f$ in the above sense. 56 : * - The inner product is the local part \f$\langle p,q\rangle|_{B_i}\f$ of the 57 : * standard Euclidean dot product in the vector space so that globally it is 58 : * \f$\langle p,q\rangle=\sum_{i}\langle p,q\rangle|_{B_i}\f$ for 59 : * \f$B=\mathop{\dot{\bigcup}}_i B_i\f$. 60 : * 61 : * In practice this means that the full vectors \f$p\f$ and \f$q\f$ can be 62 : * distributed on many elements, where each only holds local chunks \p lhs and 63 : * \p rhs of the components. Scalar multiplication and addition can be performed 64 : * locally as expected, but computing the full inner product requires a global 65 : * reduction over all elements that sums their local `inner_product`s. 66 : */ 67 : template <typename Lhs, typename Rhs> 68 1 : SPECTRE_ALWAYS_INLINE double inner_product(const Lhs& lhs, const Rhs& rhs) { 69 : return InnerProductImpls::InnerProductImpl<Lhs, Rhs>::apply(lhs, rhs); 70 : } 71 : 72 : } // namespace LinearSolver