Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <complex>
       7             : #include <cstddef>
       8             : 
       9             : #include "DataStructures/ComplexDataVector.hpp"
      10             : #include "DataStructures/DataVector.hpp"
      11             : #include "Utilities/Gsl.hpp"
      12             : 
      13             : namespace Spectral {
      14             : /// \ingroup SpectralGroup
      15             : /// Namespace for spin-weighted spherical harmonic utilities.
      16           1 : namespace Swsh {
      17             : 
      18             : /// \brief A set of labels for the possible representations of complex numbers
      19             : /// for the `ComplexDataView` and the computations performed in the
      20             : /// spin-weighted spherical harmonic transform library.
      21             : ///
      22             : /// \details The representation describes one of two behaviors:
      23             : ///  - `Interleaved`: The vectors of complex numbers will be represented by
      24             : ///  alternating doubles in memory. This causes both the real and imaginary part
      25             : ///  at a given gridpoint to be near one another, but successive real values
      26             : ///  farther. This is the native representation of complex data in the C++
      27             : ///  standard, and is the representation needed for Blaze math
      28             : ///  operations. Therefore, using this representation type in libsharp
      29             : ///  computations will cause operations which access only the real or imaginary
      30             : ///  parts individually to trace over larger memory regions. However, this
      31             : ///  representation will give rise to fewer copying operations to perform the
      32             : ///  libsharp operations.
      33             : ///
      34             : ///  - `RealsThenImags`: The vectors of complex numbers will primarily be
      35             : ///  represented by a pair of vectors of doubles, one for the real values and
      36             : ///  one for the imaginary values (the full computation cannot be performed
      37             : ///  exclusively in this representation, as it must return to a vector of
      38             : ///  `std::complex<double>` for Blaze math operations). This causes the
      39             : ///  successive real values for different gridpoints to be closer in memory, but
      40             : ///  the real and imaginary parts for a given gridpoint to be farther in
      41             : ///  memory. This is not the native representation for complex data in C++, so
      42             : ///  the data must be transformed between operations which use Blaze and the
      43             : ///  transform operations which use `RealsThenImags`. Therefore, using this
      44             : ///  representation in libsharp computations will cause operations which act on
      45             : ///  real or imaginary parts individually to have better memory locality (so
      46             : ///  likely improved cache performance, but such statements are highly
      47             : ///  hardware-dependent). However, this representation will give rise to more
      48             : ///  copying operations to perform the libsharp operations.
      49             : ///
      50             : /// \note The pair of representations is provided as a means to 'turn a dial' in
      51             : /// optimizations. It is unclear which of these representations will be
      52             : /// preferable, and it may well be the case that different representations are
      53             : /// better for different calculation types or different hardware. Therefore,
      54             : /// when optimizing code which uses libsharp, it is encouraged to profile the
      55             : /// cost of each representation for a computation and choose the one which
      56             : /// performs best.
      57           0 : enum class ComplexRepresentation { Interleaved, RealsThenImags };
      58             : 
      59             : namespace detail {
      60             : // A storage container for storing sub-vector references ("views") of a
      61             : // ComplexDataVector.
      62             : //
      63             : // This class takes as a template argument a ComplexRepresentation to use when
      64             : // returning pointers to the components of the complex data. The representation
      65             : // is either:
      66             : // - `ComplexRepresentation::Interleaved`, which indicates that the complex data
      67             : // is represented as a vector of `std::complex<double>`, and then
      68             : // ComplexDataView acts as a view (a pure reference to a subvector).
      69             : // - `ComplexRepresentation::RealsThenImags`, which indicates that the complex
      70             : // data is represented as a pair of vectors of `double`. This is no longer a
      71             : // strict view in the typical sense, as the memory layout is different from
      72             : // the data which it references. In order to minimize copies, the user must
      73             : // specify when edits to the ComplexDataView are complete by calling the
      74             : // `copy_back_to_source()` member function.
      75             : //
      76             : // Warning: For optimization reasons, mutations applied via member functions or
      77             : // manipulations of data pointed to by pointers returned by member functions
      78             : // *may or may not* be immediately applied to the source vector used to
      79             : // construct the ComplexDataView. Correct use of this class will first perform
      80             : // (potentially several) manipulations of the data via the ComplexDataView, then
      81             : // as a final step flush the data to the source vector via the member function
      82             : // `copy_back_to_source()`. At that point, the vector is guaranteed to be
      83             : // placed in the same state as the ComplexDataView. **The process of
      84             : // constructing a ComplexDataView, then mutating data in both the
      85             : // ComplexDataView and the source vector, then calling `copy_back_to_source()`
      86             : // is considered undefined behavior, and depends on the details of the memory
      87             : // layout chosen.**
      88             : template <ComplexRepresentation Representation>
      89             : class ComplexDataView {
      90             :  public:
      91             :   // The Representation determines the internal data representation
      92             :   static const ComplexRepresentation complex_representation = Representation;
      93             : 
      94             :   // The internal storage type
      95             :   using value_type = std::complex<double>;
      96             : 
      97             :   // Construct a view which starts at index `offset` of the supplied
      98             :   // vector, and extends for `size` elements
      99             :   ComplexDataView(gsl::not_null<ComplexDataVector*> vector, size_t size,
     100             :                   size_t offset = 0);
     101             : 
     102             :   // Construct a view which starts at pointer `start` and extends for
     103             :   // `size` elements. Need not be a part of a ComplexDataVector.
     104             :   ComplexDataView(std::complex<double>* start, size_t size);
     105             : 
     106             :   // For the lifetime of the data view, it points to the same portion of a
     107             :   // single vector. We disallow default move assignment, as it would exhibit
     108             :   // behavior that would contradict that. All assignment operations permitted by
     109             :   // the ComplexDataView are copying operations which act only on the data, not
     110             :   // the reference.
     111             :   ComplexDataView() = delete;
     112             :   ComplexDataView(const ComplexDataView&) = default;
     113             :   ComplexDataView(ComplexDataView&&) = default;
     114             :   ComplexDataView operator=(ComplexDataView&&) = delete;
     115             :   ~ComplexDataView() = default;
     116             : 
     117             :   // assign into the view the values of a same-sized ComplexDataVector
     118             :   ComplexDataView<Representation>& operator=(const ComplexDataVector& vector);
     119             : 
     120             :   // Assign into the view the values from a same-sized view
     121             :   ComplexDataView<Representation>& operator=(
     122             :       const ComplexDataView<Representation>& view);
     123             : 
     124             :   // Conjugate the data.
     125             :   void conjugate();
     126             : 
     127             :   // Assign into the real components of a view the values from a
     128             :   // provided `DataVector`
     129             :   ComplexDataView<Representation>& assign_real(const DataVector& vector);
     130             : 
     131             :   // Assign into the imaginary components of a view the values from a
     132             :   // provided `DataVector`
     133             :   ComplexDataView<Representation>& assign_imag(const DataVector& vector);
     134             : 
     135             :   // Gets the size of the view (the number of complex entries).
     136             :   size_t size() const { return size_; }
     137             : 
     138             :   // Gets the stride between successive real or successive imaginary components
     139             :   // in memory.
     140             :   static constexpr size_t stride() { return stride_; }
     141             : 
     142             :   // Gets the raw pointer to the start of the real data, which are separated
     143             :   // from one another by `stride()`.
     144             :   double* real_data();
     145             :   const double* real_data() const;
     146             : 
     147             :   // Gets the raw pointer to the start of the imaginary data, which are
     148             :   // separated from one another by `stride()`.
     149             :   double* imag_data();
     150             :   const double* imag_data() const;
     151             : 
     152             :   // A function for manually flushing the data back to the source. This class
     153             :   // minimizes copies wherever possible, so the manual control of bringing the
     154             :   // data back to the source vector only when necessary is important for
     155             :   // optimization.
     156             :   // Warning: Until this function is called, mutations applied via other member
     157             :   // functions or mutations applied via the pointers obtained from other member
     158             :   // functions are not guaranteed to be applied to the source vector.
     159             :   void copy_back_to_source();
     160             : 
     161             :  private:
     162             :   size_t size_;
     163             : 
     164             :   static const size_t stride_ =
     165             :       Representation == ComplexRepresentation::RealsThenImags ? 1 : 2;
     166             : 
     167             :   // In either case, we want to avoid unnecessary copies, so we track whether
     168             :   // the data has been copied from one representation to the other.
     169             :   bool real_slices_up_to_date_;
     170             :   // These two DataVectors are unused in the case of `Interleaved`
     171             :   // representation
     172             :   DataVector real_slice_;
     173             :   DataVector imag_slice_;
     174             : 
     175             :   double* data_real_;
     176             :   double* data_imag_;
     177             : 
     178             :   ComplexDataVector complex_view_;
     179             : };
     180             : 
     181             : }  // namespace detail
     182             : }  // namespace Swsh
     183             : }  // namespace Spectral