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
|
